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

calvin

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

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

ICH Q1E 统计分析方法

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

描述

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

分析方法：

假设：a1为时间（月份），a2为杂质含量

方法：回归分析

软件：R console 3.3.1

代码：如下

1. #Input stability data, Change:data in bracket
   
2. 
   
3. Stability_month<-c(0,6,18,24)
   
4. 
   
5. Stability_impurity1<-c(0.012,0.020,0.03,0.032)
   
6. 
   
7. Stability_impurity2<-c(0.008,0.012,0.023,0.035)
   
8. 
   
9. Stability_impurity3<-c(0.01,0.020,0.027,0.031)
   
10. 
    
11. Stability_content1<-c(0.995,0.99,0.991,0.978)
    
12. 
    
13. Stability_content2<-c(0.989,0.990,0.982,0.98)
    
14. 
    
15. Stability_content3<-c(0.981,0.985,0.976,0.975)
    
16. 
    
17. 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)
    
18. 
    
19. Stability_data<-stack(Stability_data,select=-Stability_month)
    
20. 
    
21. Stability_data$lot<-rep(1:3,each=4)
    
22. 
    
23. Stability_data$month=Stability_month
    
24. 
    
25. Stability_data$item=c(rep("impurity",12),rep("content",12))
    
26. 
    
27. Stability_data<-Stability_data[,-2]
    
28. 
    
29. names(Stability_data)=c("Result","lot","month","item")
    
30. 
    
31. #############factors combination analysis#############
    
32. 
    
33. #variance equivalence test
    
34. 
    
35. bartlett.test(Result~item,data=Stability_data)
    
36. 
    
37. #P value is 0.2431, cann't reject H0 hypothesis
    
38. 
    
39. #covariance test
    
40. 
    
41. summary(aov(Result~month*item,data=Stability_data))
    
42. 
    
43. #there's a significant slope, so the items should be analysis separately.
    
44. 
    
45. #############batch combination analysis#############
    
46. 
    
47. Stability_data_impurity<-Stability_data[Stability_data$item=="impurity",]
    
48. 
    
49. Stability_data_content<-Stability_data[Stability_data$item=="content",]
    
50. 
    
51. ###impurity analysis
    
52. 
    
53. #variance equivanlence test
    
54. 
    
55. bartlett.test(Result~lot,data=Stability_data_impurity)
    
56. 
    
57. #or
    
58. 
    
59. #library(car)
    
60. 
    
61. #ncvTest(fitmodel)
    
62. 
    
63. #P value is 0.8742, cann't reject H0 hypothesis
    
64. 
    
65. summary(aov(Result~lot*month,data=Stability_data_impurity))
    
66. 
    
67. #the 2-order interaction chi-square test result is P value=0.947, no significant difference of slope
    
68. 
    
69. #fit model
    
70. 
    
71. mode_impurity<-lm(Result~month,data=Stability_data_impurity)
    
72. 
    
73. #mode test
    
74. 
    
75. opar=par(no.readonly=T)
    
76. 
    
77. par(mfrow=c(2,2))
    
78. 
    
79. plot(mode_impurity)
    
80. 
    
81. #standard error test
    
82. 
    
83. #independence test
    
84. 
    
85. library(car)
    
86. 
    
87. durbinWatsonTest(mode_impurity)
    
88. 
    
89. #variance equivalence test
    
90. 
    
91. #method 1
    
92. 
    
93. ncvTest(mode_impurity)
    
94. 
    
95. #method 2
    
96. 
    
97. bartlett.test(Result~month,data=Stability_data_impurity)
    
98. 
    
99. #predict
    
100. 
     
101. Result<-predict(mode_impurity,data.frame(month=48),interval='confidence')[[3]]
     
102. 
     
103. Result
     
104. 
     
105. #visualization
     
106. 
     
107. month_interval=seq(0,48,by=0.5)
     
108. 
     
109. predict_month=data.frame(month=month_interval)
     
110. 
     
111. Line_predict<-predict(mode_impurity,predict_month,interval="confidence",level=0.95)
     
112. 
     
113. plot(Stability_data_impurity$month,Stability_data_impurity$Result,xlim=c(0,50),ylim=c(0,0.1))
     
114. 
     
115. abline(mode_impurity)
     
116. 
     
117. #the same as month_interval.
     
118. 
     
119. lines(month_interval,Line_predict[,3],col="red")
     
120. 
     
121. lines(month_interval,Line_predict[,2],col="red")
     
122. 
     
123. par(new=T)
     
124. 
     
125. plot(48,Result,xlim=c(0,50),ylim=c(0,0.1),xaxt="n")
     
126. 
     
127. segments(48,0,48,Result)
     
128. text(48,Result,paste("(48,",round(Result,4)*100,"%)"),pos=3)
     
129. #6% is lager than 5%, predict result at 36th months
     
130. predict(mode_impurity,data.frame(month=36),interval='confidence')[[3]]
     
131. #Reuslt is 4.8%
     
132. c("retest period is 36 months")

复制代码

模型诊断

预测图形

引自：Drug base （http://www.drugqm.com/?id=18）

本作品采用知识共享 署名-非商业性使用-相同方式共享 4.0 国际 许可协议进行许可。要查看该许可协议，可访问 http://creativecommons.org/licenses/by-nc-sa/4.0/ 或者写信到 Creative Commons, PO Box 1866, Mountain View, CA 94042, USA。

红色de蒲公英

想看具体分析过程

syhorchid

了解一下

birdgun

谢谢分享

一片红云

谢谢分享

xqliu

学习一下啦，谢谢提供分享。

chen5385973

谢谢分享

九三

高手呀，代码看不懂，惭愧！

DK1024

除过卧槽，牛逼之外我还能说点什么呢？

calvin

九三 发表于 2017-6-28 17:23
高手呀，代码看不懂，惭愧！

这个贴子的数据分析是基于GNU R Version 3.3.1的，其实对于药品生产企业质量管理，乃至绝大多数生产企业来说，用Minitab或JMP是最方便的，学习曲线较好。R语言的学习曲线太平坦，很难提高。不过对于数据分析自动化，可重复性报告来说，还是很有用的。FDA也用的。

我也看到你在论坛发了大量统计学的文章，还是蛮有深度的啊。以后一起探讨。

jxl2015

有些深奥，看不懂

叶子2005

谢谢分享，谢谢分享

lqthjm

谢谢分享！

李牧

看了一头雾水呀

微信e4fac340

谢谢分享，多一点说明就好了

雪少

卧槽，太厉害了