[R-sig-ME] FW: GLMM Question

[Ben Bolker]

  These are all good questions.  I'm taking the liberty of forwarding
it to r-sig-mixed-models, as that's a public venue where multiple
people can ask questions, answers are publicly viewable, and answers
are archived.

  Short answers: drop1() (my preference) or the anova method from the
'car' package are good for testing whole factors (you may also want to
see the ?pvalues help page in recent versions of the lme4 package).

  For translation to probabilities, you may want
predict(...,type="response").  You may also find the effects package
useful.

 http://glmm.wikidot.com/faq can also be helpful; if you have
particular suggestions for material to go there, you can add it
yourself or suggest it on the r-sig-mixed-models list and someone
(maybe me) might add it ...

  good luck
    Ben Bolker


On 14-03-14 11:18 AM, Brad White wrote:
> 
> 
> Dr. Bolker,
> 
> Hi, I imagine you receive quite a few questions regarding modeling 
> applications, but I have searched online to no avail.  I found
> your paper in Trends and Ecology and Evolution from 2010 very
> useful, but it didn’t answer all my questions.
> 
> 
> 
> I’m a researcher in cattle disease and we have worked in the area
> of evaluating disease risk among groups (count data) as well as
> monitoring behavioral changes (again, count data).  In the past, we
> have done events/trials modeling with GLMM to account for random
> effects of repeated measures on individual animal as well as
> housing (calves housed within pens and sometimes pens within
> operations). For a little background (probably more than you want
> or need), I attached one of our articles describing behavior and
> the figures represent the type of information I’d like to glean
> from these models.
> 
> 
> 
> I have previously used SAS GLIMMIX to model our outcomes, but
> would really like to model these with R as I think this would be
> helpful for our graduate students (and me) in the future.  Our
> projects are set up to test specific hypotheses, thus generating
> estimates and inferences are important to me. I’ve tried several
> methods in R, but can’t seem to figure out the best way to evaluate
> model output (ideally I’d like to have a test for the entire
> factor, not just the tests on individual parameters of the factor).
> Also, I typically transform the estimates from these models to
> probabilities, but am having trouble from this part also.  I would
> much appreciate some direction on places to find further resources
> on this, it appears from my online searching that this is not the
> common way R output has been used.
> 
> 
> 
> As a specific example, below are models created in glmer and
> glmmpql where success is the event happening (count of seconds) and
> failure is the total number of seconds – success.  Var1 is our
> treatment of interest (3 levels) and study day is our time factor
> of interest (again three levels).  When running these models, I get
> the estimate for Var1Hi, Var1Lo, but not for Var1CO as it is in the
> intercept. I would like to have a test for the overall effect of
> Var1 to test our hypothesis that Var1 impacts outcome: if so, then
> I would test individual treatments. Same is true for study day.
> Then I’d also like to create a table with the estimates (and SE)
> transformed to probabilities for the multiple variables.
> 
> 
> 
> Any help you can provide is much appreciated as I’ve struggled with
> this a while.
> 
> 
> 
> Thanks,
> 
> 
> 
> Brad
> 
> 
> 
> Glmer model:
> 
> 
> 
> model <- glmer(cbind(success, failure) ~ var1 + studyday + (1 |
> pen) + (1 | calfid),
> 
> family = binomial(link = "logit"), data = mydata)
> 
> summary(model)
> 
> 
> 
> Generalized linear mixed model fit by maximum likelihood
> ['glmerMod']
> 
> Family: binomial ( logit )
> 
> Formula: cbind(success, failure) ~ var1 + studyday + (1 | pen) +
> (1 |      calfid)
> 
> Data: mydata
> 
> 
> 
> AIC       BIC    logLik  deviance
> 
> 394776.0  394794.6 -197381.0  394762.0
> 
> 
> 
> Random effects:
> 
> Groups Name        Variance  Std.Dev.
> 
> calfid (Intercept) 4.352e+00 2.086e+00
> 
> pen    (Intercept) 3.237e-10 1.799e-05
> 
> Number of obs: 105, groups: calfid, 35; pen, 3
> 
> 
> 
> Fixed effects:
> 
> Estimate Std. Error z value Pr(>|z|)
> 
> (Intercept) -1.315422   0.604552   -2.18   0.0296 *
> 
> var1Hi       1.395121   0.853311    1.63   0.1021
> 
> var1Lo       1.256089   0.872415    1.44   0.1499
> 
> studyday1   -0.618371   0.001991 -310.62   <2e-16 ***
> 
> studyday2   -0.484639   0.001981 -244.64   <2e-16 ***
> 
> ---
> 
> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 
> 
> 
> Correlation of Fixed Effects:
> 
> (Intr) var1Hi var1Lo stdyd1
> 
> var1Hi    -0.708
> 
> var1Lo    -0.693  0.491
> 
> studyday1 -0.002  0.000  0.000
> 
> studyday2 -0.002  0.000  0.000  0.589
> 
> 
> 
> glmmPQL model
> 
> model <- glmmPQL(cbind(success, failure) ~ var1 + studyday, 
> random=list(~1|pen, ~1|calfid),
> 
> family = binomial(link = "logit"), data = mydata)
> 
> 
> 
> summary of glmmPQL model
> 
> Linear mixed-effects model fit by maximum likelihood
> 
> Data: mydata
> 
> AIC BIC logLik
> 
> NA  NA     NA
> 
> 
> 
> Random effects:
> 
> Formula: ~1 | pen
> 
> (Intercept)
> 
> StdDev:   0.2823529
> 
> 
> 
> Formula: ~1 | calfid %in% pen
> 
> (Intercept) Residual
> 
> StdDev: 0.0008658502 103.4659
> 
> 
> 
> Variance function:
> 
> Structure: fixed weights
> 
> Formula: ~invwt
> 
> Fixed effects: cbind(success, failure) ~ var1 + studyday
> 
> Value Std.Error DF    t-value p-value
> 
> (Intercept) -0.1944745 0.2510724 68 -0.7745754  0.4413
> 
> var1Hi       0.2526331 0.1879374 30  1.3442412  0.1889
> 
> var1Lo       0.1101132 0.1942646 30  0.5668208  0.5751
> 
> studyday1   -0.5603492 0.1997436 68 -2.8053418  0.0065
> 
> studyday2   -0.4560453 0.1985154 68 -2.2972798  0.0247
> 
> Correlation:
> 
> (Intr) var1Hi var1Lo stdyd1
> 
> var1Hi    -0.381
> 
> var1Lo    -0.373  0.497
> 
> studyday1 -0.460 -0.006  0.001
> 
> studyday2 -0.463 -0.007  0.003  0.584
> 
> 
> 
> Standardized Within-Group Residuals:
> 
> Min           Q1          Med           Q3          Max
> 
> -2.553151595 -0.664950798  0.007380211  0.706653488  2.626906196
> 
> 
> 
> Number of Observations: 105
> 
> Number of Groups:
> 
> pen calfid %in% pen
> 
> 3              35
> 
> 
> 
> 
> 
> 
> 
> Brad White, DVM, MS
> 
> Beef Production Medicine
> 
> Q211 Mosier Hall
> 
> Manhattan, KS 66506
> 
> Phone: 785-532-4243
> 
> 
>