Obecnie piszę skrypt do oceny (ograniczonej) funkcji log-wiarygodności do wykorzystania w liniowych modelach mieszanych. Potrzebuję go, aby obliczyć prawdopodobieństwo modelu z niektórymi parametrami ustalonymi na arbitralne wartości. Może ten skrypt jest również pomocny dla niektórych z was!Ocena funkcji ikelihood w liniowych modelach mieszanych (lme4)
Używam lmer() od lme4 i logLik(), aby sprawdzić, czy mój skrypt działa tak, jak powinien. A jak się wydaje, tak nie jest! Ponieważ moje wykształcenie nie zajmowało się tym poziomem statystyk, jestem trochę zagubiony znajdując błąd.
Po znajdziesz krótki przykładowy skrypt używając sleepstudy-dane:
# * * * * * * * * * * * * * * * * * * * * * * * *
# * example data
library(lme4)
data(sleepstudy)
dat <- sleepstudy[ (sleepstudy$Days %in% 0:4) & (sleepstudy$Subject %in% 331:333) ,]
colnames(dat) <- c("y", "x", "group")
mod0 <- lmer(y ~ 1 + x + (1 | group), data = dat)
# + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + #
# #
# Evaluating the likelihood-function for a LMM #
# specified as: Y = X*beta + Z*b + e #
# #
# + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
# * * * * * * * * * * * * * * * * * * * * * * * *
# * the model parameters
# n is total number of individuals
# m is total number of groups, indexed by i
# p is number of fixed effects
# q is number of random effects
q <- nrow(VarCorr(mod0)$group) # number of random effects
n <- nrow(dat) # number of individuals
m <- length(unique(dat$group)) # number of goups
Y <- dat$y # response vector
X <- cbind(rep(1,n), dat$x) # model matrix of fixed effects (n x p)
beta <- as.numeric(fixef(mod0)) # fixed effects vector (p x 1)
Z.sparse <- t([email protected]) # model matrix of random effect (sparse format)
Z <- as.matrix(Z.sparse) # model matrix Z (n x q*m)
b <- as.matrix(ranef(mod0)$group) # random effects vector (q*m x 1)
D <- diag(VarCorr(mod0)$group[1:q,1:q], q*m) # covariance matrix of random effects
R <- diag(1,nrow(dat))*summary(mod0)@sigma^2 # covariance matrix of residuals
V <- Z %*% D %*% t(Z) + R # (total) covariance matrix of Y
# check: values in Y can be perfectly matched using lmer's information
Y.test <- X %*% beta + Z %*% b + resid(mod0)
cbind(Y, Y.test)
# * * * * * * * * * * * * * * * * * * * * * * * *
# * the likelihood function
# profile and restricted log-likelihood (Harville, 1997)
loglik.p <- - (0.5) * ( (log(det(V))) + t((Y - X %*% beta)) %*% solve(V) %*% (Y - X %*% beta) )
loglik.r <- loglik.p - (0.5) * log(det(t(X) %*% solve(V) %*% X))
#check: value of above function does not match the generic (restricted) log-likelihood of the mer-class object
loglik.lmer <- logLik(mod0, REML=TRUE)
cbind(loglik.p, loglik.r, loglik.lmer)
Być może istnieją pewne LMM-eksperci tutaj kto może pomóc? W każdym razie Twoje rekomendacje są bardzo doceniane!
EDIT: BTW, funkcja prawdopodobieństwa dla LMMS można znaleźć w Harville (1977), (mam nadzieję), dostępny za pośrednictwem tego linku: Maximum likelihood approaches to variance component estimation and to related problems
Pozdrowienia, Simon
Ja ** zdecydowanie ** zalecamy, aby uzyskać wersję rozwojową 'lme4' (prawdopodobnie z github, poprzez' devtools'), która ma zdolność ('mkDevfunOnly = TRUE') zwracania funkcji odchylenia –
Dzięki! Zajrzałem do wersji 'lme4' github i zainstalowałem ją za pomocą' devtools'. Czy jest jakaś dodatkowa dokumentacja dotycząca argumentu 'devFunOnly = T' i funkcji, którą tworzy? Szczególnie interesują mnie argumenty, które muszę przekazać wynikowej funkcji odchylenia, ponieważ to znowu jest dla mnie najważniejszy krok! – SimonG
funkcja odchylenia zwrócona, gdy \ code {devFunOnly} jest \ code {TRUE} przyjmuje pojedynczy numeryczny argument wektorowy reprezentujący wektor \ code {theta} . Ten wektor definiuje funkcję wariancji-kowariancji efektów losowych w parametryzacji Cholesky'ego. Dla pojedynczego losowego efektu jest to pojedyncza wartość równa odchyleniu standardowemu efektu losowego ... –