从 lmer 获取方差分析表应用于 R 中多重插补后的数据

Question

我一直在缺失值的数据集上运行 lmer。我已经成功地估算了缺失值并成功运行了我的模型并获得了汇总拟合的摘要。但是，我需要用于报告的方差分析表和用于绘图的估计方法。我一直在使用的代码如下，是我作为一个新手从各种来源拼凑而成的。如果您发现其中有问题，请告诉我。

# remove unneeded columns before imputation
columns_to_remove = c(1,3,15:19,21:30)
dataset_with_missing_values <- dataset_with_missing_values[, -columns_to_remove] 

md.pattern(dataset_with_missing_values) # show pattern of missing data

dat <- dataset_with_missing_values
original <- dat
impmethod <- character(ncol(dat))

init <- mice(dat, maxit = 0)
meth <- init$method
predM <- init$predictorMatrix
predM["Subject", ] <- 0
impmethod["dependent_variable"] <- "2l.lmer"

set.seed(103)
imputed <- mice(dataset_with_missing_values, method = meth, predictorMatrix = predM, m = 10)
complete_imputed <- complete(imputed)

sapply(original, function(x) sum(is.na(x))) # check to see missing values have been imputed
sapply(complete_imputed, function(x) sum(is.na(x)))

mean(original$dependent_variable, na.rm = TRUE) # compare means re
mean(complete_imputed$dependent_variable, na.rm = TRUE)

fit <- with(data = imputed, 
exp = lme4::lmer(dependent_variable ~ Factor_1 * Factor_2 * Sex * Time + (1 | Subject)))

pooled_fit <- pool(fit)

plot_summs(pooled_fit) # get summary of fit

summary <- summary(pooled_fit)`

我尝试过 anova(pooled_fit) 但出现错误：

没有适用于“anova”的方法应用于“c('mipo', 'data.frame')”类的对象

aov(pooled_fit) 给出：

x$terms %||% attr(x, "term") %||% stop("没有术语组件也没有属性") 中的错误： 没有术语组件或属性

Answer 1

如果您需要提取摘要数据，例如为了绘图，你可以使用 broom 包 tidiers 来处理

lmer

和

anova

对象：

（如果您可以在代码中包含示例数据，那就更好了，所以我使用了 lme4 包中的玩具数据。）

library(lme4)
library(broom.mixed)

fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)

tidy(fm1)
#> # A tibble: 6 × 6
#>   effect   group    term                  estimate std.error statistic
#>   <chr>    <chr>    <chr>                    <dbl>     <dbl>     <dbl>
#> 1 fixed    <NA>     (Intercept)           251.          6.82     36.8 
#> 2 fixed    <NA>     Days                   10.5         1.55      6.77
#> 3 ran_pars Subject  sd__(Intercept)        24.7        NA        NA   
#> 4 ran_pars Subject  cor__(Intercept).Days   0.0656     NA        NA   
#> 5 ran_pars Subject  sd__Days                5.92       NA        NA   
#> 6 ran_pars Residual sd__Observation        25.6        NA        NA

glance(fm1)
#> # A tibble: 1 × 7
#>    nobs sigma logLik   AIC   BIC REMLcrit df.residual
#>   <int> <dbl>  <dbl> <dbl> <dbl>    <dbl>       <int>
#> 1   180  25.6  -872. 1756. 1775.    1744.         174

augment(fm1)
#> # A tibble: 180 × 14
#>    Reaction  Days Subject .fitted  .resid   .hat .cooksd .fixed   .mu .offset
#>       <dbl> <dbl> <fct>     <dbl>   <dbl>  <dbl>   <dbl>  <dbl> <dbl>   <dbl>
#>  1     250.     0 308        254.   -4.10 0.229  0.00496   251.  254.       0
#>  2     259.     1 308        273.  -14.6  0.170  0.0402    262.  273.       0
#>  3     251.     2 308        293.  -42.2  0.127  0.226     272.  293.       0
#>  4     321.     3 308        313.    8.78 0.101  0.00731   283.  313.       0
#>  5     357.     4 308        332.   24.5  0.0910 0.0506    293.  332.       0
#>  6     415.     5 308        352.   62.7  0.0981 0.362     304.  352.       0
#>  7     382.     6 308        372.   10.5  0.122  0.0134    314.  372.       0
#>  8     290.     7 308        391. -101.   0.162  1.81      325.  391.       0
#>  9     431.     8 308        411.   19.6  0.219  0.106     335.  411.       0
#> 10     466.     9 308        431.   35.7  0.293  0.571     346.  431.       0
#> # ℹ 170 more rows
#> # ℹ 4 more variables: .sqrtXwt <dbl>, .sqrtrwt <dbl>, .weights <dbl>,
#> #   .wtres <dbl>

at1 <- anova(lmer(Reaction ~ Days + (Days | Subject), sleepstudy))

tidy(at1)
#> # A tibble: 1 × 5
#>   term   npar  sumsq meansq statistic
#>   <chr> <int>  <dbl>  <dbl>     <dbl>
#> 1 Days      1 30031. 30031.      45.9

^{创建于 2024-05-03，使用 reprex v2.1.0}