16
$\begingroup$

One of the problems I've always had with mixed models is figuring out data visualizations - of the kind that could end up on a paper or poster - once one has the results.

Right now, I'm working on a Poisson mixed effects model with a formula that looks something like the following:

a <- glmer(counts ~ X + Y + Time + (Y + Time | Site) + offset(log(people))

With something fitted in glm() one could easily use the predict() to get predictions for a new data set, and build something off of that. But with output like this - how would you construct something like a plot of the rate over time with the shifts from X (and likely with a set value of Y)? I think one could predict the fit well enough just from the Fixed effects estimates, but what about the 95% CI?

Is there anything else someone can think of that help visualize results? The results of the model are below:

Random effects:
 Groups     Name        Variance   Std.Dev.  Corr          
 Site       (Intercept) 5.3678e-01 0.7326513               
            time        2.4173e-05 0.0049167  0.250        
            Y           4.9378e-05 0.0070270 -0.911  0.172 

Fixed effects:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -8.1679391  0.1479849  -55.19  < 2e-16
X            0.4130639  0.1013899    4.07 4.62e-05
time         0.0009053  0.0012980    0.70    0.486    
Y            0.0187977  0.0023531    7.99 1.37e-15

Correlation of Fixed Effects:
         (Intr) Y    time
X      -0.178              
time    0.387 -0.305       
Y      -0.589  0.009  0.085
Improve this question
$\endgroup$
6
  • 1
    $\begingroup$ (+1) @EpiGrad: Why are you concerned about the CI (i.e. about the standard error) of the predictions from the fixed-effect part of your model? $\endgroup$
    – boscovich
    Commented Oct 4, 2012 at 10:37
  • 1
    $\begingroup$ @andrea An intellectual answer, and a practical one: Intellectually, I generally favor quantifying and visualizing uncertainty when I can. Practically, because I'm pretty sure a reviewer will ask for it. $\endgroup$
    – Fomite
    Commented Oct 4, 2012 at 10:43
  • $\begingroup$ Yeah yeah, sure, but I meant something different. My comment wasn't clear enough, sorry. You write in your question "but what about the 95% CI?". My comment is: why don't you calculate the standard error of the prediction from the fixed-effect part of the model? If you're able to calculate the predicted values form the fixed-effect part, then you're able to calculate the SE too and, thus, the CI. @EpiGrad $\endgroup$
    – boscovich
    Commented Oct 4, 2012 at 10:55
  • $\begingroup$ @andrea Ah. The concern is that one of the things I'd like to predict, time, also has a random effect, which I've no idea what to do with. $\endgroup$
    – Fomite
    Commented Oct 4, 2012 at 10:56
  • $\begingroup$ Well, you want to predict counts, not time. You fix values of X, Y and time and using the fixed-effects part of your model you predict counts. It's true that time is included in your model also as a random effect (just like the intercept and Y), but it doesn't matter here because using only the fixed-effect part of your model for the prediction is like setting the random effects to 0 @EpiGrad $\endgroup$
    – boscovich
    Commented Oct 4, 2012 at 11:23

1 Answer 1

Reset to default
4
$\begingroup$

Predicting counts using the fixed-effects part of your model means that you set to zero (i.e. their mean) the random effects. This means that you can "forget" about them and use standard machinery to calculate the predictions and the standard errors of the predictions (with which you can compute the confidence intervals).

This is an example using Stata, but I suppose it can be easily "translated" into R language:

webuse epilepsy, clear
xtmepoisson seizures treat visit || subject: visit
predict log_seiz, xb
gen pred_seiz = exp(log_seiz)
predict std_log_seiz, stdp
gen ub = exp(log_seiz+invnorm(.975)*std_log_seiz)
gen lb = exp(log_seiz-invnorm(.975)*std_log_seiz)

tw (line pred_seiz ub lb visit if treat == 0, sort lc(black black black) ///
 lp(l - -)), scheme(s1mono) legend(off) ytitle("Predicted Seizures") ///
 xtitle("Visit")

The graph refers to treat == 0 and it's intended to be an example (visit is not a really continuous variable, but it's just to get the idea). The dashed lines are 95% confidence intervals.

enter image description here

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.