寫在前面
之前做群體進化琐凭，經(jīng)常會遇到需要考慮群體結(jié)構(gòu)對結(jié)果的影響设江。比如病例組和對照組之間在某個基因上存在差異，病例組在這個基因上突變了炫隶，對照組沒有，剛開始以為是導致該病的原因阎曹。后來發(fā)現(xiàn)病例組是靠近北極的人伪阶，對照組是靠近赤道的人，本來兩組人群在遺傳上的差異就賊拉噠处嫌，這個突變可能是本來就存在的差異栅贴，不一定就是和疾病相關。這時候就需要將地域信息（人種信息）加入熏迹，作為協(xié)變量對結(jié)果進行校正檐薯。（好像扯遠了...）
一般臨床上要找指示疾病的標志物，比如某個基因的表達水平高低是否與疾病相關注暗，也需要對常見的相關因素進行校正坛缕，如年齡和性別墓猎。
是否患病是0和1的游戲，所以這其實是一個binomial logistic regression的問題祷膳。除了擬合度的高低（R2,p-value）陶衅，更應該關注比值比（odd ratio)的大小。OR>1直晨，說明該因素能提高患病風險搀军；OR<1，則認為是保護因素勇皇。
實戰(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) #查看相應參數(shù)
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
實際數(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) #計算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在病例和對照組之間的比值比是5.04门烂，說明T6是風險因素，可以提高5倍的患病率（表述的好像不太準確兄淫，大概是這個意思）屯远。
SPSS
SPSS可以很方便地以表格形式輸出結(jié)果，唯一的弊端可能是表格太多捕虽，有時候不清楚看哪個慨丐。這個鏈接內(nèi)容寫得很詳細，供參考 https://www.zhihu.com/question/34502688
寫在最后
發(fā)現(xiàn)自己的統(tǒng)計泄私，真的是渣渣級別啊房揭。許多底層的基礎邏輯和思維嚴重欠缺，不說了晌端，看書去吧捅暴。
