В настоящее время я работаю с линейными смешанными моделями (LMM) и хочу провести некоторые остаточные анализы. Поскольку остатки от линейных смешанных моделей «коррелированы и не обязательно имеют постоянную дисперсию», они не способны проверить предположения, лежащие в основе модели (Fitzmaurice et al. 2011, p. 266). Таким образом, было предложено (Waternaux et al. 1989) использовать преобразованные по Холески остатки вместо того, чтобы исследовать допущения модели, см. Также Houseman et al. (2004) и Santos Nobre & da Motta Singer (2007).
То есть предположим, что у нас есть LMM вида Y i = X i β + Z i b i + 101 i , где β - вектор параметров фиксированных эффектов, b i - вектор нормально распределенных случайных эффектов, и где X i и Z i - матрицы дизайна.
Затем невязки e i = Y i - X i β коррелируют, как и их наблюдаемый аналог, использующий оцененные β-коэффициенты вместо неизвестные истинные значения. Мы обозначим эти остатки как r i . Декорреляция остатков осуществляется путем нахождения разложения Холецкого ковариационной матрицы r i . То есть, пусть Cov (e i ) = Σ i , пусть L i i '- разложение Холецкого оценки Σ i на основе r i . Преобразованные по Холецкому остатки затем определяются как r * i = L i -1 r i . Это остатки, которые я пытаюсь получить из пакета lme4 (Bates et al. 2015), но пока безуспешно.
Однако, используя данные с домашней страницы Fitzmaurice et al. (2011), мне удалось воссоздать остатки из тематического исследования в книге, используя пакет nlme (в отличие от пакета lme4). Код указан ниже.
require(foreign) # To read dta-files
require(nlme) # The nlme package used to estimate the LMM.
require(mgcv) # Needed for the extract.lme.cov function
fat <- read.dta("https://content.sph.harvard.edu/fitzmaur/ala2e/fat.dta") # We read the data
# Add spline function for time, with a knot at age of menarch, denoted time = 0 in the data.
fat$timeknot <- fat$time * (fat$time > 0)
# Estimate the LMM with nlme
nlme.model <- lme(pbf ~ time + timeknot, data=fat, random = ~ 1 + time + timeknot | id)
# Calculate the residuals
raw.residuals <- residuals(nlme.model, level=0)
# Cholesky residuals
est.cov <- extract.lme.cov(nlme.model, fat) # We extract the blocked covariance function of the residuals
Li <- t(chol(est.cov)) # We find the Cholesky transformation of the residuals. (The transform is to get the lower trangular matrix.)
cholesky.residuals <- solve(Li) %*% raw.residuals # We then calculate the transformed residuals.
hist(cholesky.residuals) # We plot the residuals
Мой вопрос, таким образом, есть ли простой способ извлечь оценочную ковариационную матрицу остатков, используя пакет lme4? Или получить преобразованные по Холецкому остатки напрямую?
Благодарен за любую помощь или указатели. С уважением,
/ Фил
Ссылки:
Бейтс Д., Мехлер М., Болкер Б. и Уокер С. (2015). Подгонка линейных моделей со смешанными эффектами с использованием lme4. Журнал статистического программного обеспечения , 67 (1), 1-48.
Fitzmaurice, G.M., Laird, N.M. & Ware, J.H. (2011). Прикладной продольный анализ , второе издание. Джон Вили и сыновья.
Houseman, E.A., Ryan, L.M. & Coll, B.A. (2004). Остатки Холецкого для оценки нормальных ошибок в линейной модели с коррелированными результатами. Журнал Американской статистической ассоциации n, 99 (466), 383-394.
Сантос Нобре, Дж., И да Мотта Сингер, Дж. (2007). Остаточный анализ для линейных смешанных моделей. Биометрический журнал: Журнал математических методов в бионауках , 49 (6), 863-875.
Waternaux, C., Laird, N.M. & Ware, J.H. (1989). Методы анализа продольных данных: концентрация свинца в крови и когнитивное развитие. Журнал Американской статистической ассоциации , 84 (405), 33-41.