02基于NHANES数据库的加权基线表及加权模型
2023-05-22 本文已影响0人
Jachin111
权重背景：
nhanes的权重由3部分组成：
基础权重basic weight：与常规抽样数据类似，对应样本被抽中的概率。
无应答调整 non-response adjustments
分层后调整 post-stratification adjustments
基础权重：
普通的问卷调查，为了使抽样样本具有一定的人群代表性，会赋予其一定的权重，相当于样本中的各个观测（本研究为每个研究对象）在总体中的重要程度，或样本中每个观测所能代表的总体中的个体的数目。
无应答调整：
之所以有这个权重，其原因是nhanes每个周期的调查中，participants参与项目的完整度不一致，且一些特定项目只有部分人员参加。
分层后调整：
权重还进行后分层，以匹配每个抽样子域的总体控制总数。这一额外调整使加权计数与美国非机构化平民人口的独立估计相同。
权重的计算:
1，基于合并周期数计算
2，基于最小子集结合研究目的选择合适权重
本质：不管怎样挑选样本都能够代表美国整体人群
基线表
library(haven)
library(Hmisc)
library(tableone)
library(survey)
bcSvy2 <- svydesign(ids=~SDMVPSU, #集群，抽样单元PSU
strata=~SDMVSTRA, #分层
weights=~WEIGHT,
nest=TRUE,
survey.lonely.psu="adjust", #抽样单元为1时不报错
data=out.put)
dput(names(out.put))
allVars <- c("RIAGENDR", "RIDAGEYR", "RIDRETH1", "DMDMARTL", "DMDEDUC2", "DMDFMSIZ",
"INDFMIN2","BMXBMI", "SMQ020", "ALQ", "DIQ010", "BPQ020", "PA", "LBXBCD",
"LBXBPB", "LBXBSE", "LBXBMN")
fvars <- c("RIAGENDR", "RIDAGEYR", "RIDRETH1", "DMDMARTL", "DMDEDUC2", "DMDFMSIZ",
"INDFMIN2","BMXBMI", "SMQ020", "ALQ", "DIQ010", "BPQ020", "PA")
Svytab2 <- svyCreateTableOne(vars=allVars, strata="PHQ9", data=bcSvy2, factorVars=fvars)
# 正态性检验(样本超5000)
data <- scale(out.put$LBXBCD)
ks.test(data,"pnorm") #<0.05，非正态分布
data <- scale(out.put$LBXBPB)
ks.test(data,"pnorm") #<0.05，非正态分布
data <- scale(out.put$LBXBSE)
ks.test(data,"pnorm") #<0.05，非正态分布
data <- scale(out.put$LBXBMN)
ks.test(data,"pnorm") #<0.05，非正态分布
bb <- c("LBXBCD", "LBXBPB", "LBXBSE", "LBXBMN")
tab3 <- print(Svytab2, factorVars=fvars, nonnormal=bb, quote=T, showAllLevels=T, smd=T, missing=F)
logistic回归分析
svy.fit <- svyglm(LUXSMED ~ ratio + factor(RIAGENDR) + factor(RIDRETH1) + factor(DMDFMSIZ) + factor(INDFMIN2) + factor(BMXBMI) + RIDAGEYR + factor(Smoke) + factor(ALQ111) + factor(diabetes) + LBXSATSI + LBXSASSI + LBXSGTSI + LBXSTR + LBXSCH, family="quasibinomial", design=bcSvy2)
tbl_regression(svy.fit, exponentiate=TRUE) %>%
add_global_p(anova_fun=gtsummary::tidy_wald_test) %>%
add_glance_table(include=c(nobs, AIC, BIC))
基于logistic回归限制性立方样条
分析逻辑：survey包无法进行加权限制性立方样条分析，需要借助rms包。先利用survey包进行加权logistic回归分析，并提取出标准化权重然后使用rms包进行分析
library(gtsummary)
library(survey)
library(haven)
library(rms)
svy.fit <- svyglm(PHQ9 ~ rcs(LBXBCD,4) + factor(RIAGENDR) + factor(RIDAGEYR) + factor(RIDRETH1) + factor(DMDMARTL) + factor(DMDEDUC2) + factor(DMDFMSIZ) + factor(INDFMIN2) + factor(BMXBMI) + factor(SMQ020) + factor(ALQ) + factor(DIQ010) + factor(BPQ020), family="quasibinomial", design=bcSvy2)
data <- svy.fit$survey.design$variables
ori.weight <- 1/(svy.fit$survey.design$prob)
mean.weight <- mean(ori.weight)
data$weights <- ori.weight/mean.weight
dd <- rms::datadist(data)
options(datadist="dd")
data$PHQ9 <- factor(data$PHQ9)
data$RIAGENDR <- factor(data$RIAGENDR)
data$RIDAGEYR <- factor(data$RIDAGEYR)
data$RIDRETH1 <- factor(data$RIDRETH1)
data$DMDMARTL <- factor(data$DMDMARTL)
data$DMDEDUC2 <- factor(data$DMDEDUC2)
data$DMDFMSIZ <- factor(data$DMDFMSIZ)
data$INDFMIN2 <- factor(data$INDFMIN2)
data$BMXBMI <- factor(data$BMXBMI)
data$SMQ020 <- factor(data$SMQ020)
data$ALQ <- factor(data$ALQ)
data$DIQ010 <- factor(data$DIQ010)
data$BPQ020 <- factor(data$BPQ020)
fit.rcs <- rms::Glm(PHQ9 ~ rcs(LBXBCD,4) + RIAGENDR + RIDAGEYR + RIDRETH1 + DMDMARTL + DMDEDUC2 + DMDFMSIZ + INDFMIN2 + BMXBMI + SMQ020 + ALQ + DIQ010 + BPQ020, data=data, family="quasibinomial", weights=weights, normwt=TRUE, control=list(eps=1e-8))
# weights=weights,normwt=TRUE
# AIC的值
extractAIC(fit.rcs)
# 非线性检验
anova(fit.rcs)
OR <- rms::Predict(fit.rcs, LBXBCD, type="predictions", fun=exp, ref.zero=TRUE)
pdf("Mn_rcs.pdf")
ggplot() +
geom_line(data=OR, aes(x=LBXBCD, y=yhat), linetype="solid", size=1, alpha=0.7, color="red") +
geom_ribbon(data=OR, aes(x=LBXBCD, ymin=lower, ymax=upper), alpha=0.1, fill="blue") +
geom_hline(yintercept=1, linetype=2, color="grey") +
scale_x_continuous("Cd(ug/L)") +
scale_y_continuous("OR (95% CI)") +
geom_text(aes(x=0.5, y=3,
label=paste0("P-overall = 0.0352", "\nP-non-linear = 0.0771")), hjust=0) +
theme_bw() +
theme(axis.line=element_line(),
panel.grid=element_blank(),
panel.border=element_blank())
dev.off()
logistic回归列线图(患者风险得分)
trl <- sample(nrow(out.put), 0.7*nrow(out.put))
train_data <- out.put[trl,]
test_data <- out.put[-trl,]
# 生成复杂抽样NHANES_design
bcSvy2 <- svydesign(data=train_data,
ids=~SDMVPSU,
strata=~SDMVSTRA,
nest=TRUE,
weights=~WTINT2YR,
survey.lonely.psu="adjust")
summary(bcSvy2)
svy.fit <- svyglm(MORTSTAT ~ factor(condition) + factor(RIAGENDR) + factor(DMDEDUC2) + factor(RIDRETH1) + RIDAGEYR + INDFMPIR, family="quasibinomial", design=bcSvy2)
# 生成预测变量，这里生成的是log(p)，如果要p的话需要转换
library(readr)
library(dplyr)
library(plyr)
library(tidyverse)
library(gtsummary)
library(tidyr)
library(survey)
library(haven)
library(rms)
library(SvyNom)
pr1 <- predict(svy.fit)
pr1 <- as.numeric(pr1)
# 以pr1为结局变量建立一个线性回归方程
f <- ols(pr1 ~ condition + RIAGENDR + DMDEDUC2 + RIDRETH1 + RIDAGEYR + INDFMPIR, sigma=1, x=TRUE, y=TRUE, data=train_data)
pr2 <- predict(f)
# 列线图
pdf("norm_logistic.pdf", width=12, height=6)
dd <- datadist(train_data)
options(datadist="dd")
ss3 = c(0.05,0.2,0.4,0.6,0.7,0.8,0.9,0.95,0.99)
nom <- nomogram(f, fun=function(x)1/(1+exp(-x)), funlabel="Risk", fun.at=ss3, lp=TRUE, vnames="labels")
plot(nom)
dev.off()
# C-index
Cindex <- rcorrcens(train_data$MORTSTAT ~ predict(svy.fit))
Cindex
# 可信区间，se等于sd的一半
se <- 0.026/2
ul <- 0.871 + 1.96*se
ll <- 0.871 - 1.96*se
# 校准曲线数据准备gg2函数
gg2 <- function(data, p, y, group=1, g, leb){
if (missing(g)){g=10}
matres <- function(data, p, y){
data <- data
y=y
sor <- order(p)
p <- p[sor]
y <- y[sor]
groep <- cut2(p, g=g)
meanpred <- round(tapply(p, groep, mean), 3)
meanobs <- round(tapply(y, groep, mean), 3)
msd <- round(tapply(y, groep, sd), 3)
se <- msd/sqrt(dim(data)[1])
matres <- as.data.frame(cbind(meanpred, meanobs, se))
matres
}
dat <- matres(data, p, y)
if (group==1) dat
if (group!=1){
dat <- cbind(dat, gro=leb)
}
dat
}
p <- exp(pr1)/(1+exp(pr1)) #或p <- predict(svy.fit, type=c("response"))
# 绘制校准曲线
pdf("calibrate_logistic.pdf", width=12, height=6)
plot1 <- gg2(train_data, p, train_data$MORTSTAT)
ggplot(plot1, aes(x=meanpred, y=meanobs)) +
geom_errorbar(aes(ymin=meanobs - 1.96*se, ymax=meanobs + 1.96*se), width=.02) +
annotate(geom="segment", x=0, y=0, xend=1, yend=1) +
expand_limits(x=0, y=0) +
scale_x_continuous(expand=c(0,0)) +
scale_y_continuous(expand=c(0,0)) +
geom_point(size=3, shape=21, fill="white") +
xlab("Prediction probability") +
ylab("Actual probability")
dev.off()
# ROC曲线
library(ResourceSelection)
library(pROC)
pr.roc <- predict(svy.fit, type=c("response"))
h2 <- hoslem.test(train_data$MORTSTAT, pr.roc, g=10)
h2 #p>0.05较好
cbind(h2$observed, h2$expected)
roccurvel <- roc(train_data$MORTSTAT ~ pr.roc)
auc(roccurvel)
pdf("train_ROC.pdf")
plot.roc(roccurvel, xlim=c(1,0), ylim=c(0,1), print.auc=TRUE, print.auc.x=0.4, print.auc.y=0.2, print.thres="best")
dev.off()
# 在验证集上生成p和log(p)
bc_best.pr1 <- predict(svy.fit, newdata=test_data)
bc_best.pr1 <- as.numeric(bc_best.pr1)
bc_best.p <- predict(svy.fit, newdata=test_data, type=c("response"))
# 验证集上求C-index
Cindex <- rcorrcens(test_data$MORTSTAT ~ predict(svy.fit, newdata=test_data))
Cindex
# 验证集校准曲线
pdf("calibrate_logistic_test.pdf", width=12, height=6)
plot2 <- gg2(test_data, bc_best.p, test_data$MORTSTAT)
ggplot(plot2, aes(x=meanpred, y=meanobs)) +
geom_errorbar(aes(ymin=meanobs - 1.96*se, ymax=meanobs + 1.96*se), width=.02) +
annotate(geom="segment", x=0, y=0, xend=1, yend=1) +
expand_limits(x=0, y=0) +
scale_x_continuous(expand=c(0,0)) +
scale_y_continuous(expand=c(0,0)) +
geom_point(size=3, shape=21, fill="white") +
xlab("Prediction probability") +
ylab("Actual probability")
dev.off()
# 验证集ROC
roccurvel2 <- roc(test_data$MORTSTAT ~ bc_best.p)
pdf("test_ROC.pdf")
auc(roccurvel2)
plot(roccurvel2, xlim=c(1,0), ylim=c(0,1), print.auc=TRUE, print.auc.x=0.4, print.auc.y=0.2, print.thres="best")
dev.off()
# 验证集列线图绘制
f1 <- ols(bc_best.pr1 ~ condition + RIAGENDR + DMDEDUC2 + RIDRETH1 + RIDAGEYR + INDFMPIR, sigma=1, x=TRUE, y=TRUE, data=test_data)
pr3 <- predict(f1)
pdf("norm_logistic_test.pdf", width=12, height=6)
dd <- datadist(test_data)
options(datadist="dd")
ss3 <- c(0.05,0.2,0.4,0.6,0.7,0.8,0.9,0.95,0.99)
nom <- nomogram(f, fun=function(x)1/(1+exp(-x)), funlabel="Risk", fun.at=ss3, lp=TRUE, vnames="labels")
plot(nom)
dev.off()
KM分析
library(readr)
library(dplyr)
library(plyr)
library(tidyverse)
library(gtsummary)
library(tidyr)
library(survey)
library(haven)
library(survival)
library(survminer)
library(jskm)
library(rms)
library(SvyNom)
library(timeROC)
# 生成复杂抽样NHANES_design
bcSvy2 <- svydesign(data=out.put,
ids=~SDMVPSU,
strata=~SDMVSTRA,
nest=TRUE,
weights=~WTINT2YR,
survey.lonely.psu="adjust")
summary(bcSvy2)
sl <- svykm(Surv(PERMTH_INT, MORTSTAT) ~ factor(condition), design=bcSvy2, se=TRUE)
pdf("survival.pdf")
svyjskm(sl, pval=TRUE, table=TRUE, showpercent=TRUE)
dev.off()
cox回归分析
ml1 <- svycoxph(Surv(PERMTH_INT, MORTSTAT) ~ factor(condition) + factor(RIAGENDR) + factor(DMDEDUC2) + factor(RIDRETH1) + RIDAGEYR + INDFMPIR, x=TRUE, design=bcSvy2)
summary(ml1)
tbl_regression(ml1, exponentiate=TRUE)
基于cox回归限制性立方样条
bcSvy2 <- svydesign(data=out.put,
ids=~SDMVPSU,
strata=~SDMVSTRA,
nest=TRUE,
weights=~WTINT2YR,
survey.lonely.psu="adjust")
summary(bcSvy2)
svy.fit <- svycoxph(Surv(PERMTH_INT, MORTSTAT) ~ rcs(RIDAGEYR,4) + factor(condition) + factor(RIAGENDR) + factor(DMDEDUC2) + factor(RIDRETH1) + INDFMPIR, x=TRUE, design=bcSvy2)
data <- svy.fit$survey.design$variables
ori.weight <- 1/(svy.fit$survey.design$prob)
mean.weight <- mean(ori.weight)
data$weights <- ori.weight/mean.weight
dd <- rms::datadist(data)
options(datadist="dd")
data$condition <- as.factor(data$condition)
data$RIAGENDR <- as.factor(data$RIAGENDR)
data$DMDEDUC2 <- as.factor(data$DMDEDUC2)
data$RIDRETH1 <- as.factor(data$RIDRETH1)
fit.rcs <- rms::cph(Surv(PERMTH_INT, MORTSTAT) ~ rcs(RIDAGEYR,4) + condition + RIAGENDR + DMDEDUC2 + RIDRETH1 + INDFMPIR, data=data, weights=weights, normwt=TRUE, control=list(eps=1e-8))
# AIC的值
extractAIC(fit.rcs)
# 非线性检验
anova(fit.rcs)
HR <- rms::Predict(fit.rcs, RIDAGEYR, type="predictions", fun=exp, ref.zero=TRUE)
ggplot() +
geom_line(data=HR, aes(x=RIDAGEYR, y=yhat), linetype="solid", size=1, alpha=0.7, color="red") +
geom_ribbon(data=HR, aes(x=RIDAGEYR, ymin=lower, ymax=upper), alpha=0.1, fill="blue") +
geom_hline(yintercept=1, linetype=2, color="grey") +
scale_x_continuous("Age") +
scale_y_continuous("HR (95% CI)") +
geom_text(aes(x=0.5, y=5,
label=paste0("P-overall < 0.0001", "\nP-non-linear = 0.7178")), hjust=0) +
theme_bw() +
theme(axis.line=element_line(),
panel.grid=element_blank(),
panel.border=element_blank())
dev.off()
HR1 <- Predict(fit.rcs, RIDAGEYR, RIAGENDR=c(1,2),
fun=exp, type="predictions",
ref.zero=TRUE, conf.int = 0.95, digits=2)
HR1$RIAGENDR <- as.factor(HR1$RIAGENDR)
ggplot() +
geom_line(data=HR1, aes(x=RIDAGEYR, y=yhat, group=RIAGENDR, col=RIAGENDR)) +
geom_ribbon(data=HR1, aes(x=RIDAGEYR, ymin=lower, ymax=upper, fill=RIAGENDR), alpha=0.3) +
geom_hline(yintercept=1, linetype=2, size=1) +
theme_classic() +
labs(title="RCS", x="age", y="HR (95%CI)")
ggplot()+
geom_line(data=HR1, aes(RIDAGEYR, yhat, color=RIAGENDR),
linetype="solid", size=1, alpha=0.7) +
geom_ribbon(data=HR1,
aes(RIDAGEYR, ymin=lower, ymax=upper, fill=RIAGENDR),
alpha=0.1)+
scale_color_manual(values = c('#0070b9','#d40e8c')) +
scale_fill_manual(values = c('#0070b9','#d40e8c')) +
theme_classic()+
geom_hline(yintercept=1, linetype=2,size=1)+
labs(title = "RCS", x="Age", y="HR (95%CI)")
cox回归列线图(患者风险得分)
out.put$condition <- as.factor(out.put$condition)
out.put$RIAGENDR <- as.factor(out.put$RIAGENDR)
out.put$DMDEDUC2 <- as.factor(out.put$DMDEDUC2)
out.put$RIDRETH1 <- as.factor(out.put$RIDRETH1)
dd <- datadist(out.put)
options(datadist="dd")
ss <- c(0.05,0.2,0.4,0.6,0.7,0.8,0.9,0.95,0.99)
mynom <- svycox.nomogram(.design=bcSvy2,
.model=Surv(PERMTH_INT, MORTSTAT) ~ factor(condition) + factor(RIAGENDR) + factor(DMDEDUC2) + factor(RIDRETH1) + RIDAGEYR + INDFMPIR,
.data=out.put, #由于传了该文件，分类变量需要提前调整为factor
pred.at=12*2,
fun.lab="Prob of 3 Yr OS")
pdf("norm.pdf", width=12, height=6)
plot(mynom$nomog)
dev.off()
# 校正曲线
pdf("calibrate.pdf", width=12, height=6)
svycox.calibrate(mynom)
dev.off()
bootit=200
library(survey)
library(rms)
data(noNA)
dd=datadist(noNA)
options(datadist="dd")
dstr2=svydesign(id=~1, strata=~group, prob=~inv_weight, fpc=~ssize, data=noNA)
mynom=svycox.nomogram(.design=dstr2, .model=Surv(survival,surv_cens)~ECOG+liver_only+Alb+Hb+Age+
Differentiation+Gt_1_m1site+lymph_only, .data=noNA, pred.at=24, fun.lab="Prob of 2 Yr OS")
cases=which(noNA$group=="long")
controls=which(noNA$group=="<24")
boot.index=matrix(NA,nrow(noNA),bootit)
for(i in 1:bootit){
boot.index[,i]=c(sample(cases,replace=TRUE),sample(controls,replace=TRUE))
}
myval=svycox.validate(boot.index,mynom,noNA)
