4

I have a transplant experiment for four locations and four substrates (taken from each location). I have determined survival for each population in each location and substrate combination. This experiment was replicated three times.

I have created a lmm as follows:

Survival.model <- lmer(Survival ~ Location + Substrate + Location:Substrate + (1|Replicate), data=Transplant.Survival,, REML = TRUE)

I would like to use the predict command to extract predictions, for example:

Survival.pred <- predict(Survival.model)

Then extract standard errors so that I can plot them with the predictions to generate something like the following plot:

enter image description here

I know how to do this with a standard glm (which is how I created the example plot), but am not sure if I can or should do this with an lmm.

Can I do this or am I as a new user of linear mixed models missing something fundamental?

I did find this post on Stack Overflow which was not helpful.

Based on a comment from RHertel, maybe I should have phrased the question: How do I plot model estimates and confidence intervals for my lmer model results so that I can get a similar plot to the one I have created above?

Sample Data:

Transplant.Survival <- structure(list(Location = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L), .Label = c("Steninge", "Molle", 
"Kampinge", "Kaseberga"), class = "factor"), Substrate = structure(c(1L, 
1L, 1L, 2L, 2L, 2L, 3L, 3L, 4L, 4L, 4L, 1L, 1L, 2L, 2L, 2L, 3L, 
3L, 3L, 4L, 4L, 4L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 
4L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L), .Label = c("Steninge", 
"Molle", "Kampinge", "Kaseberga"), class = "factor"), Replicate = structure(c(1L, 
2L, 3L, 1L, 2L, 3L, 1L, 2L, 1L, 2L, 3L, 2L, 3L, 1L, 2L, 3L, 1L, 
2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L), .Label = c("1", 
"2", "3"), class = "factor"), Survival = c(1, 1, 1, 0.633333333333333, 
0.966666666666667, 0.5, 0.3, 0.233333333333333, 0.433333333333333, 
0.966666666666667, 0.866666666666667, 0.5, 0.6, 0.266666666666667, 
0.733333333333333, 0.6, 0.3, 0.5, 0.3, 0.633333333333333, 0.9, 
0.266666666666667, 0.633333333333333, 0.7, 0.633333333333333, 
0.833333333333333, 0.9, 0.6, 0.166666666666667, 0.333333333333333, 
0.433333333333333, 0.6, 0.9, 0.6, 0.133333333333333, 0.566666666666667, 
0.633333333333333, 0.633333333333333, 0.766666666666667, 0.766666666666667, 
0.0333333333333333, 0.733333333333333, 0.3, 1.03333333333333, 
0.6, 1)), .Names = c("Location", "Substrate", "Replicate", "Survival"
), class = "data.frame", row.names = c(NA, -46L))
Community
  • 1
  • 1
Keith W. Larson
  • 1,543
  • 2
  • 19
  • 34
  • I'm not familiar with the `lme4` package but I think that you need to supply test data to `predict()`, and not just the model. The trained model is an important first step, but without input data I don't see how it could be used to predict the outcome. You could try using `predict(Survival.model, test_data)` where `test_data` is the data set for which you want to obtain a prediction. – RHertel Aug 25 '15 at 06:21
  • RHertel: I am sorry for the stupid question, but for a small study such as this, does it make sense to split the data into a training and test dataset? Would it be more appropriate then to plot my estimates with confidence intervals? – Keith W. Larson Aug 25 '15 at 06:23
  • 1
    Read the [r-sig-mixed-models FAQ](http://glmm.wikidot.com/faq). It contains an example. Also, there is [this answer](http://stackoverflow.com/a/14435982/1412059) of mine that points to a thread on the mailing lists that discusses the finer conceptional problems of this. – Roland Aug 25 '15 at 06:50
  • @Roland: thank you for the response. I have read through the pages and I understand the approach, but with such a small dataset, it seems splitting it into a training and test dataset seems problematic. I will have to rethink how I what I want to present my model results in a plot or the modelling approach all together! – Keith W. Larson Aug 25 '15 at 08:55
  • I don't see where the FAQ says that you need to split into training and test data to get confidence intervals / standard errors for predictions. – Roland Aug 25 '15 at 09:35
  • The issue is that SE's are not standard for MLMs, which is why `lme4` hasn't implemented the functions for them. Do you want the FE standard error only? The RE only? The combined standard error? Also, the structures of your variance matrices matter. – alexwhitworth Aug 25 '15 at 20:11
  • But you can use `VarCorr(Survival.model)` along with the Z and X model matrices to get the correct terms – alexwhitworth Aug 25 '15 at 20:17
  • @Roland: I understand that there are not SEs for MLMs, so I gave a simplified example of what I would like to obtain, i.e. a plot with predicted survival with some kind of SE or CI. As I want to visually present my results from my model and this seems an intuitive, if not easy approach, way to present them. I have read the FAQ, and see the example for 'lmer', but what I do not understand being new to these methods, from the example is where the 'newdat' comes from. – Keith W. Larson Aug 26 '15 at 05:04
  • @alex: can you give me an example that I can learn from to use VarCorr(Survival.model) along with the Z and X model matrices to get the correct terms? – Keith W. Larson Aug 26 '15 at 05:04
  • `newdat` is the data you want to predict for. It can be identical to the data you fitted to. – Roland Aug 26 '15 at 07:18
  • @KeithLarson There is an example in the [MLM FAQ](http://glmm.wikidot.com/faq) which was already linked above. Any of the standard texts (Goldstein, Pinheiro and Bates, etc) will also have the details. – alexwhitworth Aug 26 '15 at 16:07

1 Answers1

7

Edit: fixed bug in function / figure.

If you like to plot estimates with CI, you may want to look at the sjp.lmer function in the sjPlot package. See some example of the various plot types here.

Furthermore, the arm package provides function for computing standard Errors (arm::se.fixef and arm::se.ranef)

sjp.setTheme("forestgrey") # plot theme
sjp.lmer(Survival.model, type = "fe")

would give following plot

enter image description here

Daniel
  • 7,252
  • 6
  • 26
  • 38