求大神帮忙用R语言自己写一个GLM广义线性模型的函数

Python033

求大神帮忙用R语言自己写一个GLM广义线性模型的函数,第1张

请教如何实现广义线性模型GLM作图

1、广义线性模型GLM很简单,举个例子,药物的疗效和服用药物的剂量有关。这个相关性可能是多种多样的,可能是简单线性关系(发烧时吃一片药退烧0.1度,两片药退烧0.2度,以此类推;这种情况就是一般线性模型),也可能是比较复杂的其他关系,如指数关系(一片药退烧0.1度,两片药退烧0.4度),对数关系等等。这些复杂的关系一般都可以通过一系列数学变换变成线性关系,以此统称为广义线性模型。

2、广义线性混合模型GLMM比较复杂,GLM要求观测值误差是随机的,而GLMM则要求误差值并非随机,而是呈一定分布的。举个例子,我们认为疗效可能与服药时间相关,但是这个相关并不是简简单单的疗效随着服药时间的变化而改变。更可能的是疗效的随机波动的程度与服药时间有关。比如说,在早上10:00的时候,所有人基本上都处于半饱状态,此时吃药,相同剂量药物效果都差不多。但在中午的时候,有的人还没吃饭, 有的人吃过饭了,有的人喝了酒,结果酒精和药物起了反应,有的人喝了醋,醋又和药物起了另一种反应。显然,中午吃药会导致药物疗效的随机误差非常大。这种疗效的随机误差(而非疗效本身)随着时间的变化而变化,并呈一定分布的情况,必须用广义线性混合模型了。

我是用的pscl包,zeroinfl()函数零膨胀负二项模型(ZINB)mod <- zeroinfl(ReportedNumber~ A+B+C+D+E | F+G+H+I, data = zinb, dist = "negbin", EM = TRUE)ZINB模型由点模型和零膨胀模型两部分结合而成,ABCDE是点模型内变量, 影响因变量发生次数的多少,FGHI是零膨胀模型内变量,决定因变量是否能够发生(为0还是非0)。 http://www.ats.ucla.edu/stat/r/dae/zinbreg.htm 这个网站里讲的很清楚

我们之前提到的线性回归是利用X来预测Y,Y是连续型的数值变量。但有时候Y并不是连续型的变量,而是一种 离散型的变量 。或者说,更准确来说,是一种 定类数据 。举个《统计学习导论》书上的例子,假设现在要通过一个急救室病人的症状来预测其患病情况。我们有三种可能的诊断:中风,服药过量,癫痫发作。我们分别用数字来表示这种诊断

这里的诊断结果就是Y,症状就是X。我们也可以用前面线性回归来做,但这样就是默认其实是一个有序的输出。但实际上,中风,服药过量,癫痫发作这三类数据虽然我们用数字来代表,但其实并不是有顺序之分的,所以用线性回归并不适合。所以,我们就可以考虑用logistic回归来解决这种 分类问题

在做logistic回归的时候,我们也会将我们最后的二元结果用数字来表示,一个代表1,一个代表0。我们最后预测能得到的是 y=1 的概率!

我们来看下logistic的模型

其中 p 代表 y = 1 的概率 ,x 代表了不同的自变量, 表示了误差项。与线性回归模型对比,等式右边完全相同,实际上逻辑回归模型也是广义上的线性模型。而等式的左边形式更复杂了,引入了一些非线性的变换。大家可以看到,我们这里等式左边变成了一个对数, 我们常称为对数发生比(log-odd)或分对数(logit) 。对数里面的 就是 发生比(odd) ,取值范围可以从0到正无穷。

然后我们估计回归参数的话,就变成

有时候,也会将logistic的模型写成(我不太喜欢这么写,但有时候看助教答案是这么写的)

下面是logistic回归的一些tip

用两个作业里面的题目举个例子:

数据文件“Drivers.csv”为对45名司机的调查结果,其中四个变量的含义为:

1.请在R语言中调用logistic回归函数,计算视力状况、年龄、驾车教育与是否发生事故的logistic回归模型,并以“odds=……”的形式写出回归公式。

