Dear R users, I'm pretty new on using lmer package. My response is binary and I have fixed treatment effect (2 treatments) and random center effect (7 centers). I want to test the effect of treatment by fitting 2 models: Model 1: center effect (random) only Model 2: trt (fixed) + center (random) + trt*center interaction. Then, I want to compare these 2 models with Likelihood Ratio Test. Here are my lmer codes that I don't feel comfortable about their correctness. model1 <- try(lmer(cbind( yvect, nvect-yvect) ~ 1 + (1 | center), family = binomial, niter = 25, method = "Laplace", control = list(usePQL = FALSE) )) model2 <- try(lmer(cbind( yvect, nvect-yvect) ~ trt*center + ( 1 | center) , family = binomial, niter = 25, method = "Laplace", control = list(usePQL = FALSE) )) (I have attached outputs below) What I don't understand is; I thought in model2 I have defined center effect as a random effect. Then how come I got center effects and trt*center interactions under the fixed effects list on the output? Probably I didn't understand how to set up these models in lmer. Could anyone help me about this? Thanks a lot for your help... Emine model1 <- try(lmer(cbind( yvect, nvect-yvect) ~ 1 + (1 | center), family = binomial, niter = 25, method = "Laplace", control = list(usePQL = FALSE) ))>summary(model1)Generalized linear mixed model fit using Laplace Formula: cbind(yvect, nvect - yvect) ~ 1 + (1 | center) Family: binomial(logit link) AIC BIC logLik deviance 236.817 238.0951 -116.4085 232.817 Random effects: Groups Name Variance Std.Dev. center (Intercept) 0.088127 0.29686 number of obs: 14, groups: center, 7 Estimated scale (compare to 1) 0.2672612 Fixed effects: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.32084 0.14709 -2.1812 0.02917 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ################## model2 <- try(lmer(cbind( yvect, nvect-yvect) ~ trt*center + ( 1 | center) , family = binomial, niter = 25, method = "Laplace", control = list(usePQL = FALSE) ))>summary(model2)Generalized linear mixed model fit using Laplace Formula: cbind(yvect, nvect - yvect) ~ trt * center + (1 | center) Family: binomial(logit link) AIC BIC logLik deviance 30 39.58586 -1.547024e-07 3.094048e-07 Random effects: Groups Name Variance Std.Dev. center (Intercept) 5e-10 2.2361e-05 number of obs: 14, groups: center, 7 Estimated scale (compare to 1) 0.2672612 Fixed effects: Estimate Std. Error z value Pr(>|z|) (Intercept) -1.060869 0.065372 -16.2282 < 2e-16 *** trt2 1.118029 0.086842 12.8743 < 2e-16 *** center2 -0.325428 0.504256 -0.6454 0.51869 center3 -0.325440 0.504258 -0.6454 0.51868 center4 0.655407 0.413449 1.5852 0.11292 center5 -0.325433 0.504256 -0.6454 0.51869 center6 -0.325421 0.504255 -0.6454 0.51870 center7 0.655422 0.413448 1.5853 0.11291 trt2:center2 0.673737 0.651313 1.0344 0.30094 trt2:center3 -0.137183 0.651314 -0.2106 0.83318 trt2:center4 -0.307083 0.583845 -0.5260 0.59891 trt2:center5 -0.137203 0.651314 -0.2107 0.83316 trt2:center6 1.654555 0.712419 2.3224 0.02021 * trt2:center7 0.673705 0.651311 1.0344 0.30096 --- _________________________________________________________________ Office Live http://clk.atdmt.com/MRT/go/aub0540003042mrt/direct/01/
HI Emine, The parentheses specify a variable as random or fixed. You are using center for random intercepts, but fixed for interactions. You may want> mer(cbind( yvect, nvect-yvect) ~ trt + ( 1 | center) + (1 | > trt:center), > family = binomial, niter = 25, method = "Laplace", control = list > (usePQL > FALSE) )Cheers, Hank On Jun 1, 2007, at 2:01 PM, emine ?zg?r Bayman wrote:> Dear R users, > I'm pretty new on using lmer package. My response is binary and I > have fixed > treatment effect (2 treatments) and random center effect (7 > centers). I want > to test the effect of treatment by fitting 2 models: > > Model 1: center effect (random) only > Model 2: trt (fixed) + center (random) + trt*center interaction. > > Then, I want to compare these 2 models with Likelihood Ratio Test. > Here are > my lmer codes that I don't feel comfortable about their correctness. > > model1 <- try(lmer(cbind( yvect, nvect-yvect) ~ 1 + (1 | center), > family = binomial, niter = 25, method = "Laplace", control = list > (usePQL > FALSE) )) > > model2 <- try(lmer(cbind( yvect, nvect-yvect) ~ trt*center + ( 1 | > center) > , > family = binomial, niter = 25, method = "Laplace", control = list > (usePQL > FALSE) )) > > > (I have attached outputs below) > What I don't understand is; I thought in model2 I have defined > center effect > as a random effect. Then how come I got center effects and trt*center > interactions under the fixed effects list on the output? Probably I > didn't > understand how to set up these models in lmer. Could anyone help me > about > this? > > Thanks a lot for your help... > > Emine > > model1 <- try(lmer(cbind( yvect, nvect-yvect) ~ 1 + (1 | center), > family = binomial, niter = 25, method = "Laplace", control = list > (usePQL > FALSE) )) > >> summary(model1) > Generalized linear mixed model fit using Laplace > Formula: cbind(yvect, nvect - yvect) ~ 1 + (1 | center) > Family: binomial(logit link) > AIC BIC logLik deviance > 236.817 238.0951 -116.4085 232.817 > Random effects: > Groups Name Variance Std.Dev. > center (Intercept) 0.088127 0.29686 > number of obs: 14, groups: center, 7 > > Estimated scale (compare to 1) 0.2672612 > > Fixed effects: > Estimate Std. Error z value Pr(>|z|) > (Intercept) -0.32084 0.14709 -2.1812 0.02917 * > --- > Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > > > ################## > model2 <- try(lmer(cbind( yvect, nvect-yvect) ~ trt*center + ( 1 | > center) > , > family = binomial, niter = 25, method = "Laplace", control = list > (usePQL > FALSE) )) > >> summary(model2) > Generalized linear mixed model fit using Laplace > Formula: cbind(yvect, nvect - yvect) ~ trt * center + (1 | center) > Family: binomial(logit link) > AIC BIC logLik deviance > 30 39.58586 -1.547024e-07 3.094048e-07 > Random effects: > Groups Name Variance Std.Dev. > center (Intercept) 5e-10 2.2361e-05 > number of obs: 14, groups: center, 7 > > Estimated scale (compare to 1) 0.2672612 > > Fixed effects: > Estimate Std. Error z value Pr(>|z|) > (Intercept) -1.060869 0.065372 -16.2282 < 2e-16 *** > trt2 1.118029 0.086842 12.8743 < 2e-16 *** > center2 -0.325428 0.504256 -0.6454 0.51869 > center3 -0.325440 0.504258 -0.6454 0.51868 > center4 0.655407 0.413449 1.5852 0.11292 > center5 -0.325433 0.504256 -0.6454 0.51869 > center6 -0.325421 0.504255 -0.6454 0.51870 > center7 0.655422 0.413448 1.5853 0.11291 > trt2:center2 0.673737 0.651313 1.0344 0.30094 > trt2:center3 -0.137183 0.651314 -0.2106 0.83318 > trt2:center4 -0.307083 0.583845 -0.5260 0.59891 > trt2:center5 -0.137203 0.651314 -0.2107 0.83316 > trt2:center6 1.654555 0.712419 2.3224 0.02021 * > trt2:center7 0.673705 0.651311 1.0344 0.30096 > --- > > _________________________________________________________________ > > Office Live http://clk.atdmt.com/MRT/go/aub0540003042mrt/direct/01/ > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting- > guide.html > and provide commented, minimal, self-contained, reproducible code.Dr. Hank Stevens, Assistant Professor 338 Pearson Hall Botany Department Miami University Oxford, OH 45056 Office: (513) 529-4206 Lab: (513) 529-4262 FAX: (513) 529-4243 http://www.cas.muohio.edu/~stevenmh/ http://www.muohio.edu/ecology/ http://www.muohio.edu/botany/ "E Pluribus Unum"
On 6/1/07, emine ?zg?r Bayman <eozgur10 at hotmail.com> wrote:> Dear R users, > I'm pretty new on using lmer package. My response is binary and I have fixed > treatment effect (2 treatments) and random center effect (7 centers). I want > to test the effect of treatment by fitting 2 models: > > Model 1: center effect (random) only > Model 2: trt (fixed) + center (random) + trt*center interaction. > > Then, I want to compare these 2 models with Likelihood Ratio Test. Here are > my lmer codes that I don't feel comfortable about their correctness. > > model1 <- try(lmer(cbind( yvect, nvect-yvect) ~ 1 + (1 | center), > family = binomial, niter = 25, method = "Laplace", control = list(usePQL > FALSE) )) > > model2 <- try(lmer(cbind( yvect, nvect-yvect) ~ trt*center + ( 1 | center) > , > family = binomial, niter = 25, method = "Laplace", control = list(usePQL > FALSE) ))As you have seen, that model includes center as both a fixed effect and as a grouping factor for a random effect. You don't want to do that. I think the model you want is either cbind(yvect, nvect-yvect) ~ trt + (trt|center) which allows for correlation between the random effect for the intercept and the random effect for treatment within center, or cbind(yvect, nvect-yvect) ~ trt + (1|center/trt) which is equivalent to cbind(yvect, nvect-yvect) ~ trt + (1|center) + (1|center:trt) This model has a random effect for the center and a random effect for each center:trt combination but the variance of the latter random effects are assumed to be constant. Thus there are a total of three variance components (random effects for center, for center:trt combinations and for the per-observation noise term) in this model. The previous model will have 1 + (ntrt * (ntrt + 1))/2 variances and covariances to be estimated (where ntrt is the number of levels of the treatment factor). I hope this helps.> (I have attached outputs below) > What I don't understand is; I thought in model2 I have defined center effect > as a random effect. Then how come I got center effects and trt*center > interactions under the fixed effects list on the output? Probably I didn't > understand how to set up these models in lmer. Could anyone help me about > this? > > Thanks a lot for your help... > > Emine > > model1 <- try(lmer(cbind( yvect, nvect-yvect) ~ 1 + (1 | center), > family = binomial, niter = 25, method = "Laplace", control = list(usePQL > FALSE) )) > > >summary(model1) > Generalized linear mixed model fit using Laplace > Formula: cbind(yvect, nvect - yvect) ~ 1 + (1 | center) > Family: binomial(logit link) > AIC BIC logLik deviance > 236.817 238.0951 -116.4085 232.817 > Random effects: > Groups Name Variance Std.Dev. > center (Intercept) 0.088127 0.29686 > number of obs: 14, groups: center, 7 > > Estimated scale (compare to 1) 0.2672612 > > Fixed effects: > Estimate Std. Error z value Pr(>|z|) > (Intercept) -0.32084 0.14709 -2.1812 0.02917 * > --- > Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > > > ################## > model2 <- try(lmer(cbind( yvect, nvect-yvect) ~ trt*center + ( 1 | center) > , > family = binomial, niter = 25, method = "Laplace", control = list(usePQL > FALSE) )) > > >summary(model2) > Generalized linear mixed model fit using Laplace > Formula: cbind(yvect, nvect - yvect) ~ trt * center + (1 | center) > Family: binomial(logit link) > AIC BIC logLik deviance > 30 39.58586 -1.547024e-07 3.094048e-07 > Random effects: > Groups Name Variance Std.Dev. > center (Intercept) 5e-10 2.2361e-05 > number of obs: 14, groups: center, 7 > > Estimated scale (compare to 1) 0.2672612 > > Fixed effects: > Estimate Std. Error z value Pr(>|z|) > (Intercept) -1.060869 0.065372 -16.2282 < 2e-16 *** > trt2 1.118029 0.086842 12.8743 < 2e-16 *** > center2 -0.325428 0.504256 -0.6454 0.51869 > center3 -0.325440 0.504258 -0.6454 0.51868 > center4 0.655407 0.413449 1.5852 0.11292 > center5 -0.325433 0.504256 -0.6454 0.51869 > center6 -0.325421 0.504255 -0.6454 0.51870 > center7 0.655422 0.413448 1.5853 0.11291 > trt2:center2 0.673737 0.651313 1.0344 0.30094 > trt2:center3 -0.137183 0.651314 -0.2106 0.83318 > trt2:center4 -0.307083 0.583845 -0.5260 0.59891 > trt2:center5 -0.137203 0.651314 -0.2107 0.83316 > trt2:center6 1.654555 0.712419 2.3224 0.02021 * > trt2:center7 0.673705 0.651311 1.0344 0.30096 > --- > > _________________________________________________________________ > > Office Live http://clk.atdmt.com/MRT/go/aub0540003042mrt/direct/01/ > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. >