本文是《Quality by Design for ANDAs: An Example for Immediate-Release Dosage Forms》第一个处方的R语言解决方案。
第一个处方研究评估原料药粒径分布、MCC/Lactose比例、崩解剂用量对制剂CQAs的影响。
第二处方研究用于理解颗粒外加硬脂酸镁和滑石粉对片剂质量和可生产性的影响。处方研究在实验规模进行（1.0kg，5000片）。
处方研究1的目的是选择MCC/Lactose的比例以及崩解剂的用量，并理解这些变量是否与原粒药粒径分布有交互作用。这个研究也用来确立处方稳健性。这些处方因素对响应变量的影响的带3个中心点的23全因子DOE列于表21。
崩解剂加到颗粒内，用量范围为1-5%。MCC/Lactose的比例为1：2，1：1，2：1。原料药的重量占10%。滑石粉的用量占2.5%。硬脂酸镁的用量占1%。片重为200.0mg。
用三种滚压力滚压物料。调整压片力（允许的范围是11-13.0kP）以研究片剂硬度为12.0kP时的溶出。选择12.0kP的片剂硬度以研究处方变量对溶出的影响因为更高的硬度是溶出的最差条件。如果固定压片力则可以混杂片剂硬度的影响。
表21列出了考察的因子和响应。
表22列出了溶出，含量均一性，粉末流动系数和10KN压片力的片剂硬度。
library(FrF2)
study1<-FrF2(nruns=8,nfactors=3, ncenter=3,factor.names=list(A=c("-1","+1"),B=c("-1","+1"),C=c("-1","+1")), replications=1,randomize=FALSE)
> study1
y1<-c(99.5,77.0,98.7,86.0,99.0,76.0,99.0,84.0,91.0,89.4,92.0)
y3<-c(4.1, 2.9, 4.0, 3.2, 5.0, 3.8, 5.1, 4.0, 4.0, 3.9, 4.1)
y5<-c(6.16,8.46,6.09,8.46,4.97,7.56,4.77,7.25,6.62,6.66,6.46)
y7<-c(9.1,8.3, 9.1, 8.6, 13.5,12.5,12.9,13.2,10.6,10.9,11.3)
study1 <- add.response( study1, y1)
study1 <- add.response( study1, y3)
study1 <- add.response( study1, y5)
study1 <- add.response( study1, y7)
mod1 <- lm( y1 ~A*B, data = study1)
anova(mod1)
将Residuals分解为弯曲性、纯误差、Lack of Fit:
SS Pure Error=((91.0)^2+(92.0)^2+(89.4)^2-(91.0+92.0+89.4)^2/3)=3.44
纯误差的自由度为3-1=2
MS Pure Error =3.44/2=1.72
SS Curvature=8*3*((84.0+98.7+77.0+99.0+99.0+99.5+86.0+76.0)/8-(91.0+92.0+89.4)/3)^2/(8+3)= 1.77
弯曲性的自由度为1
MS Curvature=1.77/1=1.77
F Curvature=MS Curvature/MS Residuals = 1.74
SS Lack of fit=7.88-3.44-1.77=2.67
Lack of fit的自由度为4
MS lack of fit= SS Lack of fit/4=0.67
F lack of fit=MS lack of fit/MS Pure Error=0.39
Pr(>F) Curvature=1-pf(1.74,1,6)=0.2352431
Pr(>F) lack of fit=1-pf(0.39,4,2)=0.8079788
最终的方差分析表如下：
Response: y1
Df Sum Sq Mean Sq F value Pr(>F)
A 1 669.78 669.78 657.72 2.316471e-07 ***
B 1 32.81 32.81 32.22 0.001287723 **
A:B 1 39.60 39.60 38.89 0.00078691 ***
Curvature 1 1.77 1.77 1.74 0.2352431
Residuals 6 6.11 1.02 - -
lack of fit 4 2.67 0.67 0.39 0.8079788
Pure Error 2 3.44 1.72
表23是溶出度的ANOVA分析结果。
如表23所示，溶出度的弯曲效应并不显著。用所有的数据（包括中心点）拟合模型。
library(FrF2)
study1<-FrF2(nruns=8,nfactors=3, ncenter=3,factor.names=list(A=c("-1","+1"),B=c("-1","+1"),C=c("-1","+1")), replications=1,randomize=FALSE)
study1
y1<-c(99.5,77.0,98.7,86.0,99.0,76.0,99.0,84.0,91.0,89.4,92.0)
y3<-c(4.1, 2.9, 4.0, 3.2, 5.0, 3.8, 5.1, 4.0, 4.0, 3.9, 4.1)
y5<-c(6.16,8.46,6.09,8.46,4.97,7.56,4.77,7.25,6.62,6.66,6.46)
y7<-c(9.1,8.3, 9.1, 8.6, 13.5,12.5,12.9,13.2,10.6,10.9,11.3)
study1 <- add.response( study1, y1)
study1 <- add.response( study1, y3)
study1 <- add.response( study1, y5)
study1 <- add.response( study1, y7)
mod1 <- lm( y1 ~A*B, data = study1)
anova(mod1)
> anova(mod1)
Analysis of Variance Table
Response: y1
Df Sum Sq Mean Sq F value Pr(>F)
A 1 669.78 669.78 595.188 4.951e-08 ***
B 1 32.81 32.81 29.152 0.0010093 **
A:B 1 39.60 39.60 35.194 0.0005802 ***
Residuals 7 7.88 1.13
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
将Residuals分解为纯误差、Lack of Fit:
SS Pure Error=((91.0)^2+(92.0)^2+(89.4)^2-(91.0+92.0+89.4)^2/3)=3.44
纯误差的自由度为3-1=2
MS Pure Error =3.44/2=1.72
SS Lack of fit=7.88-3.44=4.44
Lack of fit的自由度为7-2=5
MS lack of fit= SS Lack of fit/5=0.888
F lack of fit=MS lack of fit/MS Pure Error=0.516
Pr(>F) lack of fit=1-pf(0.516,5,2)=0.7618
library(daewr)
fullnormal(coef(mod1)[-1], alpha=.025)
effects <-coef(mod1)
effects <-effects[2:4]
effects <-effects[ !is.na(effects) ]
library(daewr)
halfnorm(effects, names(effects), alpha=.25)
library(BsMD)
LenthPlot(mod1, main = "Lenth Plot of Effects")
A.num <-study1$A
levels(A.num) <- c(10,30)
B.num <- study1$B
levels(B.num) <- c(1,5)
A.num <- as.numeric(as.character(A.num))
B.num <- as.numeric(as.character(B.num))
mod1 <- lm( y1 ~A.num*B.num, data = study1)
library(rsm)
contour(mod1, ~ A.num+B.num)
persp(mod1, ~ A.num+B.num, zlab="y1", contours=list(z="bottom"))
如图10的半正态图和表24的未调整模型的ANOVA结果所示，影响溶出度的显著因子为A（原料药粒径）,B（崩解剂用量），AB（原料粒径与崩解剂用量交互作用）。
图11展示原粒药粒径与崩解剂用量对溶出度的影响。
A <- factor(c(1,-1, 1,-1,-1,-1,1, 1))
B <- factor(c(1, 1,-1,-1, 1,-1,1,-1))
C <- factor(c(1,-1,-1, 1, 1,-1,1, 1))
y1<-c(84.0, 98.7, 77.0, 99.0, 99.0, 99.5, 86.0,76.0)
y3<-c(4.0, 4.0, 2.9, 5.0, 5.1, 4.1, 3.2, 3.8 )
y5<-c(7.25, 6.09, 8.46, 4.97, 4.77, 6.16, 8.46,7.56)
y7<-c(13.2, 9.1, 8.3, 13.5, 8.3, 9.1, 8.6, 12.5)
study1<- data.frame (A=A,B=B,C=C,y1=y1, y3=y3, y5=y5, y7=y7)
write.csv( plan, file = "E://data//study1.csv", row.names = FALSE )
study1 <- read.csv("E://data//study1.csv")
mod1 <- lm( y1 ~A*B, data = study1)
anova(mod1)
> anova(mod1)
Analysis of Variance Table
Response: y1
Df Sum Sq Mean Sq F value Pr(>F)
A 1 669.78 669.78 1003.416 5.92e-06 ***
B 1 32.80 32.80 49.146 0.002180 **
A:B 1 39.60 39.60 59.333 0.001529 **
Residuals 4 2.67 0.67
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
with(study1, (interaction.plot(A,B,y1, type = "b", pch = c(18,24,22), leg.bty = "o", main = "Interaction Plot of A and B", xlab = "A",ylab = "y1")))
#以下将因子变量转为自然变量
A.num <-study1$A
levels(A.num) <- c(10,30)
B.num <- study1$B
levels(B.num) <- c(1,5)
A.num <- as.numeric(as.character(A.num))
B.num <- as.numeric(as.character(B.num))
mod1 <- lm( y1 ~A.num*B.num, data = study1)
contour(mod1, ~ A.num + B.num)
persp(mod1, ~ A.num + B.num, zlab=" y1", contours=list(z="bottom"))
影响崩解时间的显著因子
崩解剂是唯一影响片剂崩解时间的显著因子。但是所有的批次都在小于4分钟内快速的崩解。
影响片剂含量的显著因子
所有的批次含量均在标准范围内(95.0-105.0% w/w) ，没有因子对含量有显著影响。