Оценка функции оценки в MLE для модели lme в nlme и R - PullRequest
Новая
       11

Оценка функции оценки в MLE для модели lme в nlme и R

0 голосов
/ 27 октября 2018

В качестве проверки работоспособности для более поздних целей я пытаюсь понять, как связаны компоненты объекта lme.Я пытаюсь оценить функции оценки в MLE, полученных из объекта lme, созданного с использованием gamm, чтобы убедиться, что они равны нулю.Я представил минимальный рабочий пример.Мои функции оценки были получены из уравнения (7) в https://arxiv.org/pdf/1612.04911.pdf.

library(lme4)
library(mgcv)
library(lasso2) # for Prostate data set
data(Prostate)

model = gamm(lpsa ~lcavol + s(lweight), data=Prostate)$lme

Здесь я вычисляю объекты, необходимые для оценки функций оценки 

X = model.matrix(formula(model), data=model$data)
YminXbeta = c(model$data$y - X %*% fixed.effects(model))
varcomps = VarCorr(model)
varcomps = as.numeric(varcomps[nrow(varcomps) - (1:0),1])
Sigmainv = solve(extract.lme.cov2(model, data=model$data)$V)
Sigmainvsq = Sigmainv %*% Sigmainv

ZVZT = (model$data$Xr %*% (getVarCov(model, 
 type='random.effects')/varcomps[1]) %*% t(model$data$Xr))
SigmainvZVZT = Sigmainv %*% ZVZT
SigmainvZVZTsq = SigmainvZVZT %*% SigmainvZVZT

Теперь я оцениваю оценки наMLE 

# the scores
# the score for the mean parameter
sbeta = t(X) %*% Sigmainv %*% YminXbeta
#                 [,1]
# X(Intercept)    5.329071e-14
# Xlcavol         1.190159e-13
# Xs(lweight)Fx1 -1.110223e-15

# the score for the error variance parameter
ssigmasq = (-sum(diag(Sigmainv)) + YminXbeta %*% Sigmainvsq %*% 
YminXbeta)/2
# 3.664974e-08

# the score for the random effect parameter
stausq = (-sum(diag(SigmainvZVZT)) + YminXbeta %*% SigmainvZVZT 
%*% Sigmainv %*% YminXbeta)/2
# -7.507903

Обратите внимание, что компонент дисперсии случайных эффектов близок к нулю 

varcomps[1] # 2.665509e-09

, но я не уверен, что это сделает производную ненулевой.

Почемуфункция оценки для члена компонента отклонения не близка к нулю?Я делаю ошибку или неправильно понимаю, что за объекты в объекте lme?

...