선형 혼합 모델에서 ikelihood 함수를 계산합니다 (lme4)

Question

현재 선형 혼합 모델에서 사용하기 위해 (제한된) 대수 우도 함수를 평가하는 스크립트를 작성하고 있습니다. 임의의 값으로 고정 된 매개 변수를 사용하여 모델의 가능성을 계산해야합니다. 아마이 스크립트는 여러분 중 일부에게 도움이 될 것입니다!

나는 을 lme4과 logLik()에서 사용하여 내 스크립트가 제대로 작동하는지 확인합니다. 그리고 보이는 것처럼, 그렇지 않습니다! 제 교육 배경이이 정도의 통계에별로 관련이 없으므로 실수를 찾는 것을 조금 잃어 버렸습니다.

# * * * * * * * * * * * * * * * * * * * * * * * * 
    # * 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(mod0@Zt)       # 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) 

어쩌면 사람을 도울 수 있습니다 여기에 몇 가지 LMM-전문가가 :

에 따라, 당신은 sleepstudy 데이터를 사용하여 간단한 예제 스크립트를 찾을 수 있나요? 어쨌든 당신의 권고는 크게 평가됩니다!

편집 : Maximum likelihood approaches to variance component estimation and to related problems

감사합니다, 사이먼

Answer 1

솔루션 (월의 같은 : BTW, LMMS에 대한 우도 함수는 Harville (1977), (희망)이 링크를 통해 액세스 할 수에서 찾을 수 있습니다 2013)는 lme4의 개발 버전을 설치하고 devFunOnly 인수를 사용했습니다.

개발 버전이이 기능과 함께, lme4 on CRAN 이후 14 월 2014에 사용할 수 있음을

하고 reference guide은 원래의 질문에 패키지 저자의 (벤 Bolker) 주석을 칭찬 설명을 제공합니다.