R语言学习DAY04:回归分析

Python014

R语言学习DAY04:回归分析,第1张

R本身是一门统计语言,主要用于统计分析,前面的语法部分算是基础,接下来开始进入统计模型应用。首先从最常用的回归分析说起。

有关线性回归分析模型的基本假定需要注意:1)关于随机干扰项的高斯-马尔科夫定理;2)关于自变量的:不存在共线性;3)关于模型的:模型设定正确。

用 glm 函数建立广义线性模型,用参数 family 指定分布类型,logistic模型指定为binomial

用 predict 函数进行预测, predict(model, data, type = 'response'

此外,还可以用 mlogit 包中的 mlogit 函数做多分类变量logistic回归, rms 包中的 lrm 函数做顺序变量logistic回归, glmnet 包中的 glmnet 函数做基于正则化的logistic回归

在医学研究中,我们经常构建回归模型来分析自变量和因变量之间的关系。事实上,大多数的回归模型有一个重要的假设就是自变量和因变量呈线性关联,这个条件实际很难满足。常见的解决方法是将连续变量分类,但类别数目和节点位置的选择往往带有主观性,并且分类往往会损失信息。因此,一个更好的解决方法是拟合自变量与因变量之间的非线性关系,限制性立方(Restricted cubic spline,RCS)就是分析非线性关系的最常见的方法之一。

近年来在Lancet、BMJ等杂志经常见到利用限制性立方样条来拟合非线性关系。

什么是立方样条?

回归样条(regression spline)本质上是一个分段多项式, 但它一般要求每个分段点上连续并且二阶可导,这样可以保证曲线的平滑性。而限制性立方样条是在回归样条的基础上附加要求:样条函数在自变量数据范围两端的两个区间内为线性函数。

在利用限制性立方样条绘制曲线关系时,通常需要设置样条函数节点的个数(k)和位置(ti)。绝大多数情况下, 节点的位置对限制性立方样条的拟合影响不大, 而节点的个数则决定曲线的形状, 或者说平滑程度。当节点的个数为2时, 得到的拟合曲线就是一条直线,大多数研究者推荐的节点为3-5个。

在《Regression Modeling Strategies》这本书中,Harrell建议节点数为4时,模型的拟合较好,同时可以兼顾曲线的平滑程度和避免过拟合造成的精度降低。而当样本量较大时,例如因变量为未删失的连续变量并且大于100时,5个节点是更好的选择。小样本(如n<30)可以选择3个节点。以下是Harrell推荐的节点数和相应的节点位置,大家可以参考。

案例说明(模拟数据)

目前SAS、STATA、R等软件都可以进行限制性立方样条分析。基于画图的方便,我们以R语言为例进行说明。首先参照rms包,生成一个模拟数据集,包括性别(sex),年龄(age)以及生存时间(time)和结局变量(death)。

若想分析年龄和生存率之间关系,传统的方法可以在Cox回归中将年龄作为连续变量处理,也可以对年龄进行分组,这样的做法都无法更直观的呈现年龄与死亡风险之间的关联。以下我们用限制性立方样条来分析年龄与死亡风险之间的关系:

可以看到age整体是有意义的(包括线性或者非线性关联),然后看P-Nonlinear =0.0168<0.05,这里我们可以说年龄与死亡风险之间存在非线性关联。

如果自变量与关注的结局变量存在非线性关系,如何在文章中对结果更详细的描述呢,建议大家可以参照上文中提到的Lancet的文章。

rm(list=ls())##清除之前数据

library(rms)  ###加载rms包#

library(foreign)

library(survival)

library(rmda)

##install.packages("DynNom")

library(DynNom)

setwd("C://Users//zhanglingyu//Desktop//测试")  #设置工作目录

rt<- read.table('测试.txt',header=T,sep="\t")

head(rt)  #查看前5行的数据

str(rt)

##接着对所有变量数据进行打包

ddist <- datadist(rt)  #将数据打包

options(datadist='ddist')

model1<- lrm(rt$fustat ~ Age+Gender+Tau+EMPG+Glu+HCY+VB12+ALB, data =  rt)

summary(model1)

par(mgp=c(1.6,0.6,0),mar=c(2,2,2,2))  ##设置画布

nomogram <- nomogram(model1,fun=function(x)1/(1+exp(-x)), ##逻辑回归计算公式

                    fun.at = c(0.001,0.01,0.05,seq(0.1,0.9,by=0.1),0.95,0.99,0.999),#风险轴刻度

                    funlabel = "Risk of Alzheimer Disease", #风险轴便签

                    lp=T,  ##是否显示系数轴

                    conf.int = F, ##每个得分的置信度区间,用横线表示,横线越长置信度越

                    abbrev = F#是否用简称代表因子变量

)

plot(nomogram,

    #1.变量与图形的占比

    xfrac=.35,

    #2.变量字体加粗

    cex.var=1,

    #3.数轴:字体的大小

    cex.axis=0.8,

    #4.数轴:刻度的长度

    tcl=-0.5,

    #5.数轴:文字与刻度的距离

    lmgp=0.3,

    #6.数轴:刻度下的文字,1=连续显示,2=隔一个显示一个

    label.every=1,

    #7.1个页面有几个数轴(这个可以压缩行间距)

    naxes=13,

    #8.垂直线的颜色.

    col.grid=gray(c(0.8, 0.95)),

    #9.线性预测轴名字

    lplabel="Linear Predictorlp",

    #10变量分数名字

    points.label='Points',

    #11总分名字

    total.points.label='Total Points',

    force.label=F#没啥用TRUE强制标记的每个刻度线都绘制标签,我也没研究明白

)

#运行以展示Nomogram每个变量的分数

model1

##生成改良诺模图

plot(nomogram,col.grid = c("Tomato2","DodgerBlue"))