2
$\begingroup$

Let's say you have a lmer model that test the drug effect of a set of rats with a set of rats (Control):

lme1 <- lmer(lVolume ~ Days*Drug + (Days|Drug))

where Days is the times of the drug applied, lVolume the log(Volume), and Drug a 2-factor level variable (Control and DrugA).

Dat0 <- list(
    list("Rat1", Volume=c(85,90 , 95, 120, 140, 170, 175, 185,190,240,250,300,320), Days = c(0,7,14,21,28,35,42,49,56,63,70,77,84), Drug = "Control"),
    list("Rat2", Volume=c(100,105,130,190,210,250,360,460,475,500,570,680,781), Days = c(0,7,14,21,28,35,42,49,56,63,70,77,84), Drug = "Control"),
    list("Rat3", Volume=c(120,125,150,155,280,281,350,360,390,430,440,550,670), Days = c(0,7,14,21,28,35,42,49,56,63,70,77,84), Drug = "Control"),
    list("Rat4", Volume=c(130,120,70,40,39,15,13,11,1,1,1,1,1), Days = c(0,7,14,21,28,35,42,49,56,63,70,77,84), Drug = "DrugA"),
    list("Rat5", Volume=c(140,110,85,81,80,79,85,65,60,90,110,105,110), Days = c(0,7,14,21,28,35,42,49,56,63,70,77,84), Drug = "DrugA"),
    list("Rat6", Volume=c(80,78,55,60,65,70,90,70,55,65,60,75,80), Days = c(0,7,14,21,28,35,42,49,56,63,70,77,84), Drug = "DrugA"))

DF <- purrr::map_dfr(dat0,
                       ~ data.frame(Rat = .[[1]],
                                    lVolume =.$log(Volume), Days = 
                                .$Days, Drug = .$Drug))

Then, you look at the estimated marginal means of such interaction model with

# average marginal effect of interaction
emmeans(lme1, c("Days", "Drug"))

Days      Drug   emmean    SE    df lower.CL upper.CL
 29.5 Control     6.044 0.176  6.04    5.615     6.47
 29.5   DrugA     0.617 0.484 49.46   -0.355     1.59

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

##Also
contrast(emmeans(lme1,  ~ Days*Drug), method="pairwise") 

contrast                                            estimate  SE   df   t.ratio p.value
 Control 29.5166666666667 - DrugA 29.5166666666667    5.43 0.515 41.7  10.544  <.0001

Degrees-of-freedom method: kenward-roger 


#p value
emmeans(lme1, specs = Drug ~ Days*Drug, adjust = "none")%>%
  summary(infer = TRUE)

Drug       Days emmean    SE    df lower.CL upper.CL t.ratio p.value
Control    29.5  6.044 0.176  6.04    5.615     6.47  34.421  <.0001
DrugA      29.5  0.617 0.484 49.46   -0.355     1.59   1.275  0.2082

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

My first question here is, how do you interpret these emmeans results in the first part and in the second part how to interpret the p-value estimated?

I plotted the emmeans, where the y axis is the lVolume, x axis the Days, red line is the estimated marginal mean for the Control, and the blue line is the estimated marginal mean for the DrugA.

My second question here is how can I compute the area under the curve between these two lines ? Is there a way to compute it with the emmeans() in R ? And also how would you interpret these two lines in terms of the drug effect over the Control part?

enter image description here

I tried to do:

probs <- predict(lme1, DF, type="response")
plot(roc(DF$lSize, probs), print.auc = TRUE)

But I got this extrange result:

Setting levels: control = 1.9505558410977, case = 2.74351132379233
Setting direction: controls < cases
Area under the curve: 1
Warning message:
In roc.default(response, predictor, auc = TRUE, ...) :
  'response' has more than two levels. Consider setting 'levels' explicitly or using 'multiclass.roc' instead

enter image description here

Thanks a lot for any hints/help/suggestions to understand this.

