0

I have a mixed model with three factors: one categorical and two continuous. When trying to plot the fitted lines over the original data, it passes much far below. Why would that happen? Here is the model output:

Linear mixed-effects model fit by maximum likelihood
  Data: x 
   AIC      BIC    logLik
143.5392 155.3614 -64.76961

Random effects:
Formula: ~1 | col
     (Intercept) Residual
StdDev: 2.737332e-05 1.221771

Fixed effects:  ctmax ~ ur * hl + ratectmax 
            Value Std.Error DF   t-value p-value
(Intercept) 23.855561  2.868299 31  8.316971  0.0000
urhigh       8.717054  4.542828 31  1.918861  0.0642
hl           0.125498  0.035472 31  3.537982  0.0013
ratectmax    7.283523  2.284297 31  3.188518  0.0033
urhigh:hl   -0.080326  0.049294 31 -1.629521  0.1133

Here is the script to plot:

 plot(x$ctmax~x$hl,pch=24,bg=c("light grey","purple","transparent","red","yellow")[x$col], 
 ylim=c(20,48), family = "Times", cex.lab = ".85", cex.axis = ".85", xlab="heating rate   (°C/min)", ylab="Temperature (°C)", family = "Times")
abline(23.85, 0.12,col="black",lwd=1.8)

Here is the plot enter image description here

This method also plots nothing, despite it worked perfectly for a model with two continous factors only.

lines(predict(m7), col="black", lwd=1, lty=1)
Agus camacho
  • 868
  • 2
  • 9
  • 24
  • You are trying to plot 4 dimensional data on a 2 dimensional plot. With the type of model you've chosen. The other variables you aren't drawing affect the line but you are ignoring their effects. It's not clear to me what you are trying to accomplish because mathematically this doesn't quite make sense. If you just want marginal effects, them maybe this can help: https://cran.r-project.org/web/packages/ggeffects/vignettes/introduction_randomeffects.html – MrFlick Apr 12 '21 at 18:01
  • I was expecting to see that the lines representing the predictive model pass over the observed values. In that way, i see if the model's output, which accounts for all dimensions, actually predicts the observations made. However, nor the predicted coefficients or "predict" brought them back. ggpredict in the link you suggested brought the following error: ggpredict(m7,terms="hl") Error: Must extract column with a single valid subscript. x The subscript `terms[2]` can't be `NA`. Run `rlang::last_error()` to see where the error occurred. – Agus camacho Apr 13 '21 at 07:54

1 Answers1

0

In the end, this behaved as I expected: plotting the model's predicted values for each level of the original predictive factor, gave me much more reasonable lines.

x$pred=predict(m7)
points(x$ctmax[!x$ur%in%"50"]~x$hl[!x$ur%in%"50"],pch=24,family = "Times", cex.lab = ".85", cex.axis = ".85", xlab="Hydration level (%)", ylab=" Temperature (°C)", ylim=c(25,50),bg="light grey")
with(x,abline(lm(x$pred[!x$ur%in%"50"]~hl[!x$ur%in%"50"]), col="light grey"))
# dry treatment
points(x$ctmax[!x$ur%in%"85"]~x$hl[!x$ur%in%"85"],pch=24,family = "Times", cex.lab   = ".85", cex.axis = ".85", xlab="Hydration level (%)", ylab=" Temperature (°C)", ylim=c(25,50),bg="black")
with(x,abline(lm(x$pred[!x$ur%in%"85"]~hl[!x$ur%in%"85"]), col="black"))

enter image description here

Agus camacho
  • 868
  • 2
  • 9
  • 24