[R-sig-ME] Simulating mixed-effects longitudinal logistic data in R

[Simon Harmel]

Many thanks all!  I basically want to generate 100 bernoulli outcomes whose
probability of being 1 is longitudinally modeled based on 4 time points
(0:3), membership to 2 groups (Control vs. treatment), and gender (male vs.
female).

Thanks a lot!

On Sun, Jul 12, 2020 at 5:02 PM Ben Bolker <bbolker using gmail.com> wrote:

>    Or:
>
>
> (1)  see ?lme4::simulate.merMod(), i.e. set up your covariates, then
> simulate(formula, newdata=data_frame, newparams=..., family=binomial,
> weights=1)
>
> or
>
> (2) after setting everything up as in your example, use
>
> response <- rbinom(length(y), prob=plogis(y), size=1)
>
> On 7/12/20 5:58 PM, D. Rizopoulos wrote:
> > Have a look in this example:
> https://drizopoulos.github.io/GLMMadaptive/articles/GLMMadaptive_basics.html
> >
> > Best,
> > Dimitris
> >
> > ��
> > Dimitris Rizopoulos
> > Professor of Biostatistics
> > Erasmus University Medical Center
> > The Netherlands
> > ________________________________
> > From: R-sig-mixed-models <r-sig-mixed-models-bounces using r-project.org> on
> behalf of Simon Harmel <sim.harmel using gmail.com>
> > Sent: Sunday, July 12, 2020 10:44:17 PM
> > To: r-sig-mixed-models <r-sig-mixed-models using r-project.org>
> > Subject: [R-sig-ME] Simulating mixed-effects longitudinal logistic data
> in R
> >
> > Good afternoon,
> >
> > I know to simulate a linear mixed-effects longitudinal data in R, I can
> use
> > the procedure below.
> >
> > But I wonder how to generate a similar design, but for a logistic model
> in
> > R?
> >
> > I appreciate your suggestions,
> > Simon
> >
> > #------------- linear mixed-effects longitudinal data in R-------
> >
> > library(MASS)
> >
> > growth.sim <- function(J, n.time, gammas, G, sigma = 1, seed = NULL) {
> >
> >     X <- cbind(1, seq_len(n.time) - 1) # time indicators for each
> individual
> >     X <- X[rep(seq_len(n.time), J), ]  # Repeat each row n.time times
> >
> >    st.id <- seq_len(J)                 # individual id
> >    st.id <- rep(st.id, each = n.time)  # repeat each ID n.time times
> >
> >
> > set.seed(seed) ## what should go below, possibly correlated Beta?
> > uj <- MASS::mvrnorm(J, mu = rep(0, 2), Sigma = G) # Generate u0 and u1
> > random effects for intercepts and slopes
> >
> > ## I think no error term is needed for logistic version?
> > eij <- rnorm(J * n.time, sd = sigma)              # Generate error
> > term for observations
> >
> >
> > betaj <- matrix(gammas, nrow = J, ncol = 2, byrow = TRUE) + uj #
> > Compute beta_j's
> >
> > y <- rowSums(X * betaj[st.id, ]) + eij          # Compute outcome
> > based on Yij = sum(Xij*Bj) + eij:
> >
> > dat <- data.frame(st.id, time = X[ , 2], y) # Output a data frame
> > return(dat)                                 # Return data}
> > *##==== Example of use: =====*
> > growth.sim(10, 4, gammas = c(300, 25),
> >         G = matrix(c(0.1, 0,
> >                      0, 0.01), nrow = 2))
> >
> >          [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-mixed-models using r-project.org mailing list
> >
> https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-mixed-models&data=02%7C01%7Cd.rizopoulos%40erasmusmc.nl%7Cf46c256b636848749fc608d826a46b4e%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C0%7C637301834909767585&sdata=GjOGQ%2Fv5aMwH8PK%2FuM%2FEeHkElDFwB6%2BaD8UZy2wVdx4%3D&reserved=0
> >
> >       [[alternative HTML version deleted]]
> >
> >
> > _______________________________________________
> > R-sig-mixed-models using r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>

	[[alternative HTML version deleted]]