Şu anda doğrusal karışık modellerde kullanım için (kısıtlı) log-olabilirity fonksiyonunu değerlendirmek üzere bir betik yazıyorum. Bir modelin rasgele değerlere sabitlenmiş bazı parametreler ile olasılığını hesaplamak için buna ihtiyacım var. Belki de bu betik sizin için de faydalıdır!Lineer karışık modellerde (lme4) ikelilik fonksiyonunun değerlendirilmesi
Komut dosyanızın gerektiği gibi çalışıp çalışmadığını kontrol etmek için lme4 ve logLik()'dan lmer() ve logLik() kullanıyorum. Ve göründüğü gibi değil! Eğitim geçmişim, bu istatistik düzeyiyle gerçekten ilgili değildi, ben de bu hataları bulmakta biraz kaydım.
# * * * * * * * * * * * * * * * * * * * * * * * *
# * 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)
ardından, sleepstudy-verileri kullanarak kısa örnek senaryosunu bulacak? Her durumda, önerileriniz büyük ölçüde takdir edilmektedir!
düzenleme: Maximum likelihood approaches to variance component estimation and to related problems
Saygılar, Simon
I ** güçlü bir imleme fonksiyonunu döndürme yeteneğine sahip 'lme4' (muhtemelen github, 'devtools' aracılığıyla), geliştirme sürümünü almanızı öneririm –
Teşekkürler! Lme4'ün github versiyonuna baktım ve bunu “devtools” kullanarak kurdum. 'DevFunOnly = T' argümanı ve ürettiği işlev hakkında daha fazla belge var mı? Özellikle ortaya çıkan sapma fonksiyonuna beslemek zorunda olduğum argümanlar ile ilgileniyorum, çünkü bu yine benim için en önemli adım! \ Kod {devFunOnly} \ kod {DOĞRU} \ kodu {teta} vektörü temsil eden tek bir sayısal vektör argüman alır olduğunda – SimonG
sapma fonksiyonu döndü. Bu vektör, Cholesky parametrelendirmesindeki rastgele efektlerinin varyans-kovaryans fonksiyonunu tanımlar. Tek bir rastgele etkisi için, bu rastgele efekti standart sapmasına eşit bir değerdir ... –