image.png

#分类变量（非正态分布）：广义线性模型
#相关性不能告知因果关系，也不能预测曲线关系
#相关系数r越小，线性关系越弱
library(rcompanion)
library(rms)
library(ResourceSelection)
library(VGAM)
library(magrittr)

Framingham <- read.csv(file = "R语言入门必学 资料/2.20.Framingham.csv")
any(is.na(Framingham)) #anyNA(Framingham) 判断是否有缺失
Framingham <- na.exclude(Framingham) #剔除缺失值

Framingham$sex <- factor(Framingham$sex)
attach(Framingham)

1 经典Logistic回归(二项分类)

1.1 建模

mod1 <- glm(chdfate~sex, family = binomial(link = "logit"))
#单因素 
#family指定概率分布的类型 binomial:二项分布; logit: logistic回归
#得到AIC值

glm(chdfate~sex,family = gaussian())
#线性正态分布，等价于lm()

summary(mod1)
#得到P 

mod2 <- glm(chdfate~.,data = Framingham,family = binomial())
summary(mod2)

mod3 <- glm(chdfate~.,data = Framingham[,-3],family=binomial()) #去掉第3列的因素

anova(mod2,mod3,test = "Chisq") 
#返回模型的方差分析表
#比较两个模型是否显著差异

1.2 逐步法筛选自变量

  mod.none <- glm(chdfate~1,family = binomial)
  mod.full <- glm(chdfate~.,data = Framingham,family = binomial())
  step(mod.none,scope = formula(mod.full),direction = "both",trace = F) %>% summary()

1.3 模型提供的信息

OR

##OR：比值比，疾病与暴露之间关联强度指标
exp(coef(mod3))
#coef()提取模型系数,coefficients（r2?）
#exp():计算e的幂
nagelkerke(mod3)
#计算r2（拟合优度）
#最大似然检验的p值（Likelihood.ratio.test）：p<0.05：与空模型相比，有显著性差异
hoslem.test(chdfate,fitted(mod3))
#计算P值（观察值/因变量，fitted预测值：各个样本为阳性结果的概率？）
#P>0.05，观察值与预测值没有差异，模型预测得好

1.4 作图

列线图

ddist <- datadist(Framingham)
options(datadist="ddist")
#把数据的基本统计特征放在拟合模型中
mod4 <- lrm(chdfate ~sex + sbp+scl+age+bmi,data = Framingham)
#nomogram输入模式
#除此之外，仍用glm(),lrm()结果难读懂
mod4$coefficients
nom <- nomogram(mod4,fun=plogis, funlabel = "CHD Fate", lp=F)
plot(nom)

矫正曲线

mod4 <- lrm(chdfate ~sex + sbp+scl+age+bmi,data = Framingham,x=T,y=T)
#x,y:保存自变量矩阵和因变量  mod4$x,mod4$y

cal <- rms::calibrate(mod4, method="boot",B=1000, bw =T,rule="p",sls = 0.05)
#rule:p值作为判断准则；aic:aic作为准则
plot(cal)
#apparent 自己做的模型，实线是bootstrap预测的模型

2 多分类logistic回归

##例子：好 中 差
Framingham$scl.f <- cut(scl,quantile(scl),labels = c("A","B","C","D"),ordered_result = T, include.lowest = T)
#ordered_result:有序变量
#include.lowest:闭区间
attach(Framingham)

mod.p <- vglm(scl.f~sex+sbp+age+bmi,family = cumulative(reverse=T,parallel = T))
#reverse=F：计算样本落在小于等于某一层的概率
#各自变量的斜率不变，parallel=T

summary(mod.p)

2.2 平行性检验

mod.np <- vglm(scl.f~sex+sbp+age+bmi,family = cumulative(reverse=T,parallel = F))
lrtest(mod.p,mod.np)
#p>0.05，表示符合平行假设

mod.n <- vglm(scl.f~1,family = cumulative(reverse=T,parallel = F))
lrtest(mod.p,mod.n)
#检测跟空模型是否有差异

#若p<0.05, ....