写在前面
之前做群体进化，经常会遇到需要考虑群体结构对结果的影响。比如病例组和对照组之间在某个基因上存在差异，病例组在这个基因上突变了，对照组没有，刚开始以为是导致该病的原因。后来发现病例组是靠近北极的人，对照组是靠近赤道的人，本来两组人群在遗传上的差异就贼拉哒，这个突变可能是本来就存在的差异，不一定就是和疾病相关。这时候就需要将地域信息（人种信息）加入，作为协变量对结果进行校正。（好像扯远了...）
一般临床上要找指示疾病的标志物，比如某个基因的表达水平高低是否与疾病相关，也需要对常见的相关因素进行校正，如年龄和性别。
是否患病是0和1的游戏，所以这其实是一个binomial logistic regression的问题。除了拟合度的高低（R2,p-value），更应该关注比值比（odd ratio)的大小。OR>1，说明该因素能提高患病风险；OR<1，则认为是保护因素。
实战呗
例子文件
###
library(epiDisplay)
library(carData)
data(Wells, package="carData")
head(Wells)
# switch arsenic distance education association
#1 yes 2.36 16.826 0 no
#2 yes 0.71 47.322 0 no
#3 no 2.07 20.967 10 no
glm1 <- glm(switch~arsenic+distance+education+association,
family=binomial, data=Wells)
summary(glm1) #查看相应参数
logistic.display(glm1) #计算OR值
#Logistic regression predicting switch : yes vs no
# crude OR(95%CI) adj. OR(95%CI) P(Wald's test) P(LR-test)
#arsenic (cont. var.) 1.46 (1.35,1.58) 1.6 (1.47,1.73) < 0.001 < 0.001
#distance (cont. var.) 0.9938 (0.9919,0.9957) 0.9911 (0.989,0.9931) < 0.001 < 0.001
#education (cont. var.) 1.04 (1.02,1.06) 1.04 (1.02,1.06) < 0.001 < 0.001
#association: yes vs no 0.86 (0.75,1) 0.88 (0.76,1.03) 0.106 0.106
#Log-likelihood = -1953.913
#No. of observations = 3020
#AIC value = 3917.826
pseudo_r2=1-(glm$deviance/glm$null.deviance) #McFadden's pseudo-R2;不同模型用不同的pseudo-R2
关于pseudo r2，可参考以下描述：
The denominator of the ratio can be thought of as the sum of squared errors from the null model--a model predicting the dependent variable without any independent variables.
了解不同的pseudo_r2 https://web.archive.org/web/20130701052120/http://www.ats.ucla.edu:80/stat/mult_pkg/faq/general/Psuedo_RSquareds.htm
实际数据
data<-read.table("adjust_gender_age.inp",header=T)
head(data)
#GROUP AGE SEX T6 T7 T8 T9
#1 0 48 1 0.06 0.56 0.38 0.74
#2 0 45 1 0.07 0.95 0.34 0.77
#3 0 24 1 0.00 0.65 0.40 0.83
#(T6)binomial logistic regression:要求y必须用0和1来分类 (0<=y<=1)
glm<- glm(GROUP~AGE+SEX+T6,
family=binomial, data=data)
pseudo_r2=1-(glm$deviance/glm$null.deviance)
name=cbind(T6,pseudo_r2)
write.table(name,file="glm.out",quote=FALSE,sep="\t",append=TRUE,col.names=FALSE,row.names=FALSE)
or=logistic.display(glm) #计算OR
write.table(or$table,file="glm.out",quote=FALSE,sep="\t",append=TRUE,col.names=FALSE)
输出结果如下：
#T6 0.0128045811182196
#AGE (cont. var.) 0.9995 (0.994,1.0051) 0.9996 (0.9941,1.0053) 0.899 0.899
#SEX (cont. var.) 0.6 (0.52,0.71) 0.61 (0.52,0.72) < 0.001 < 0.001
#T6 (cont. var.) 6.62 (1.55,28.29) 5.04 (1.17,21.6) 0.03 0.025
从结果可以看到，在校正AGE和SEX后，T6在病例和对照组之间的比值比是5.04，说明T6是风险因素，可以提高5倍的患病率（表述的好像不太准确，大概是这个意思）。
SPSS
SPSS可以很方便地以表格形式输出结果，唯一的弊端可能是表格太多，有时候不清楚看哪个。这个链接内容写得很详细，供参考 https://www.zhihu.com/question/34502688
写在最后
发现自己的统计，真的是渣渣级别啊。许多底层的基础逻辑和思维严重欠缺，不说了，看书去吧。
网友评论