R语言-医药研究 logistic regression

科学技术献祖国 5 医药研究 logistic regression 介绍R 语言的应用实际上都是在介绍算法，算法的本质属于脑筋急转弯，所以算法是不能申请组专利保护的。然而算法一旦固化到电子产品当中，产生了可以摸得到的商品,这种商品就可以注册版权，受知识产权的保护。同样的新发现的药物也是可以申请专利受法律保护，这就是现在药物发掘的工作如火如荼的原因。那些制药公司都希望尽可能多的申请药物专利从中谋利。在这里我们将简单介绍新的药品是如何发现的。 药品实验都是在生物活体上进行的，使用的实验对象可以是：小白鼠，狗，或者是其它的单细胞，多细胞生物。关键的一点是使用的实验对象必须是完全相同的，实验结果才能有可比性。举一个例子，双胞胎个体就是完全相同的生物个体，因为它们是从单一的受精卵分裂出来的。但实验中我们有时需要用到上百个完全相同的实验个体，通过使用双胞胎不可能得到足够的实验对象。另外一种办法，可以通过选种的方法来生成众多的完全相同的实验个体。例如说我们可以有一百个小白鼠，它们每一个都有着完全相同的基因组成。因为这些小白鼠是通过选种培育出来的，它们都有着共同的单一祖先。在下面的实验中， 我们使用200个这样具有相同基因组成的小白鼠。 介绍疼痛实验，一般是使用一束强光照射到小白鼠的尾巴上。如果小白鼠的尾巴突然间抽动就说明小白鼠有了疼痛的感觉，这种实验是用来研究哺乳类动物对疼痛承受能力的方法。 吗啡是一种止疼药物，使用吗啡后哺乳类动物都会迟钝对疼痛的感觉。研究人员新近发现了另外一种药物被称作代尔9号也可以用作麻醉剂。研究人员想用吗啡与代尔9号相比较来明白哪种药物更有效。研究人员也想知道当这两种药物同时被使用的时候，它们的药性可不可以相互促进。发现这种相互促进的药性是十分重要的，因为单一使用某种药物可能会毒性非常的大。如果有另外一种药物能使它们之间互相促进，就可以在两种药物都使用非常小的剂量的状况下达到比较好的医治效果。这种药物互相作用在研究杀虫剂的课题中也经常被用到。经常是单一的药物不可能对害虫进行有效的杀灭。我们就需要寻找其它的药物与其互相作用来提高杀灭的效果。所以世面上的药物都是复合了几种化学成分，很少有单一的药物。 下面我们来逐一介绍这个实验所需要的R 命令 flick<-read.table("~mcdougal/zBook/zData/c71.html", header=TRUE, skip=30, nrows=20) flick 上面的命令生成了一个叫做flick 数据对象，使用read.table 命令对数据所在的因特网链接进行扫描。因为因特网上的数据中有很多不需要的部分，所以要使用skip 的选项来跳过前30行。然后使用nrows 的选项只读取其中的20 行。在下面的截图中可以看到数据输出后的样子 这个数据有20 行7列每一行代表10 个小白鼠，第一列代表的是行数从1 到20，第二列代表的是使用的药物，A 代表吗啡B 代表代尔9号AB 代表吗啡和代尔9号同时使用O 代表不使用任何药物。第三列代表的是吗啡的用量单位是mg/kg。 第四列代表的是代尔9号的用量单位是mg/kg. 第五列代表的是小白鼠尾巴没有反应的个数。第六列代表的是每组所使用小白鼠的个数。第七列代表的是每一组小白鼠尾巴没有反应的百分数。 每一组小白鼠都要根据上面的药物剂量注射吗啡和代尔9号，然后对所有的小白鼠的尾巴进行强光照射。如果小白鼠已经被麻醉了它的尾巴就不会有抽动的反应。我们用这种方法来测试这两种麻醉剂的药效。 比方说数据的第一行代表的是第一组小白鼠，使用的药物符号为O 表示没有使用吗啡也没有使用代尔9号。所以在吗啡和代尔9号用量上都是0.0mg/kg，也就是注射的生理盐水。第一行的第五列数据值为0，说明在没有使用吗啡和代尔9号的时候所有的小白鼠都对热的刺激有了反应。这一组小白鼠有10 个。最后对热没有反应的小白鼠的百分数为0. 再比方说数据的第二行代表了第二组小白鼠，也有10只。这组小白鼠只注射了1mg/kg的吗啡，它们全部对热刺激有了抽动的反应，所以没有反应的百分数为0。 flick.glm<-glm(prop~morphine*del9,w=rep(10,20), data=flick,family=binomial) summary(flick.glm) qqnorm(resid(flick.glm),main="model 1") vcov(flick.glm) 在上面的命令中我们使用glm命令来建立逻辑回归模型。也就是用吗啡和代尔9号的用量来预测小白鼠被麻醉的百分数。其中family 选项设置为binomial 说明使用逻辑回归建模而不是普通的回归。因为小白鼠被麻醉的数据是百分数，对百分数的回归模型必须是逻辑回归。具体逻辑回归的细节就不在这里进行介绍，读者可以自行学习。Summary 命令帮助我们输出模型的p-value 和相关参数。Qqnorm命令生成模型误差的qq图。关于qq图的解释请读者参考 R 语言样品比较应用举例。上面vcov命令输出的是模型参数之间方差协方差矩阵。下面是输出结果截图 模型的系数有4个：第一个是截距 值为-3.2832.第二个参数也就是 是吗啡的斜率代表吗啡的药性值为0.7604 p-value 为4.93e-11 几乎为0，就是说吗啡的药效是显著的。同样代尔9号的斜率 值为0.9703 p-value 也是非常的小说明代尔9号的药效也是显著的。 第四个参数 代表的是药物互相作用的斜率值为0.2588 p-value 为0.0548 稍微的大于0.05，说明这两个药物之间的互相作用并不是十分的明显。但是这只是我们的第一个模型，下面我们还要介绍一个改进过的模型。 临床医生在使用麻醉剂的时候会用到50%有效剂量（ED50）这个概念。这个概念的意思是，用多少剂量的药物可以使50%的病人被麻醉，或者说多大剂量的药物可以刚好对50%的病人有效。我们用ED50 代表这个剂量，计算这个剂量的公式如下： ED50= 从上面知道， 的数值为-3.2832 ， 的数值为0.7604 所以吗啡的ED50 就是：3.2832/0.7604 等于4.317728 mg/kg 。下面是如何使用R 语言来进行ED50 的计算。 mor.ED50=-(-3.2832/0.7604) del.ED50=-(-3.2832/0.9703) mor.ED75=(log(.75/(1-.75))-(-3.2832))/0.7604 del.ED75=(log(.75/(1-.75))-(-3.2832))/0.9703 mor.ED50 #吗啡ED50 del.ED50 #代尔9号ED50 mor.ED75 #吗啡ED75 del.ED75 #代尔9号ED75 我们同样可以计算ED75 也就是使75%的病人有效的剂量，计算公式为： 来看输出结果截图： ED75 高于ED50 说明要使更多的病人被麻醉就需要使用更高的药物剂量。从上面的结果可以看出代尔9号的药效比较快，需要用的剂量比较低。 上面介绍ED50 的计算都是吗啡和代尔9号独立使用时的ED50。下面来介绍但这两种药同时使用时ED50 的计算。我们需要事先固定代尔9号的用量然后计算相应所需的吗啡用量。公式如下： Dose1= 从上面的模型以及估计出了所有 参数的具体数值。只需要将这些参数值插入上面的公式中，我们把代尔9号的用量固定在1.25mg/kg 。因为通过以前的一些实验中我们知道人体对代尔9号的最大承受能力就在1.25mg/kg左右。下面就是所使用的R 语句和输出： dose.mor=(-(-3.2832)-0.9703*1.25)/(0.7604+0.2588*2.35) dose.mor 可以看到将代尔9号的用量固定在1.25 mg/kg，同时只需要用1.512754 mg/kg 的吗啡即可达到50% 病人进行麻醉的效果。 下面来介绍改建了的模型。在药物的使用过程中经常是一开始使用较小的剂量就可以有很明显的效果，量药升高但药效的增加但却不明显，也就是说用药的剂量与药效之间的关系不是线性的。所以我们要使用LOG 变形来改进我们的模型，同时出现的问题是LOG（ 0） 是不被定义的。比方说在我们的数据中，当吗啡不被使用的时候吗啡的剂量就是0，就无法用LOG 来对吗啡的剂量进行变形。所以我们需要将所有吗啡的剂量在数值上都加1使之都大于0 就可以进行LOG 变形了。下面是具体的R 命令。 #log+1 model log.morphine=log(flick$morphine+1) #对吗啡用量LOG+1 变形 log.del9=log(flick$del9+1) #对代尔9号用量LOG+1 变形 prop=flick$prop #提出小白鼠被麻痹的百分比 flick3<-data.frame(cbind(log.morphine,log.del9,prop)) #生成建模所需数据矩阵 flick.glm3<-glm(prop~log.morphine*log.del9,w=rep(10,20), data=flick3,family=binomial) #进行LOG+1变形后的逻辑回归模型 summary(flick.glm3) #输出模型参数，所有的 qqnorm(resid(flick.glm3),main="model 3") #使用QQ 图检验正态分布 值得注意的是在新的模型中，这两种药物的互相作用产生了0.043 的p-value 小于0.05 也就是说我们有相当的证据来证明这两种药物之间存在着显著的互相作用关系。因为我们使用了更符合药物研究特性的模型，我们才能有机会检测到这种互相关系。下面我们来继续使用模型参数，计算ED50 和ED75， 代码如下： mor.ED50.3=exp(-(-6.6772/4.1324))-1 继续阅读