我可以使用 mi 包汇集估算的随机效应模型估计值吗?

2024-02-12

看来mi在过去几年的某个时候,软件包进行了相当大的重写。

以下教程详细概述了“旧”的处理方式:http://thomasleeper.com/Rcourse/Tutorials/mi.html http://thomasleeper.com/Rcourse/Tutorials/mi.html

“新”的做事方式(坚持 Leeper 的模拟演示)看起来像这样:

#load mi
library(mi)
#set seed
set.seed(10)
#simulate some data (with some observations missing)
x1 <- runif(100, 0, 5)
x2 <- rnorm(100)
y <- 2*x1 + 20*x2 + rnorm(100)
mydf <- cbind.data.frame(x1, x2, y)
mydf$x1[sample(1:nrow(mydf), 20, FALSE)] <- NA
mydf$x2[sample(1:nrow(mydf), 10, FALSE)] <- NA

# Convert to a missing_data.frame
mydf_mdf <- missing_data.frame(mydf)

# impute
mydf_imp <- mi(mydf_mdf)

尽管函数名称已更改,但这实际上与“旧”的处理方式非常相似。

最大的变化(从我的角度来看)是替换了以下“旧”功能

lm.mi(formula, mi.object, ...)

glm.mi(formula, mi.object, family = gaussian, ...)

bayesglm.mi(formula, mi.object, family = gaussian, ...)

polr.mi(formula, mi.object, ...)

bayespolr.mi(formula, mi.object, ...)

lmer.mi(formula, mi.object, rescale=FALSE, ...)

glmer.mi(formula, mi.object, family = gaussian, rescale=FALSE, ...).

以前,用户可以使用这些函数之一为每个估算数据集计算模型,然后使用mi.pooled() (or coef.mi()如果我们遵循 Leeper 的例子)。

在当前版本中mi(我安装了 v1.0),最后这些步骤似乎已合并为一个函数,pool(). The pool()函数似乎读取在上述插补过程中分配给变量的族和链接函数,然后使用以下方法估计模型bayesglm使用如下所示的指定公式。

# run models on imputed data and pool the results
summary(pool(y ~ x1 + x2, mydf_imp))

## 
## Call:
## pool(formula = y ~ x1 + x2, data = mydf_imp)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.98754  -0.40923   0.03393   0.46734   2.13848  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.34711    0.25979  -1.336    0.215    
## x1           2.07806    0.08738  23.783 1.46e-13 ***
## x2          19.90544    0.11068 179.844  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.7896688)
## 
##     Null deviance: 38594.916  on 99  degrees of freedom
## Residual deviance:    76.598  on 97  degrees of freedom
## AIC: 264.74
## 
## Number of Fisher Scoring iterations: 7

看起来我们即将恢复模拟 Beta 值(2 和 20)。换句话说,它的行为符合预期。

让我们采用一组稍微大一些的数据,并简单地模拟随机效应,只是为了获得分组变量。

mydf2 <- data.frame(x1 = rep(runif(100, 0, 5), 20)
                   ,x2 = rep(rnorm(100, 0, 2.5), 20)
                   ,group_var = rep(1:20, each = 100)
                   ,noise = rep(rnorm(100), 20))

mydf2$y <- 2*mydf2$x1 + 20*mydf2$x2 + mydf2$noise

mydf2$x1[sample(1:nrow(mydf2), 200, FALSE)] <- NA
mydf2$x2[sample(1:nrow(mydf2), 100, FALSE)] <- NA

# Convert to a missing_data.frame
mydf2_mdf <- missing_data.frame(mydf2)

show(mydf2_mdf)

## Object of class missing_data.frame with 2000 observations on 5 variables
## 
## There are 4 missing data patterns
## 
## Append '@patterns' to this missing_data.frame to access the corresponding pattern for every observation or perhaps use table()
## 
##                 type missing method  model
## x1        continuous     200    ppd linear
## x2        continuous     100    ppd linear
## group_var continuous       0   <NA>   <NA>
## noise     continuous       0   <NA>   <NA>
## y         continuous       0   <NA>   <NA>
## 
##             family     link transformation
## x1        gaussian identity    standardize
## x2        gaussian identity    standardize
## group_var     <NA>     <NA>    standardize
## noise         <NA>     <NA>    standardize
## y             <NA>     <NA>    standardize

Since missing_data.frame()似乎正在翻译group_var作为连续的,我使用change()函数来自mi重新分配给"un"对于“无序分类”,然后按上述方式进行。

mydf2_mdf <- change(mydf2_mdf, y = "group_var", what = "type", to = "un"  )

# impute
mydf2_imp <- mi(mydf2_mdf)

现在,除非1.0版本mi删除了以前版本的功能(即可用的功能lmer.mi and glmer.mi),我假设在公式中添加随机效应应该指出pool()到适当的lme4功能。然而,最初的错误消息表明情况并非如此。

