[R-sig-ME] Simulating Data from a Generalized Linear Mixed Model

[Dimitris Rizopoulos]

for simple models you can use something like the following:

sleepstudy$Reaction2 <- sleepstudy$Reaction > 
median(sleepstudy$Reaction)
(fm1 <- lmer(Reaction2 ~ Days + (Days | Subject), family = "binomial", 
sleepstudy))
############
n <- length(unique(sleepstudy$Subject))
ni <- as.vector(tapply(sleepstudy$Subject, sleepstudy$Subject, 
length))
p <- fm1 at nc
X <- fm1 at X
Z <- X[, 1:2] # keep the relevant columns of X
V <- matrix(VarCorr(fm1)[[1]]@x, p)
b <- mvrnorm(n, rep(0, p), V)
mu <- X %*% fixef(fm1) + rowSums(Z * b[rep(1:n, ni), ])
y <- rbinom(sum(ni), 1, plogis(mu))
# y <- rnorm(sum(ni), mu, attr(VarCorr(fm1), "sc")) for linear mixed 
model


I hope it helps.

Best,
Dimitris

----
Dimitris Rizopoulos
Ph.D. Student
Biostatistical Centre
School of Public Health
Catholic University of Leuven

Address: Kapucijnenvoer 35, Leuven, Belgium
Tel: +32/(0)16/336899
Fax: +32/(0)16/337015
Web: http://med.kuleuven.be/biostat/
     http://www.student.kuleuven.be/~m0390867/dimitris.htm


----- Original Message ----- 
From: "Chuck Cleland" <ccleland at optonline.net>
To: <r-sig-mixed-models at r-project.org>
Sent: Friday, August 24, 2007 2:42 PM
Subject: [R-sig-ME] Simulating Data from a Generalized Linear Mixed 
Model


> Hello:
>  Does anyone have suggestions for how one can simulate data for a
> generalized linear mixed model?  For example, below are summaries of 
> two
> models fit by lmer() for which I would like to simulate new data. 
> If
> the summaries do not provide sufficient information to generate new
> observations, what other information is needed?
>
> thanks for any help,
>
> Chuck
>
>> summary(fm1)
> Linear mixed-effects model fit by REML
> Formula: Y ~ V7 * TIME + (TIME | id)
>   Data: df.uni
>  AIC  BIC logLik MLdeviance REMLdeviance
> 8795 8836  -4390       8763         8781
> Random effects:
> Groups   Name Variance Std.Dev. Corr
> id            0.23671  0.48653
>               0.52332  0.72341  -0.079
> Residual      0.51130  0.71506
> number of obs: 2611, groups: id, 480
>
> Fixed effects:
>            Estimate Std. Error t value
> (Intercept)  0.02499    0.04311   0.580
> V7          -0.08775    0.06027  -1.456
> TIME         0.15238    0.04826   3.158
> V7:TIME      0.29308    0.06717   4.364
>
> Correlation of Fixed Effects:
>        (Intr) V7     TIME
> V7      -0.715
> TIME    -0.106  0.076
> V7:TIME  0.076 -0.107 -0.718
>
>> summary(fm2)
> Generalized linear mixed model fit using Laplace
> Formula: Y ~ V7 * TIME + (TIME | id)
>   Data: df.uni
> Family: binomial(logit link)
>  AIC  BIC logLik deviance
> 2278 2319  -1132     2264
> Random effects:
> Groups Name Variance Std.Dev. Corr
> id          0.63468  0.79667
>             2.96036  1.72057  0.125
> number of obs: 2601, groups: id, 480
>
> Estimated scale (compare to  1 )  0.6274443
>
> Fixed effects:
>            Estimate Std. Error z value Pr(>|z|)
> (Intercept)  0.29549    0.12484   2.367   0.0179 *
> V7           0.08769    0.17198   0.510   0.6101
> TIME        -0.67834    0.14820  -4.577 4.71e-06 ***
> V7:TIME      0.45752    0.20165   2.269   0.0233 *
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Correlation of Fixed Effects:
>        (Intr) V7     TIME
> V7      -0.726
> TIME    -0.260  0.189
> V7:TIME  0.191 -0.253 -0.735
>
> -- 
> Chuck Cleland, Ph.D.
> NDRI, Inc.
> 71 West 23rd Street, 8th floor
> New York, NY 10010
> tel: (212) 845-4495 (Tu, Th)
> tel: (732) 512-0171 (M, W, F)
> fax: (917) 438-0894
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> 


Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm