vlambda博客
学习文章列表

R语言进阶之广义线性回归

广义线性回归是一类常用的统计模型,在各个领域都有着广泛的应用。今天我会以逻辑回归和泊松回归为例,讲解如何在R语言中建立广义线性模型。


在R语言中我们通常使用glm()函数来构建广义线性模型,glm实际上是generalized linear model(广义线性模型)的首字母缩写,它的具体形式如下所示:


glm(formula, family=familytype(link=linkfunction), data=)# formula就是我们的模型形式,family是我们指定的具体回归类型(见下表)


Family
Default Link Function
binomial
(link = "logit")
gaussian
(link = "identity")
Gamma
(link = "inverse")
inverse.gaussian
(link = "1/mu^2")
poisson
(link = "log")
quasi
(link = "identity", variance =  "constant")
quasibinomial
(link = "logit")
quasipoisson
(link = "log")


你可以通过 help(glm)来查看其它的模型选项, 使用 help(family)来查看每一族的连接函数。在这里我主要和大家讲解一下逻辑(logistic)回归和泊松(poisson)回归这两个模型。


第一部分 逻辑回归


逻辑回归主要应用于因变量(y)是二分类变量而自变量(x)是连续型变量的情形,当然这里的自变量和因变量也可以都是分类变量。由于逻辑回归本身的假设条件并没有那么严格,因此它的应用范围比判别分析要更广。关于判别分析的知识,我会在后续内容中和大家详细介绍。


这里我们使用鸢尾花(iris)数据集,将setosa这一类去掉后鸢尾花的种类(Species)就是一个二分类变量,将virginica设置为0,versicolor设置为1,使用花瓣和花萼数据来预测鸢尾花的种类。


# 逻辑回归mydata <- iris[which(iris$Species!= "setosa"),] # 去除setosa这类mydata$type[which(mydata$Species== "virginica")] <- 0 # virginica编号为0mydata$type[which(mydata$Species== "versicolor")] <- 1 # versicolor编号为1 attach(mydata) # 固定数据集fit <-glm(type~Sepal.Length+Sepal.Width+Petal.Length+Petal.Width,data=mydata,family=binomial())# 指定分布为二项分布summary(fit) # 输出拟合结果confint(fit) # 输出95%置信区间exp(coef(fit)) # 取拟合系数的自然指数exp(confint(fit)) # 取95%置信区间的自然指数predict(fit, type="response") # 输出拟合值residuals(fit, type="deviance") # 输出残差



从输出结果来看,花瓣长度是可以较好区分这两类鸢尾花的,但是这个模型是原始和粗糙的,我们应该通过回归诊断的方式来修正此模型,使之更加精确,关于回归诊断请参见,这里就不赘述。


当然我们也可以用anova(fit1,fit2,test="Chisq")来比较模型的优劣,这个在入门阶段也已经介绍过了,不明白的可以参考往期内容。

 

第二部分 泊松回归


泊松回归主要用于因变量(y)是计数资料而自变量(x)是连续型变量的时候,当然自变量(x)也可以是分类变量。


这里我先和大家介绍一下数据的信息,这个数据主要包括三部分信息:treatment代表对患者采取的治疗措施,分成1、2、3三类,1代表被认可的有效药,2代表新药A,3是指新药B;outcome是指患者治疗之后的结局,同样可分成1、2、3三类,1代表病情好转,2代表病情迁延不愈(没恶化),3代表病情恶化;counts是指采取不同治疗措施的不同结局的患者个数,是一个计数资料。注意这里不使用安慰剂作为空白对照的原因主要是考虑到伦理学问题,原则上要使患者利益最大化。


# 泊松回归# counts是计数值# outcome是指患者治疗后可能的结局# treatment是指对患者采取的治疗措施counts <- c(18,17,15,20,10,20,25,13,12) # 生成计数值outcome <- gl(3,1,9) # 生成结局treatment <- gl(3,3) # 生成治疗方案print(d.AD <- data.frame(treatment, outcome, counts))glm.D93 <- glm(counts ~outcome + treatment, family = poisson()) # 指定泊松回归模型summary(glm.D93) # 输出回归结果




这里我们主要看一下相关系数(coefficients),只有outcome2的p值显著并且其效应量值(estimate)是负值,由此可见这三种药之间的效果可能差异不大,并且都能使患者受益。那么只能说这两个新药和现行药的疗效差不多,并不是新药的效果更好。


当然,如果拟合模型的残差比自由度大很多,这个时候最好使用quasipossion()


关于广义线性回归模型的应用就先分享到这里,希望大家持续关注【生信与临床】!