權(quán)重背景：

nhanes的權(quán)重由3部分組成：
基礎(chǔ)權(quán)重basic weight：與常規(guī)抽樣數(shù)據(jù)類似，對應(yīng)樣本被抽中的概率基跑。
無應(yīng)答調(diào)整 non-response adjustments
分層后調(diào)整 post-stratification adjustments

基礎(chǔ)權(quán)重：
普通的問卷調(diào)查利职，為了使抽樣樣本具有一定的人群代表性吉拳，會賦予其一定的權(quán)重，相當(dāng)于樣本中的各個觀測（本研究為每個研究對象）在總體中的重要程度，或樣本中每個觀測所能代表的總體中的個體的數(shù)目紊婉。
無應(yīng)答調(diào)整：
之所以有這個權(quán)重，其原因是nhanes每個周期的調(diào)查中杖挣，participants參與項目的完整度不一致肩榕，且一些特定項目只有部分人員參加。
分層后調(diào)整：
權(quán)重還進行后分層惩妇，以匹配每個抽樣子域的總體控制總數(shù)。這一額外調(diào)整使加權(quán)計數(shù)與美國非機構(gòu)化平民人口的獨立估計相同筐乳。

權(quán)重的計算:

1歌殃，基于合并周期數(shù)計算
2，基于最小子集結(jié)合研究目的選擇合適權(quán)重
本質(zhì)：不管怎樣挑選樣本都能夠代表美國整體人群

基線表

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)


# 正態(tài)性檢驗(樣本超5000)
data <- scale(out.put$LBXBCD)
ks.test(data,"pnorm")  #<0.05氓皱，非正態(tài)分布
data <- scale(out.put$LBXBPB)
ks.test(data,"pnorm")  #<0.05，非正態(tài)分布
data <- scale(out.put$LBXBSE)
ks.test(data,"pnorm")  #<0.05勃刨，非正態(tài)分布
data <- scale(out.put$LBXBMN)
ks.test(data,"pnorm")  #<0.05波材，非正態(tài)分布

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包無法進行加權(quán)限制性立方樣條分析，需要借助rms包身隐。先利用survey包進行加權(quán)l(xiāng)ogistic回歸分析廷区，并提取出標(biāo)準(zhǔn)化權(quán)重然后使用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回歸列線圖(患者風(fēng)險得分)

trl <- sample(nrow(out.put), 0.7*nrow(out.put))
train_data <- out.put[trl,]
test_data <- out.put[-trl,]

# 生成復(fù)雜抽樣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)

# 生成預(yù)測變量，這里生成的是log(p)贾铝，如果要p的話需要轉(zhuǎn)換
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為結(jié)局變量建立一個線性回歸方程
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

# 可信區(qū)間隙轻，se等于sd的一半
se <- 0.026/2
ul <- 0.871 + 1.96*se
ll <- 0.871 - 1.96*se


# 校準(zhǔn)曲線數(shù)據(jù)準(zhǔn)備gg2函數(shù)
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"))


# 繪制校準(zhǔn)曲線
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
# 驗證集校準(zhǔn)曲線
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)

# 生成復(fù)雜抽樣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回歸列線圖(患者風(fēng)險得分)

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,  #由于傳了該文件埠帕，分類變量需要提前調(diào)整為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)