R语言（二）——简单线性模型中的指数变换

目录

一、数据

1.数据信息

2.数据处理

二、简单线性回归

三、指数变换

四、生存分析数据的Cox回归模型

---

一、数据

1.数据信息

口咽癌数据（pharynx.csv）是针对口咽若干位置癌细胞的临床实验。分成两组，一组仅使用放疗（TX=1），另一组使用放疗和化疗（TX=2）。来自Kalbfleisch and Prentice(1980)，原本可从http://www.umass.edu/statdata/statdata/stat-nonlin.html下载，但现在似乎已经不可以了，我把资源上传了，在评论区大家可以自取。【评论区有小伙伴提供了github的下载地址，不知道和我的数据有没有区别，不过那边下载不用积分，大家可以尝试一下】

2.数据处理

 数据中的CASE和ENTRY.DT属于无用数据，进行删除。

部分观测值有缺失，进行删除（数据中用9表示缺失）

COND中部分数据过少，进行合并

w=read.csv("pharynx.csv")
w=w[,-c(1,11)]  #去掉CASE和日期
w=w[w$COND!=9&w$GRADE!=9,] #去掉两列NA
w$COND[w$COND==3|w$COND==4]=2
w$COND[w$COND==0]=1
#存储
write.csv(w,"pharynx1.csv",quote=F,row.names=F)

二、简单线性回归

除TIME和AGE外，其他为定性（分类）变量，需要把它们将数值型转换为因子型。

TIME为我们的因变量，进行简单线性回归

u=read.csv("pharynx1.csv")
x=1:11;(x=x[-c(5,11)])  #定性变量的下标
for(i in x)u[,i]=factor(u[,i]) #把定性变量从数值型转换成因子型
a=lm(TIME~.,data=u);summary(a)
shapiro.test(a$res)

结果显示R²只有0.577，正态性检验值1.569e-06。不满足正态性

三、指数变换

尝试对TIME进行指数变换，利用MASS（）包的boxcox（）函数可以找到回归中变换的参数

library(MASS)
b=boxcox(TIME~.,data=u)
I=which(b$y==max(b$y))
b$x[I]   

得到不同值的对数似然曲线，其中最大的为0.38，我们可以取0.4。对变换后的因变量做逐步回归

a=lm(TIME^0.4~INST+SEX+TX+AGE+COND+T.STAGE+N.STAGE+STATUS,data=u)
b=step(a)
summary(b) #anova(b)
shapiro.test(b$res)

此时R²为0.5618，调整的R²为0.5377。残差的正态性检验P值0.1473，虽不能否认，但依然很难有效地说明。

四、生存分析数据的Cox回归模型

生存分析中，研究的主要对象是寿命超过某一时间的概率，也可以描述其他一些事件发生的概率，如产品的失效、出狱犯人第一次犯罪、失业人员第一次找到工作、青少年第一次吸毒等。

利用本数据，如果仅考虑变量TX而不考虑其他各因素的影响，可以画出根据数据得到的综合生存函数图。画图使用了survival包。

library(survival)
fit=survfit(Surv(TIME,as.numeric(STATUS))~TX,data=u)
#口咽癌数据的综合生存函数
plot(fit,lty=1:2,ylab="S(t)",xlab="t",main='Survival Functions')
legend("top",c("TX=1","TX=2"),lty=1:2)

由于生存函数不是一成不变的，可能是一些自变量的函数。

Cox比例危险回归模型把生存函数的某种变换形式描述成自变量的线性函数。

如果用X代表自变量向量，而用危险函数或其对数作为响应变量，则Cox比例危险回归模型可写成

或

library(survival)
fit=coxph(Surv(TIME,as.numeric(STATUS))~.,data=u)
#口咽癌数据的综合生存函数
plot(survfit(fit))
summary(fit)

从输出来看，许多变量并不显著，可以减少一些变量。

进一步用程序包mfp来尝试关于Cox比例危险模型的多重分数多项式模型。

其类似于多项式回归，但更加灵活，且自动进行逐步回归及变量选择

#生存函数变化的的Cox模型
library(survival)
fit=coxph(Surv(TIME,as.numeric(STATUS))~.,data=u)
#口咽癌数据的综合生存函数
plot(survfit(fit))
summary(fit)

此时似然比检验的p值很小，R²为0.999，拟合效果比之前好