0
$\begingroup$

I have some rather basic questions that I suspect I know the answer to but would like to confirm.

I have a simple experiment measuring the number of a certain insect growth stage in the proportional presence of a dietary treatment. I have modeled the data using a binomial glmm. I would now like to plot a probability curve of the fixed effect (Treatment) using the sjPlot package. Following this tutorial, I am using "fe.slope" as the plot type. However, in the tutorial the y-axis is numbered from 0 - 1, as befits a probability curve, whereas my y-axis has a proportional scale from 0% - 100%. I assume it has something to do with the proportional aspect of my glmm. Can I simply change the y-axis scale manually to go from 0 - 1 without affecting the nature of the plot? Further, my Treatments are named 'T1, T2, T3, T4 and T5', yet the x-axis is labelled 0 - 4. Can I also manually change these without affecting the plot's true meaning? Finally, what do the points on this plot represent? Are they residuals?

Plot:

enter image description here

Code:

library(lme4)
library(sjPlot)

df<-as.data.frame(df)

df$OLRE <- seq_len(nrow(df)) #create observation level random effect

df[,2:3] %<>% lapply(function(x) as.integer(x)) #change to integers

m1<-glmer(Stage/Total~Treatment + (1|OLRE), data=df, family='binomial',weights=Total)

sjp.glmer(m1, type = "fe.slope",show.ci=TRUE)

Data:

dput(df)
structure(list(Treatment = c("0_diatoms", "0_diatoms", "0_diatoms", 
"25_diatoms", "25_diatoms", "25_diatoms", "50_diatoms", "50_diatoms", 
"50_diatoms", "75_diatoms", "75_diatoms", "75_diatoms", "100_diatoms", 
"100_diatoms", "100_diatoms"), Stage = c(45L, 50L, 50L, 58L, 
40L, 72L, 31L, 45L, 40L, 63L, 61L, 50L, 23L, 33L, 30L), Total = c(100, 
100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 
100), OLRE = 1:15), .Names = c("Treatment", "Stage", "Total", 
"OLRE"), row.names = c(NA, -15L), class = "data.frame")
$\endgroup$
  • $\begingroup$ Why do you divide Stage by Total when deriving your outcome variable and also including Total in the weights= option of your model? In other words, what is it that you are trying to model? What does Stage represent, what does Total represent and what does OLRE represent? $\endgroup$ – Isabella Ghement Mar 24 '18 at 19:50
  • $\begingroup$ Hi @IsabellaGhement, thanks for the comment. See this post, which has the same design as mine. OLRE refers to observation level random effect. $\endgroup$ – J.Con Mar 25 '18 at 23:27
3
$\begingroup$

The plot you produced doesn't make sense for your data, since Treatment is supposed to be a categorical variable (not an integer variable as per your plot).

