在 R 中模拟多级模型 (lme/lmer) 的能力：broom::tidy

问题描述

去年八月左右，我编写了以下代码来模拟多级模型的功率。 它工作了很长时间，但现在似乎 broom::tidy() 已被多级模型弃用。 我不明白这究竟如何影响我的代码，但这是发生了什么：

#Power simulation for mixed effects models

library(nlme)
library(dplyr)
library(tidyr)
library(purrr)

# 1. First,determine your parameters. The below numbers will create 4 dfs with different Ns 
#(eg N=5,10,50,& 200) and cross each level with the other level of the other parameters specified.

params <- expand.grid(
  N = c(50,135,1000),#this is the N for each sample df
  B = c(.1,.3,.5),#the effect size B,a linear coefficient of x predicting y
  items = c(5,20)) #the number of repeating within-subject time points,items,etc
n_samples <- 100 #how many simulations you want to run

#2. Then,run this entire chunk with cmd + shift + enter

#Don't modify anything below
###############################################################
#create a sample
sim_sample <- function(N,B,items){
tid <- rep(1:items,N) #time id,not really needed but maybe useful
subj <- rep(1:N,each=items) #subject id
x = sample(1:6,items*N,replace = T) #a fixed effect predictor
y <- B*x + rnorm(n = items*N) #y is defined as being predicted by x + random error
smaller.df <- data.frame(tid,subj,x,y) #combine them into a df

lmeout <- summary(lme(y ~ x,random=~1|subj,data=smaller.df))$tTable %>% 
  data.frame()
broom::tidy(lmeout) #build a model,extract coefs,put coefs into a table
}

#another function to simulate all the parameters specified above
sim_all_params <- function(params) {
  sim_results <- purrr::pmap(params,sim_sample)
  params %>%
    mutate(sim_results = sim_results) %>%
    tidyr::unnest(sim_results)
}

#execute the simulation and re-run it the pre-specified amount of times
#sim_all_params(params)
full_results <- purrr::rerun(.n = n_samples,sim_all_params(params)) %>%
  bind_rows(.id = 'sample_num') %>% 
  as.data.frame(full_results)
#save the full results in a table
table <- full_results %>%
  #dplyr::filter(rownames == 'x') %>% #this select only the rows for the x predictor in the case that you have multiple predictors
  dplyr::group_by(N,items) %>%
  dplyr::summarise(power = sum(p.value < .05)/n()) #compute power for each row
#print the table
table

模拟代码不再适用于 lme() 功率估计值始终在约 50% 左右；如果我增加样本量，功效不会增加 我尝试使用 tibble() 和 data.frame() 而不是 broom::tidy() 但我得到了奇怪的结果，如 power = 1.5 或 2.5。

我真的不明白交易是什么以及为什么 broom::tidy() 的这种变化会对这段代码产生如此大的影响。

另一方面，将以下代码用于非多级设计（例如 t 检验）时，功效模拟可以正常工作：

#Power simulation for t test

library(nlme)
library(dplyr)
library(tidyr)
library(purrr)

# 1. First,& 200) and cross each level with the other level of the other parameters specified.

params <- expand.grid(
  N = c(20,75,150,300),#this is the n in *each* group
  d = c(.2,.4,#standardised mean difference (ie Cohen's d)
    int = c(1,1,1) #this is the intercept (ie mean)
)
n_samples <- 100 #how many simulations you want to run
#2. Then,run this entire chunk with cmd + shift + enter


#Don't modify anything below
###############################################################
#create a sample
sim_sample <- function(N,d,int){
d1 <- data.frame(x = rnorm(N,int,1),condition = 0)
d2 <- data.frame(x = rnorm(N,(int+d),condition = 1)
smaller.df <- rbind(d1,d2)

tout <- t <- broom::tidy(t.test(x ~ condition,smaller.df)) #build a model,extract coefs
#broom::tidy(tout) #put coefs into a table
}

#another function to simulate all the parameters specified above
sim_all_params <- function(params) {
  sim_results <- purrr::pmap(params,sim_sample)
  
  params %>%
    mutate(sim_results = sim_results) %>%
    tidyr::unnest(sim_results)
}

#execute the simulation and re-run it the pre-specified amount of times
#sim_all_params(params)
full_results <- purrr::rerun(.n = n_samples,sim_all_params(params)) %>%
  bind_rows(.id = 'sample_num')
#save the full results in a table
table <- full_results %>%
  group_by(N,d) %>%
  summarise(power = sum(p.value < .05)/n()) #compute power for each row
#print the table
table

谁能帮帮我

1. 修复第一个代码
2. 解释为什么它不再有效？

谢谢！

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）