# run models on imputed data and pool the results
summary(pool(y ~ x1 + x2 + (1|group_var), mydf2_imp))
## Warning in Ops.factor(1, group_var): '|' not meaningful for factors
## Warning in Ops.factor(1, group_var): '|' not meaningful for factors
## Error in if (prior.scale[j] < min.prior.scale) {: missing value where TRUE/FALSE needed

按照我的警告消息并从我的因子中提取整数确实可以让我得到一个估计,但结果表明pool()仍在估计固定效应模型bayesglm并保持我尝试的随机效应不变。

summary(pool(y ~ x1 + x2 + (1|as.numeric(as.character(group_var))), mydf2_imp))

## 
## Call:
## pool(formula = y ~ x1 + x2 + (1 | as.numeric(as.character(group_var))), 
##     data = mydf2_imp)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.93633  -0.69923   0.01073   0.56752   2.12167  
## 
## Coefficients:
##                                               Estimate Std. Error  t value
## (Intercept)                                  1.383e-01  2.596e+02    0.001
## x1                                           1.995e+00  1.463e-02  136.288
## x2                                           2.000e+01  8.004e-03 2499.077
## 1 | as.numeric(as.character(group_var))TRUE -3.105e-08  2.596e+02    0.000
##                                             Pr(>|t|)    
## (Intercept)                                        1    
## x1                                            <2e-16 ***
## x2                                            <2e-16 ***
## 1 | as.numeric(as.character(group_var))TRUE        1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.8586836)
## 
##     Null deviance: 5384205.2  on 1999  degrees of freedom
## Residual deviance:    1713.9  on 1996  degrees of freedom
## AIC: 5377
## 
## Number of Fisher Scoring iterations: 4

我的问题是:

  1. 是否可以使用以下方法轻松生成汇总随机效应估计mi?, and
  2. 如果是,怎么办?