I would re-fit the model by forcing Treatment to be a factor. I think OLRE should be a factor too - usually, the random effect goes hand in hand with a grouping variable (i.e., a factor in R's parlance).

So maybe try something like this:

df$Outcome <- df$Stage/df$Total

df$Treatment <- factor(df$Treatment)

df$OLRE <- factor(df$OLRE)


m1 <- glmer(Outcome ~ Treatment + (1|OLRE), data=df,
            family='binomial',weights=Total)


set.seed(17)

p1 <- sjp.glmer(m1, type = "pred.fe", vars = c("Treatment"), 
                show.ci = FALSE, show.scatter=FALSE)
p1

p2 <- sjp.glmer(m1, type = "pred", vars = c("Treatment"), 
                show.ci = FALSE, show.scatter=FALSE)
p2

According to the help file for sjp.glmer, which is available at https://www.rdocumentation.org/packages/sjPlot/versions/2.4.0/topics/sjp.glmer, the options type = "pred.fe" and "type = "pred" will help plot the predicted values against the response, with the predicted values being derived only from the fixed effects (for type = "pred.fe") or conditional on random intercept (for "type = "pred").

The beauty with saving the two plots is that you can then inspect the output produced by R when typing the plot names and "steal" that output to create your own customized graph. For example, the output for the plot p1 is as follows:

> p1
$data
x            y ci.low ci.high resp.y grp
1  1 0.4832236     NA      NA   0.45   1
2  1 0.4832236     NA      NA   0.50   1
3  1 0.4832236     NA      NA   0.50   1
4  3 0.5684214     NA      NA   0.58   1
5  3 0.5684214     NA      NA   0.40   1
6  3 0.5684214     NA      NA   0.72   1
7  4 0.3855458     NA      NA   0.31   1
8  4 0.3855458     NA      NA   0.45   1
9  4 0.3855458     NA      NA   0.40   1
10 5 0.5807438     NA      NA   0.63   1
11 5 0.5807438     NA      NA   0.61   1
12 5 0.5807438     NA      NA   0.50   1
13 2 0.2850919     NA      NA   0.23   1
14 2 0.2850919     NA      NA   0.33   1
15 2 0.2850919     NA      NA   0.30   1

$plot

attr(,"class")
[1] "sjPlot"        "sjpglm.ppresp"
Warning message:
In eval(family$initialize) : non-integer #successes in a binomial glm!

So the component p1$data of the plot gives you access to the data used to create the plot.

For example, for your first treatment,0_diatoms, you can compute the predicted probability by hand as follows and compare it against the one produced by sjp.glmer for type = "pred.fe":

intercept <- summary(m1)$coeff["(Intercept)","Estimate"]

exp(intercept)/(1 + exp(intercept))

The probability comes out to be 0.4832236 and this is listed under p1$data (in the column titled y) for each of the three observations corresponding to the first treatment group (i.e., 0_diatoms). The observations themselves are listed under the column titled resp.y of p1$data.

For the second treatment, the predicted probability computed by hand would be:

intercept <- summary(m1)$coeff["(Intercept)","Estimate"]

slope25_diatoms <- summary(m1)$coeff["Treatment25_diatoms","Estimate"]

exp(intercept + slope25_diatoms)/(1 + exp(intercept + slope25_diatoms))

hence 0.5684214. This is listed under the column y in p$data for each of the observations in your 25_diatoms treatment group.

The predicted probabilities for the first two treatment groups that are produced by sjp.glmer with the option type="pred" can be replicated by hand with the commands below:

intercept <- summary(m1)$coeff["(Intercept)","Estimate"]

ranef(m1)$OLRE[1, ]
ranef(m1)$OLRE[2, ]
ranef(m1)$OLRE[3, ]

exp(intercept + ranef(m1)$OLRE[1, ])/
  (1 + exp(intercept + ranef(m1)$OLRE[1, ]))

exp(intercept + ranef(m1)$OLRE[2, ])/
  (1 + exp(intercept + ranef(m1)$OLRE[2, ]))

exp(intercept + ranef(m1)$OLRE[3, ])/
  (1 + exp(intercept + ranef(m1)$OLRE[3, ]))

and

intercept <- summary(m1)$coeff["(Intercept)","Estimate"]

slope25_diatoms <- summary(m1)$coeff["Treatment25_diatoms","Estimate"]

exp(intercept + slope25_diatoms + ranef(m1)$OLRE[4, ])/
 (1 + exp(intercept + slope25_diatoms + ranef(m1)$OLRE[4, ]))

exp(intercept + slope25_diatoms + ranef(m1)$OLRE[5, ])/
 (1 + exp(intercept + slope25_diatoms + ranef(m1)$OLRE[5, ]))

exp(intercept + slope25_diatoms + ranef(m1)$OLRE[6, ])/
 (1 + exp(intercept + slope25_diatoms + ranef(m1)$OLRE[6, ]))

These probabilities come out to be 0.4646612, 0.492606 and 0.492606 for the first treatment (i.e., 0_diatoms) and 0.5748393, 0.4744047 and 0.6512543 for the second treatment (i.e., 25_diatoms). You will recognize then in the first 6 rows of column y contained in the R output corresponding to p2$data:

> p2$data
x            y ci.low ci.high resp.y grp
1  1 0.4646612     NA      NA   0.45   1
2  1 0.4926060     NA      NA   0.50   1
3  1 0.4926060     NA      NA   0.50   1
4  3 0.5748393     NA      NA   0.58   1
5  3 0.4744047     NA      NA   0.40   1
6  3 0.6512543     NA      NA   0.72   1
7  4 0.3447117     NA      NA   0.31   1
8  4 0.4209869     NA      NA   0.45   1
9  4 0.3934527     NA      NA   0.40   1
10 5 0.6078542     NA      NA   0.63   1
11 5 0.5968777     NA      NA   0.61   1
12 5 0.5358618     NA      NA   0.50   1
13 2 0.2574995     NA      NA   0.23   1
14 2 0.3081921     NA      NA   0.33   1
15 2 0.2927040     NA      NA   0.30   1

The above should provide enough insight into what sjp.glmer is actually plotting to enable you to generate your own plots. In particular, you can plot the predicted probabilities on the 0 to 1 scale (rather then converting them to percentage points, as sjp.glmer does). You can also force sjp.glmer to return uncertainty intervals around the predicted probability by using the show.ci = TRUE option.

I am not sure why sjp.glmer jitters the response values in the column resp.y before plotting them when you include the option show.scatter=TRUE. I would have thought that using show.scatter=FALSE would have turned off the random jittering, but it actually removes the response observations altogether from the plot. So either use set.seed(17) (or some other random seed) before creating the plot and then invoke the option show.scatter=TRUE for sjp.glmer, or create your own plot from p1$data or p2$data (whichever you need), but superimpose the values stored in y.resp on your plot, which are the actual response values. If you don't set the random seed, sjp.glmer will come up with a random configuration of jittered response values every time you invoke it for the same model.

Disclaimer: I do not know enough about your study, data and model to determine if your current modeling approach is sensible. Assuming it is, the code I provided here will make sense. Otherwise, it might not make sense.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you so much for taking the time to write this very comprehensive answer. $\endgroup$ – J.Con Mar 26 '18 at 22:42
  • 1
    $\begingroup$ You're welcome! I hope it gave you enough background info to enable you to get the plot you want. $\endgroup$ – Isabella Ghement Mar 26 '18 at 22:43
  • $\begingroup$ Yes, definitely. I can now confirm if my plots are correct using the manual procedures you provided. I have now 'stolen' the info and made a graph with ggplot()! $\endgroup$ – J.Con Mar 26 '18 at 22:45
  • 1
    $\begingroup$ That's great! There's nothing like understanding what it is that you are plotting! 😀 $\endgroup$ – Isabella Ghement Mar 26 '18 at 22:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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