Графики остатков модели смешанных эффектов с использованием функции ggplot2

Я пытаюсь построить график остаточных эффектов модели смешанных эффектов, используя функцию ggplot2. Однако после выполнения поиска я нашел некоторые доступные функции, но мне кажется, что для функции nlme они не работают.

Графики, которые я собираюсь построить, соответствуют приведенному ниже примеру:

Данные здесь.

Данные: https://drive.google.com/file/d/19mykz4B7jkTilbtwPQb3NUI09YZwohhs/view?usp=sharing

Вычислительные подпрограммы, которые я пробовал изначально, приведены ниже, см. ошибки, которые появляются при выполнении функции в ggplot2.

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

setwd("C:\\Users\\Desktop")
datanew1 = read.table("dadosnew.csv", header = T, sep=";", dec = ",")
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 

Другой способ, которым я пытался построить график, был с функцией ниже, но мне это не удалось.

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.


person user55546    schedule 25.01.2021    source источник
comment
Похоже, вы упустили какой-то код из шагов импорта/очистки данных, потому что я не могу понять, какие столбцы находятся в этом наборе данных (один столбец имеет заголовок с содержимым Variable2; Variable3; Response; Variable1 ;DummyVariable, а два других не помечены). Однако я подозреваю, что ваша проблема заключается в том, что nmle создает собственный объект S3 для хранения своей информации, и вам придется заглянуть в него, чтобы увидеть, где находятся данные для построения. ggplot2, кажется, изо всех сил пытается привести ваши данные к прямоугольному формату - проверьте, что такое str (completemodel).   -  person Dubukay    schedule 26.01.2021
comment
@Dubukay Там я обновил код необходимой информацией и недостающим пакетом сплайнов. Обратите внимание, что bs, weights, form, random — это функции пакетных сплайнов и nlme, а не переменные!   -  person user55546    schedule 26.01.2021
comment
Ах, европейский формат csv! Я должен был подумать, чтобы попробовать это. Похоже, мои подозрения оправдались — данные хранятся в объекте модели S3. Ваши имена столбцов (.fitted, .resid) предполагают, что данные должны быть сначала пропущены через fortify, но это не определено для объекта NMLE, поэтому возникает ошибка. Однако подгонка и невязки доступны с помощью completemodel$fitted и completemodel$residuals, если вы хотите рассчитать их самостоятельно.   -  person Dubukay    schedule 26.01.2021
comment
Не волнуйся. О, это слишком плохо!! Вы мне тогда говорите, что эти графики невозможно сделать в ggplot2, используя остатки, полученные настройкой функции nlme? Спасибо за информацию!!   -  person user55546    schedule 26.01.2021
comment
@Dubukay Я только что опубликовал решение.   -  person user55546    schedule 26.01.2021


Ответы (1)


[РЕШЕНИЕ]

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

datanew1 = read.table("E:/Downloads/dadosnew.csv", header = T, sep=";", dec = ",")

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")'

person user55546    schedule 24.02.2021