欢迎您注册蒲公英
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
本帖最后由 calvin 于 2017-6-25 21:41 编辑
药品稳定性研究从研发阶段的长期、加速,以及生产阶段的承诺稳定性考察,持续稳定性考察。都积累了大量数据,对于这些数据,有什么用,我们又该怎么分析呢?
经过思考,打消了把代码写成自定义函数的想法~下面是原始代码。
ICH Q1E 统计分析方法
下载链接:http://pan.baidu.com/s/1jGB6tuu 密码:4g35
描述
分析方法: 假设:a1为时间(月份),a2为杂质含量 方法:回归分析 软件:R console 3.3.1 代码:如下 - #Input stability data, Change:data in bracket
- Stability_month<-c(0,6,18,24)
- Stability_impurity1<-c(0.012,0.020,0.03,0.032)
- Stability_impurity2<-c(0.008,0.012,0.023,0.035)
- Stability_impurity3<-c(0.01,0.020,0.027,0.031)
- Stability_content1<-c(0.995,0.99,0.991,0.978)
- Stability_content2<-c(0.989,0.990,0.982,0.98)
- Stability_content3<-c(0.981,0.985,0.976,0.975)
- Stability_data=data.frame(Stability_month=Stability_month,Stability_impurity1=Stability_impurity1,Stability_impurity2=Stability_impurity2,Stability_impurity3=Stability_impurity3,Stability_content1=Stability_content1,Stability_content2=Stability_content2,Stability_content3=Stability_content3)
- Stability_data<-stack(Stability_data,select=-Stability_month)
- Stability_data$lot<-rep(1:3,each=4)
- Stability_data$month=Stability_month
- Stability_data$item=c(rep("impurity",12),rep("content",12))
- Stability_data<-Stability_data[,-2]
- names(Stability_data)=c("Result","lot","month","item")
- #############factors combination analysis#############
- #variance equivalence test
- bartlett.test(Result~item,data=Stability_data)
- #P value is 0.2431, cann't reject H0 hypothesis
- #covariance test
- summary(aov(Result~month*item,data=Stability_data))
- #there's a significant slope, so the items should be analysis separately.
- #############batch combination analysis#############
- Stability_data_impurity<-Stability_data[Stability_data$item=="impurity",]
- Stability_data_content<-Stability_data[Stability_data$item=="content",]
- ###impurity analysis
- #variance equivanlence test
- bartlett.test(Result~lot,data=Stability_data_impurity)
- #or
- #library(car)
- #ncvTest(fitmodel)
- #P value is 0.8742, cann't reject H0 hypothesis
- summary(aov(Result~lot*month,data=Stability_data_impurity))
- #the 2-order interaction chi-square test result is P value=0.947, no significant difference of slope
- #fit model
- mode_impurity<-lm(Result~month,data=Stability_data_impurity)
- #mode test
- opar=par(no.readonly=T)
- par(mfrow=c(2,2))
- plot(mode_impurity)
- #standard error test
- #independence test
- library(car)
- durbinWatsonTest(mode_impurity)
- #variance equivalence test
- #method 1
- ncvTest(mode_impurity)
- #method 2
- bartlett.test(Result~month,data=Stability_data_impurity)
- #predict
- Result<-predict(mode_impurity,data.frame(month=48),interval='confidence')[[3]]
- Result
- #visualization
- month_interval=seq(0,48,by=0.5)
- predict_month=data.frame(month=month_interval)
- Line_predict<-predict(mode_impurity,predict_month,interval="confidence",level=0.95)
- plot(Stability_data_impurity$month,Stability_data_impurity$Result,xlim=c(0,50),ylim=c(0,0.1))
- abline(mode_impurity)
- #the same as month_interval.
- lines(month_interval,Line_predict[,3],col="red")
- lines(month_interval,Line_predict[,2],col="red")
- par(new=T)
- plot(48,Result,xlim=c(0,50),ylim=c(0,0.1),xaxt="n")
- segments(48,0,48,Result)
- text(48,Result,paste("(48,",round(Result,4)*100,"%)"),pos=3)
- #6% is lager than 5%, predict result at 36th months
- predict(mode_impurity,data.frame(month=36),interval='confidence')[[3]]
- #Reuslt is 4.8%
- c("retest period is 36 months")
复制代码模型诊断
预测图形
引自:Drug base (http://www.drugqm.com/?id=18)
|