[R-sig-ME] Simulating mixed-effects longitudinal logistic data in R
Uanhoro, James
u@nhoro@1 @end|ng |rom buckeyem@||@o@u@edu
Mon Jul 13 00:25:26 CEST 2020
You can keep the same exact code with some changes:
- If you want to use logistic regression to recover the population
parameters, set: eij <- rlogis(J * n.time, scale = 1)
- If you use the normal distribution (rnorm), you will need probit
regression to recover population parameters: eij <- rnorm(J * n.time,
sd = 1)
- Note that the sd/scale has to be 1 as the estimated parameters will
be re-scaled by the value of sd(eij)
- Your y then will have to be as below to create a binary variable: y
<- ((rowSums(X * betaj[st.id, ]) + eij) > 0) + 0. Basically, you
dichotomize the generated y at 0. If you dichotomize at some other
value, the estimated intercept will shift accordingly.
- Finally, your population parameters have to make sense of the logit
scale. Intercept values of 300 and coefficients values of 25 are
unbelievable log-odds. You will generate all 1s with these population
parameters.
This is the latent variable approach to setting up the binary outcome
model.
The other answers do not use an error term to set up the model. For
these responses, you generate probabilities using either pnorm(LP)
(probit) or plogis(LP) (logistic) where LP is the linear predictor,
then pass the probabilities to rbinom() to obtain a binary variable.
James.
On Sun, 2020-07-12 at 17:11 -0500, Simon Harmel wrote:
> 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://urldefense.com/v3/__https://drizopoulos.github.io/GLMMadaptive/articles/GLMMadaptive_basics.html__;!!KGKeukY!jq0K6FHI0BtJ8LoX0xkhk9pv-4RgkiFLMeqSS8izqfz60XgQKODdgFKHcwYytJ5tqjwbaith3Ks$
> >
> > >
> > > 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://urldefense.com/v3/__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__;JSUlJSUlJSUlJSUlJSUlJSUlJQ!!KGKeukY!jq0K6FHI0BtJ8LoX0xkhk9pv-4RgkiFLMeqSS8izqfz60XgQKODdgFKHcwYytJ5tqjwbWhOO65s$
> >
> > >
> > > [[alternative HTML version deleted]]
> > >
> > >
> > > _______________________________________________
> > > R-sig-mixed-models using r-project.org mailing list
> > >
https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models__;!!KGKeukY!jq0K6FHI0BtJ8LoX0xkhk9pv-4RgkiFLMeqSS8izqfz60XgQKODdgFKHcwYytJ5tqjwbOi_tX1Q$
> > >
> >
> > [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-mixed-models using r-project.org mailing list
> >
https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models__;!!KGKeukY!jq0K6FHI0BtJ8LoX0xkhk9pv-4RgkiFLMeqSS8izqfz60XgQKODdgFKHcwYytJ5tqjwbOi_tX1Q$
> >
> >
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
>
https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models__;!!KGKeukY!jq0K6FHI0BtJ8LoX0xkhk9pv-4RgkiFLMeqSS8izqfz60XgQKODdgFKHcwYytJ5tqjwbOi_tX1Q$
>
More information about the R-sig-mixed-models
mailing list