假设,我有一个数据集,类似于
require(nlme)
?Orthodont
我的模型是
fm2 <- lme(distance ~ age + Sex, data = Orthodont, random = ~ 1)
如何使用模型拟合对象fm2生成多个数据集,样本大小为 300、400、500,...?
我在 r-sig-mixed-models 帮助上阅读了这个很好的答案,但它似乎不完整。
假设,我有一个数据集,类似于
require(nlme)
?Orthodont
我的模型是
fm2 <- lme(distance ~ age + Sex, data = Orthodont, random = ~ 1)
如何使用模型拟合对象fm2生成多个数据集,样本大小为 300、400、500,...?
我在 r-sig-mixed-models 帮助上阅读了这个很好的答案,但它似乎不完整。
注意:使用simulate.lme 的模拟数据不匹配原始数据结构或模型拟合的元素(例如方差、效应大小......),也不为实验设计测试从头创建数据。
require(nlme)
?nlme::simulate.lme
fit <- lme(distance ~ age + Sex, data = Orthodont, random = ~ 1)
orthSim <- simulate.lme(fit, nsim = 1)
这会产生一个模拟拟合(使用可能的替代模型)。
这要归功于@Momo 对我的一个问题的回答: 是否有一种通用的方法可以从公式或分析中模拟数据?
如果您需要模拟数据,则需要从模拟.lme 函数创建一个新函数。
simulate.lme.data<-edit(simulate.lme)
在最后一个括号之前添加以下行
return(base2)
然后,您可以根据需要创建尽可能多的数据:
orthSimdata <- simulate.lme.data(fit, nsim = 1)
请注意,这是来自我对模拟.lme 中未注释代码的(可能是错误的)解释。
尽管这很有用,但这似乎只是在现有数据中添加高斯噪声。
这不能用于直接从头模拟数据。我目前通过添加我的实验设计数据框的因子水平的数值来创建起始数据(例如 response=as.numeric(factor1)+as.numeric(factor2)+as.numeric(factor1)*as.numeric (因子 1)+rnorm(sd=2)...)。
这是一种从 中获取所有值的方法fm2。您可以向函数添加更多参数,以允许您更改模拟中的值。
library(nlme)
fm2 <- lme(distance ~ age + Sex, data = Orthodont, random = ~ 1)
simfun <- function(n) {
# n is the number of subjects, total rows will be 4*n
# sig.b0 is the st. dev. of the random intercepts
# it might be easier to just copy from the output
sig.b0 <- exp(unlist(fm2$modelStruct$reStruct))*fm2$sigma
b0 <- rnorm(n, 0, sig.b0)
sex <- rbinom(n, 1, 0.5) # assign sex at random
fe <- fixef(fm2)
my.df <- data.frame( Subject=rep(1:n, each=4),
int = fe[1] + rep(b0, each=4),
Sex=rep(sex,each=4), age=rep( c(8,10,12,14), n ) )
my.df$distance <- my.df$int + fe[2] * my.df$age +
fe[3]*my.df$Sex + rnorm(n*4, 0, fm2$sigma)
my.df$int <- NULL
my.df$Sex <- factor( my.df$Sex, levels=0:1,
labels=c('Male','Female') )
my.df
}
Orthodont2 <- simfun(100)
我可能会随机抽样并从您的数据中的主题中替换,直到我有正确的样本量。这是引导方法。它比识别变量的多元分布然后从中采样更简单。此外,bootstrap 不会对数据的多元结构做出额外的假设。
首先在你的大型模拟研究中设置参与者的数量
nits=300
在小型研究中获得独特的参与者
sub=unique(Orthodont$Subject)
随机抽取唯一的参与者并替换
subs=sample(sub,nits,rep=T)
制作一个空的数据框
df=Orthodont[-(1:dim(Orthodont)[1]),]
循环遍历样本大小并将其绑定在一起。
for( i in 1:nits) {
df=rbind(df,Orthodont[which(Orthodont$Subject==subs[i]),])
}
最后一个 for 循环很慢,有一种更好的编写方式。
现在您可以在更大的数据集上运行回归并观察置信区间变小。