寫在前面
之前做群體進(jìn)化，經(jīng)常會(huì)遇到需要考慮群體結(jié)構(gòu)對(duì)結(jié)果的影響。比如病例組和對(duì)照組之間在某個(gè)基因上存在差異，病例組在這個(gè)基因上突變了，對(duì)照組沒有，剛開始以為是導(dǎo)致該病的原因。后來發(fā)現(xiàn)病例組是靠近北極的人，對(duì)照組是靠近赤道的人，本來兩組人群在遺傳上的差異就賊拉噠，這個(gè)突變可能是本來就存在的差異，不一定就是和疾病相關(guān)。這時(shí)候就需要將地域信息（人種信息）加入，作為協(xié)變量對(duì)結(jié)果進(jìn)行校正。（好像扯遠(yuǎn)了...）
一般臨床上要找指示疾病的標(biāo)志物，比如某個(gè)基因的表達(dá)水平高低是否與疾病相關(guān)，也需要對(duì)常見的相關(guān)因素進(jìn)行校正，如年齡和性別。
是否患病是0和1的游戲，所以這其實(shí)是一個(gè)binomial logistic regression的問題。除了擬合度的高低（R2,p-value），更應(yīng)該關(guān)注比值比（odd ratio)的大小。OR>1，說明該因素能提高患病風(fēng)險(xiǎn)；OR<1，則認(rèn)為是保護(hù)因素。
實(shí)戰(zhàn)唄
例子文件
###
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) #查看相應(yīng)參數(shù)
logistic.display(glm1) #計(jì)算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
關(guān)于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
實(shí)際數(shù)據(jù)
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) #計(jì)算OR
write.table(or$table,file="glm.out",quote=FALSE,sep="\t",append=TRUE,col.names=FALSE)
輸出結(jié)果如下：
#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
從結(jié)果可以看到，在校正AGE和SEX后，T6在病例和對(duì)照組之間的比值比是5.04，說明T6是風(fēng)險(xiǎn)因素，可以提高5倍的患病率（表述的好像不太準(zhǔn)確，大概是這個(gè)意思）。
SPSS
SPSS可以很方便地以表格形式輸出結(jié)果，唯一的弊端可能是表格太多，有時(shí)候不清楚看哪個(gè)。這個(gè)鏈接內(nèi)容寫得很詳細(xì)，供參考 https://www.zhihu.com/question/34502688
寫在最后
發(fā)現(xiàn)自己的統(tǒng)計(jì)，真的是渣渣級(jí)別啊。許多底層的基礎(chǔ)邏輯和思維嚴(yán)重欠缺，不說了，看書去吧。
