Graphs of the mixed effects model residuals using the ggplot2 function

Question

I am trying to graph the residual effects of the mixed effects model using the ggplot2 function. However, after performing a search I found some functions available but what seems to me is that for the function nlme they are not working.

The graphs that I intend to make are those of the example below:

The computational routines that I tried initially are below, see the errors that are appearing when executing the function in ggplot2.

library(splines)
library(ggplot2)
library(nlme)
library(gridExtra)

datanew1$DummyVariable = as.factor(datanew1$DummyVariable)
datanew1$Variable2 = as.factor(datanew1$Variable2)
datanew1$Variable3 = as.factor(datanew1$Variable3)
#############################################################################
############################## Model ########################################
#############################################################################
model <-  lme(Response~(bs(Variable1, df=3)) + DummyVariable,
                         random=~1|Variable2/Variable3, datanew1, method="REML")
completemodel <- update(model, weights = varIdent(form=~1|DummyVariable))

p1 <- qplot(.fitted, .resid, data = completemodel) +
  geom_hline(yintercept = 0) +
  geom_smooth(se = FALSE)

Erro: `data` must be a data frame, or other object coercible by `fortify()`, not an S3 object with class lme
Run `rlang::last_error()` to see where the error occurred.

p2 <- qplot(sample =.stdresid, data = completemodel, stat = "qq") + geom_abline()
grid.arrange(p1,p2)

Erro: `data` must be a data frame, or other object coercible by `fortify()`, not an S3 object with class lme
Run `rlang::last_error()` to see where the error occurred.
Além disso: Warning message:
`stat` is deprecated

Another way in which I tried to make the graph was with the function below, but I was not successful.

ggplot(completemodel, aes(.fitted, .resid)) + geom_point()

Erro: `data` must be a data frame, or other object coercible by `fortify()`, not an S3 object with class lme
Run `rlang::last_error()` to see where the error occurred.

It seems like you've left some code out of the data import/cleanup steps, because I can't understand what the columns are in this data set (a single column has a header, with the contents "Variable2;Variable3;Response;Variable1;DummyVariable", and the other two are unlabeled). However, I suspect your problem is that `nmle` is creating a custom S3 object to hold its info, and you'll have to peer into that to see where the data to plot is. ggplot2 seems to be struggling to coerce your data into a rectangular format - check what str(completemodel) is. — Dubukay, Jan 25 '21 at 22:22
@Dubukay There, I updated the code with the necessary information and the missing splines package. See that "bs", "weights", "form", "random" are functions of packages splines and nlme and not variables! — user55546, Jan 25 '21 at 23:02
Ah, the European csv format! I should've thought to try that. It looks like my suspicions were correct - the data is stored within the S3 model object. Your column names (.fitted, .resid) suggest that the data should be run through `fortify` first, but that's not defined for the NMLE object, thus the error. The fitted and residuals are available with completemodel$fitted and completemodel$residuals if you'd like to calculate them yourself, however. — Dubukay, Jan 25 '21 at 23:34
Don't worry. Oh, that’s too bad!! Are you telling me then that these graphs are not possible to be made in ggplot2 using the residues obtained by adjusting the nlme function? Thanks for the information!! — user55546, Jan 26 '21 at 00:03

user55546 · Accepted Answer · 2021-10-15T14:25:17.690

[SOLUTION]

library(splines)
library(ggplot2)
library(nlme)
library(gridExtra)

datanew1$DummyVariable = as.factor(datanew1$DummyVariable)
datanew1$Variable2 = as.factor(datanew1$Variable2)
datanew1$Variable3 = as.factor(datanew1$Variable3)


model <-  lme(Response~(bs(Variable1, df=3)) + DummyVariable,
              random=~1|Variable2/Variable3, datanew1, method="REML")
completemodel <- update(model, weights = varIdent(form=~1|DummyVariable))

df_model <- broom.mixed::augment(completemodel)
#> Registered S3 method overwritten by 'broom.mixed':
#>   method      from 
#>   tidy.gamlss broom
df_model[".stdresid"] <- resid(completemodel, type = "pearson")

p1 <- ggplot(df_model, aes(.fitted, .resid)) + 
  geom_point() +
  geom_hline(yintercept = 0) +
  geom_smooth(se=FALSE)

p2 <- ggplot(df_model, aes(sample = .stdresid)) +
  geom_qq() +
  geom_qq_line()

grid.arrange(p1,p2)
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

Graphs of the mixed effects model residuals using the ggplot2 function

1 Answers1