只是为了提供一种替代方案,有一个包非常关注混合效应模型的 MI 以及汇集从中获得的结果(mitml, 在这里找到它 https://cran.r-project.org/web/packages/mitml/index.html).

使用该包非常简单。它依赖于包pan and jomo用于插补,但它也可以处理来自其他 MI 包的输入(?as.mitml.list).

来自混合效应模型的汇集估计大部分是自动化的,并包含在testEstimates功能。

require(mitml)
require(lme4)

data(studentratings)

# impute example data using 'pan'
fml <- ReadDis + SES ~ ReadAchiev + (1|ID)
imp <- panImpute(studentratings, formula=fml, n.burn=1000, n.iter=100, m=5)

implist <- mitmlComplete(imp, print=1:5)

# fit model using lme4
fit.lmer <- with(implist, lmer(SES ~ (1|ID)))

# pool results using 'Rubin's rules'
testEstimates(fit.lmer, var.comp=TRUE)

Output:

# Call:

# testEstimates(model = fit.lmer, var.comp = TRUE)

# Final parameter estimates and inferences obtained from 5 imputed data sets.

#              Estimate Std.Error   t.value        df   p.value       RIV       FMI 
# (Intercept)    46.988     1.119    41.997   801.800     0.000     0.076     0.073 

#                         Estimate 
# Intercept~~Intercept|ID   38.272 
# Residual~~Residual       298.446 
# ICC|ID                     0.114 

# Unadjusted hypothesis test as appropriate in larger samples. 
本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容,请联系:hwhale#tublm.com(使用前将#替换为@)

我可以使用 mi 包汇集估算的随机效应模型估计值吗? 的相关文章

  • 使用变量在 r 中像 aes_string 一样选择轴

    我试图提供一个带有列名的变量来创建一个plotly图表 类似于ggplot2 aes string 不知怎的 我被困住了 plot ly iris x Sepal Length y Sepal Width works as expected
  • RStudio 在临时目录中从 Rmarkdown 创建 PDF 文件

    我使用 RStudio 和 Rmarkdown 来创建报告 由于某种原因 使用 Knit 按钮会导致它仅在某个临时目录中创建 从命令输出来看 似乎 pandoc 本身被指示这样做 我觉得这很奇怪 usr lib rstudio bin pa
  • 将 data.frame 的列中的值替换为另一个 data.frame 中的值

    我的情况是 我有一个数据框 其中有一列填充了整数 1 到 6 我想用更具描述性的标签替换这些整数 这些标签在另一个充当 键 的数据框中提供 V1 V2 1 1 LABEL1 2 2 LABEL2 3 3 LABEL3 4 4 LABEL4
  • 在 mts 对象上使用 Apply 系列函数

    在 mts 对象上使用 apply 或 sapply 会在发送到函数时删除其时间序列属性 我应该如何在 mts 对象中的每个时间序列上应用相同的函数 带有 ts 输入和 ts 输出 并返回它 最好是 mts 我的意思是除了使用 for 循环
  • 单击并按住 R 中的按钮闪亮?

    我希望能够通过单击 R 闪亮按钮来更改参数的值 所以我需要按钮 一个用于增加值 一个用于减少值 我想在按住按钮的同时保持值以一定的速度减少 增加 通过释放按钮的点击 动作应该停止 到目前为止我还没有找到这个选项actionButtons在
  • 在 for 循环中绘制的多个 ggplot2 绘图的网格

    作为一个新的 ggplot2 用户 我对可能性的数量感到有点迷失 并且很难在网上找到我认为简单问题的简单答案 我想在同一张纸上显示 ggplot2 的多个图 但知道这些图来自 for 循环 以下示例无法编译 仅用于说明 for i in c
  • 在 Shiny 中使用 readlines(prompt = )

    我有一个代码 使用以下方式获取输入readlines prompt 功能 您能告诉我 Shiny 中的哪个输入函数足以将此代码适应 Shiny 应用程序吗 我需要一个交互功能 我无法使用简单的输入selectInput 因为我有很多read
  • glm() 模型的交叉验证

    我正在尝试对我之前在 R 中构建的一些 glm 模型进行 10 倍交叉验证 我对cv glm 函数在boot包 尽管我已经阅读了很多帮助文件 当我提供以下公式时 library boot cv glm data glmfit K 10 这里
  • udunits2 R 安装:找不到 udunits2.h

    我正在尝试在 R 中安装 udunits2 以满足对ggforce包裹 但是 安装程序在检查 udunits2 时始终失败 我已经尝试过中的说明this https stackoverflow com questions 47059517
  • 仅在具有重复块名称的另一个 Rmarkdown 文档中运行一个 Rmarkdown 文档中的代码

    我正在 Rmarkdown 中编写一系列相互补充的报告 我想将上一份报告的结果纳入我目前正在编写的报告中 我看到其他建议使用的问题purl从 Rmarkdown 文档中提取 R 代码然后运行它 所以我尝试了以下操作 r read previ
  • mclapply 用户时间大于已用时间

    我正在尝试使用mclapply的功能parallel封装在R 该函数通过计算对数似然距离将值分配给序列矩阵 这是一个 CPU 密集型操作 所结果的system time价值观令人困惑 gt system time mclapply work
  • 如何用日语创建 ggplot2 标题?

    我正在准备日语演示文稿 并希望图像的标题和图例名称为日语 我可以让文本在 RStudio 中渲染得很好 但是当渲染图像时 日语字符仅显示为方框 x 10 10 y x x df data frame x y ggplot df aes x
  • 我可以调整scale_color_brewer的下限吗?

    我已经订购了我想使用 color Brewer 的分类数据 但我很难看到非常低的值 有没有办法去掉这些较低的值或设置范围的下限 ggplot data frame x 1 6 y 10 15 w letters 1 6 aes x y co
  • 当子集长度为零时,如何简洁地处理子集?

    从向量中排除元素x x lt c 1 4 3 2 我们可以减去位置向量 excl lt c 2 3 x excl 1 1 2 这也是动态工作的 excl lt which x which max x gt quantile x 25 1 2
  • 如何将 mcmc.list 转换为 bugs 对象?

    我正在使用rjagsR 库 功能coda samples产生一个mcmc list 例如 来自example coda samples library rjags data LINE LINE recompile LINE out lt c
  • 通过 RCpp 返回 NA

    新手 RCpp 问题在这里 How can I make a NumericVector returnNA到R 例如 假设我有一个 RCpp 代码 它分配NA到向量的第一个元素 RCpp export NumericVector myFun
  • 使用 dplyr 的 select 引用变量名[重复]

    这个问题在这里已经有答案了 通常我会想要选择变量的子集 其中该子集是函数的结果 在这个简单的例子中 我首先获取与宽度特征相关的所有变量名称 library dplyr library magrittr data iris width var
  • 在 ggplot 中过滤管道 df

    我正在使用 dplyr 管道来清理我的 df 然后直接输入到 ggplot 中 但是 我只想一次只绘制一组 因此我需要过滤到该组 问题是 我希望比例保持不变 就好像所有群体都存在一样 是否可以在 ggplot 命令中进一步过滤管道 df 例
  • 通过 r 中的组变量进行汇总

    我有一个数据框如下 head newStormObject FATALITIES INJURIES PROPVALDMG CROPVALDMG EVTYPE total 1 0 15 2 5e 05 0 TORNADO 15 2 0 0 2
  • R - 如何为数据范围内的缺失值绘制条形图零点?

    假设我对 1 到 10 之间的整数的 200 个点有 10 个观察值 mysample sample rep seq 1 10 20 10 我想用条形图绘制它 barplot table mysample barplot https i s

随机推荐