$\endgroup$
4
  • $\begingroup$ From the graph, the area you want is essentially that of a triangle, since the lines nearly coincide at the left. That area is 0.5*b*h by 3rd grade math, where h is the height which is equal to the span of Days, and b is the base which is proportional to the estimate of 5.43 you have at 29.5 days. Fr,i'm that, you can figure out what multiple of 5.43 is the area of that triangle; but you already have a test of it since the area is proportional to the difference that you already have tested. $\endgroup$
    – Russ Lenth
    Apr 3, 2022 at 0:30
  • $\begingroup$ @RussLenth I tried to compute auc(DF$lSize, probs) but what do I have AUC = 1 ? Setting levels: control = 1.9505558410977, case = 2.74351132379233 Setting direction: controls < cases Area under the curve: 1 Warning message: In roc.default(response, predictor, auc = TRUE, ...) : 'response' has more than two levels. Consider setting 'levels' explicitly or using 'multiclass.roc' instead Where probs <- predict(lme1, DF, type="response") $\endgroup$
    – Rosa Maria
    Apr 3, 2022 at 1:32
  • $\begingroup$ The thing is that I have several models, and I thought that auc() could be a good test to get this area by using the model. $\endgroup$
    – Rosa Maria
    Apr 3, 2022 at 1:34
  • $\begingroup$ Please note that your two lines do not comprise an ROC curve because you don't have a plot of sensitivity vs. specificity. That said, sometimes people do use areas below or between curves as an overall measure of an effect. $\endgroup$
    – Russ Lenth
    Apr 3, 2022 at 3:38

1 Answer 1

2
$\begingroup$

The bottom line of this answer is that the difference of the means is proportional to the area between the lines, and the constant of proportionality is just the range of the Days. So you might as well use the difference of the means and let it go at that.

If you do

EMM <- emmeans(model, "Drug", cov.reduce = range)

then a reference grid is created consisting of the four combinations of the two drugs at the lowest and highest values of Days, then the predictions are averaged over the two days. Then

pairs(EMM)

... takes the difference between those two means. Note that if the two lines in your picture have end values $y_{11}$ and $y_{12}$ for the first line and $y_{21}$ and $y_{22}$, then this result is equal to $\frac12(y_{11} + y_{12} - y_{21} - y_{22})$. Meanwhile, it is easy to show (see below) that the area between the two lines is $\frac b2(y_{11} + y_{12} - y_{21} - y_{22})$ where $b = t_{max} - t_{min}$, the range of the Days variable. Thus, the AUC is in fact the constant $b$ times the result of pairs(EMM). Thus, comparing those two means is essentially the same as computing the area between the lines.

To show that the area is as shown, first consider the case where the lines do not intersect; then we have a trapezoid of area $b$ times the average of the left difference and the right difference, which is the same as $b$ times the difference of the averages. Now, if the upper line is lowered by $\epsilon$, then the area is decreased by $b\epsilon$ and it keeps decreasing linearly even when the lines cross (where we then have the difference of areas of two triangles).

If you have curves instead of lines, then instead of cov.reduce = range, use at = list(Days = c(...)) and put in several equally-spaced values spanning the minimum to the maximum.

$\endgroup$
11
  • $\begingroup$ First of all thanks for the answer, it is very much appreciated. Following this, to make I did it correctly, pairs(EMM) = contrast estimate SE df t.ratio p.value Control - DrugA 6.96 0.645 48.8 10.793 <.0001 this is the estimates for pairs(EMM) is 6.96, then Days = 29.51667, therefore AUC = 29.5* 6.96 = 205.32 . $\endgroup$
    – Rosa Maria
    Apr 3, 2022 at 9:18
  • $\begingroup$ If you are going to compare several models using this AUC idea, then you had better use exactly the same Days range for all of them, e.g. at = list(Days = c(0, 30)) $\endgroup$
    – Russ Lenth
    Apr 3, 2022 at 12:55
  • 1
    $\begingroup$ Actually, you can include the scale factor in pairs: emm = emmeans(model, "Drug", at = list(Days = c(0,30)) then pairs(emm, scale = 30) $\endgroup$
    – Russ Lenth
    Apr 3, 2022 at 14:13
  • $\begingroup$ Yes, the thing is that not all the lines intersect, there are some that overlap, so I am interpreting this as more AUC more effect of the DrugA. It is the same model lme1 but different scenarios of the rats. So, I also was thinking is it possible to extract the coefficients that I get forn the intercept and the slope from each line, and then compute the finite integral, and just rest them to check the AUC under this (I need to be very careful with the intervals of the inteegral)? The only problem is that I don't know how to extract those individual intercepts and slopes from the lmer, do you? $\endgroup$
    – Rosa Maria
    Apr 3, 2022 at 14:27
  • $\begingroup$ Use fixef(model). But why belabor it, when you already have code that does the required computation? $\endgroup$
    – Russ Lenth
    Apr 3, 2022 at 14:38

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.