R语言里面调用logistic回归函数是用 glm 函数。glm函数其实是拟合广义线性模型的函数,logistic回归只是其中一种。所以你要加上family=binomial,代表逻辑斯蒂回归

然后我们的logistic回归模型就是

但题目要求的是odds的形式,那么我们再改写下

2.指出(1)得到的模型中哪些因素对是否发生事故有显著性影响。如果存在对是否发生事故没有显著性影响的因素,请去除这些因素后重新计算logistic回归模型,并以“p=……”的形式写出回归公式

我们根据前面结果的p-value,就发现视力是由显著性影响的,其他因素没有显著性影响。

然后我们去掉这些因素,只留下视力这个因素,再次拟合一个logistic回归模型

这回是以p的形式

我们还可以比较下,这两个方程在统计学上是不是有差异的,就是跟线性回归一样用anova函数

发现两个模型是没有差异的。

3.A是一名参加过驾车教育,但视力有问题的50岁老司机;B是一名没有参加过驾车教育,但视力良好的20岁新手。现在A、B都想在某保险公司投保,但按公司规定,被保险人必须满足“明年出事故的概率不高于40%”的条件才能予以承保。请预测A、B两者明年出事故的概率,并告诉保险公司谁可以投保。

这里就是用我们构建的模型去预测结果。但有一个问题就是拿哪个模型去预测呢。 答案里面是用了精简以后的(去掉了不显著的变量)模型。

这题的数据比较多,且比较麻烦。我就不放数据了,写出来只是为了讲下参考变量和训练测试集等几个概念。

这题的大致目的就是用10个自变量去做出诊断,肿瘤是否是良性的(M = malignant(恶性的), B = benign(良性的))。总的来说,数据集有357个良性肿瘤,212个恶性肿瘤,即共有569个数据点。

1. Use all mean features to construct a logistic regression model

因为原始数据集是包括了mean,var等等。这里只要求用mean部分的数据。所以我们就先提取了mean那部分数据集,然后还是用glm。

2. Then try to reduce the number of features from your last model, construct another regression model,and you will need to write down the equation of your logistic regression model(Tips: Logit P = α+β1X1+β2X2+..+βpXp)

我们可以把显著性的变量挑出来,再次构建一个新的logistic回归模型

回归模型为

3. Use proper test to test the difference between two models

用anova就可以了

4. You may split the data properly, use part of them to train your regression model and use another part to make predictions. Lastly, you may try to calculate the accuracy of your model.(Tips: To split the data, you can use the first 398 rows as training data, use the last 171 rows as prediction data.The predict function return a value between 0 and 1, 0.~0.5 belong to the first class, and 0.5~1 belong to second class in binary classification problems)

这样一套下来,大家可能会感觉有些奇怪。在题目1的时候,1代表出事故,0代表没出事故。然后 里面的p代表的是y=1的概率,即出事故的概率。这一切都很顺理成章。但是在题目2的时候,M代表恶性的,B代表良性的。那为啥 里面的p代表的就是M呢。或者说为啥M代表的就是1,B代表的就是0呢。

事实上,对于二元变量,glm会确定你的 响应变量 里面谁是 reference level 。或者说,会确定谁是那个 0 。那么,glm是怎么确定的呢。

通过上面的讲述,我们就发现了B是reference level,即是0。而M是对立的那个level,即是1。而我们logistic输出的是y=1的概率,即y=恶性肿瘤的概率。然后也刚好对应的我们前面的

只有y=M的概率足够大,才定义为M。这里我们设定的 足够大 是0.5。你也可以认为大于0.9才算出是M,这样结果就会不太一样。

我们也可以把我们的响应变量直接变成0和1

那么,我们该如果更改我们的reference level呢,一种方法我觉得应该可以是像上面那样,响应变量直接变成0和1,非常直观,但个人觉得不太符合逻辑……。另一种就是下面那样

参考文章

https://stats.stackexchange.com/questions/207427/confused-with-the-reference-level-in-logistic-regression-in-r

https://stackoverflow.com/questions/17772775/change-reference-group-using-glm-with-binomial-family

https://stackoverflow.com/questions/23282048/logistic-regression-defining-reference-level-in-r