26
2015
10

药品稳定性实验数据分析及可视化

药品稳定性研究从研发阶段的长期、加速,以及生产阶段的承诺稳定性考察,持续稳定性考察。都积累了大量数据,对于这些数据,有什么用,我们又该怎么分析呢?

 

经过思考,打消了把代码写成自定义函数的想法~下面是原始代码。

 

ICH Q1E 统计分析方法

 

下载链接:http://pan.baidu.com/s/1jGB6tuu 密码:4g35

 

描述

  1. 研发阶段有效期制定

    研发阶段,稳定性实验数据除了根据药典及一些指导原则规定的,外括法来制定有效期外,还应对其进行统计分析,以获得更加合理的有效期,并且,随着数据的增加,可以作为延长有效期的依据。 

  2. 生产阶段的质量评估

    生产阶段的稳定性数据可以用来评估产品质量稳定性,以评估在有效期内药品是合格的。

 

分析方法:

假设:a1为时间(月份),a2为杂质含量

方法:回归分析

软件:R console

代码:如下

 

#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")

 

 

模型诊断

 

222.jpeg

 

预测图型

 

regression.jpeg

111.jpeg

 

« 上一篇 下一篇 »

评论列表:

1.smile  2016-01-18 12:00:44 回复该评论
学习,谢谢Jonie分享

发表评论:

◎欢迎参与讨论,请在这里发表您的看法、交流您的观点。