机器算法验证 - 国家一级流行率的荟萃分析 - 吾爱随笔录

我正在对流行数据进行荟萃分析。目的是获得国家一级的流行率估计值。主要问题是疾病与年龄高度相关，纳入研究的样本年龄高度异质。大多数研究只有中位年龄，所以我不能使用类似 SMR 的技巧。我想我可以使用元回归来解决这个问题，包括将年龄作为固定效应并引入研究层面和国家层面的随机效应。

这个想法（我从Fowkes 等人那里得到）是使用这个模型对 15 到 60 岁的每个 5 岁年龄组（使用该组的中位年龄）进行特定国家/地区的流行率预测，并应用这些预测与所选国家中每个群体的实际人口规模，以获得总感染人口并据此计算 15-60 岁人口的年龄调整患病率。

我尝试了几种方法来使用带有包的 Rmeta和mgcv. 我得到了一些令人满意的结果，但我对我的结果没有那么自信，希望得到一些反馈。

首先是一些模拟数据，然后是我的不同方法的描述：

data<-data.frame(id_study=c("UK1","UK2","UK3","FRA1","FRA2","BEL1","GER1","GER2","GER3"),
                 country=c("UK","UK","UK","FRANCE","FRANCE","BELGIUM","GERMANY","GERMANY","GERMANY"),
                 n_events=c(91,49,18,10,50,6,9,10,22),
                 n_total=c(3041,580,252,480,887,256,400,206,300),
                 study_median_age=c(25,50,58,30,42,26,27,28,36))

标准随机效应荟萃分析包meta。

我过去常常metaprop()在不考虑年龄的情况下对每个国家的流行率进行初步估计，并获得权重。正如预期的那样，异质性非常高，所以我使用了随机效应模型的权重。

 meta <- metaprop(event=n_events,n=n_total,byvar=country,sm="PLOGIT",method.tau="REML",data=data)
 summary(meta)
 data$weight<-meta$w.random

我使用 meta 在不考虑年龄的情况下对流行率进行了初步估计，并获得了权重。正如预期的那样，异质性非常高，所以我使用了随机效应模型的权重。

广义的加法模型，包括年龄和包mgcv。

使用gam()BIC 和 GCV 数（此处未显示）选择模型参数（k 和 sp）。

 model <- gam( cbind(n_events,n_total-n_events) ~ s(study_median_age,bs="cr",k=4,sp=2) + s(country,bs="re"), weights=weight, data=data, family="binomial"(link=logit), method="REML")
 plot(model,pages=1,residuals=T, all.terms=T, shade=T)

如前所述，每个年龄组的预测都是从该模型中获得的。CI 直接使用获得predict.gam()，即使用参数的贝叶斯后验协方差矩阵。例如考虑英国：

 newdat<-data.frame(country="UK",study_median_age=seq(17,57,5))
 link<-predict(model,newdat,type="link",se.fit=T)$fit
 linkse<-predict(model,newdat,type="link",se.fit=T)$se
 newdat$prev<-model$family$linkinv(link)
 newdat$CIinf<-model$family$linkinv(link-1.96*linkse)
 newdat$CIsup<-model$family$linkinv(link+1.96*linkse)
 plot(newdat$prev~newdat$study_median_age, type="l",ylim=c(0,.12))
 lines(newdat$CIinf~newdat$study_median_age, lty=2)
 lines(newdat$CIsup~newdat$study_median_age, lty=2)

结果令人满意，表明患病率随着年龄的增长而增加，具有一致的置信区间。我使用国家人口结构获得了该国的总患病率（未显示，我希望它足够清楚）。

但是，我认为我需要包括研究级别的随机效应，因为存在很高的异质性（即使我没有在元回归后计算异质性）。

使用 package引入研究级随机效应gamm4。

由于mgcv模型无法处理那么多随机效应参数，我不得不切换到gamm4.

 model2 <- gamm4(cbind(n_events,n_total-n_events) ~ s(study_median_age,bs="cr",k=4) + s(country,bs="re"), random=~(1|id_study), data=data, weights=weight, family="binomial"(link=logit))
 plot(model2$gam,pages=1,residuals=T, all.terms=T, shade=T)

 link<-predict(model2$gam,newdat,type="link",se.fit=T)$fit
 linkse<-predict(model2$gam,newdat,type="link",se.fit=T)$se
 newdat$prev2<-model$family$linkinv(link)
 newdat$CIinf2<-model$family$linkinv(link-1.96*linkse)
 newdat$CIsup2<-model$family$linkinv(link+1.96*linkse)
 plot(newdat$prev2~newdat$study_median_age, type="l",col="red",ylim=c(0,0.11))
 lines(newdat$CIinf2~newdat$study_median_age, lty=2,col="red")
 lines(newdat$CIsup2~newdat$study_median_age, lty=2,col="red")
 lines(newdat$prev~newdat$study_median_age, type="l",ylim=c(0,.12))
 lines(newdat$CIinf~newdat$study_median_age, lty=2)
 lines(newdat$CIsup~newdat$study_median_age, lty=2)

由于研究级别的随机效应在适合的部分，我不必处理它。

正如你所看到的，我得到了相当不同的结果，年龄和患病率之间的关系更加平滑，置信区间也完全不同。在全数据分析中甚至更加不同，其中 CI 在模型中更广泛，包括研究级 RE，以至于有时几乎没有信息（流行率在 0 到 15% 之间，但如果是这样的话是...）。此外，当排除异常值时，研究级 RE 模型似乎更稳定。

所以，我的问题是：

我是否正确地从 metaprop() 函数中提取了权重并进一步使用它们？
我是否正确构建了我的模型gam()和gamm4()模型？我读了很多关于这个，但我不习惯这个模型之王。
我应该使用这些模型中的哪一个？

我真的很感激一些帮助，因为我的老师和我的同事都不能。进行系统审查真的很苛刻，而且在分析中挣扎非常令人沮丧......提前谢谢你！