Đánh giá chức năng ikelihood trong các mô hình hỗn hợp tuyến tính (lme4)

Question

Tôi hiện đang viết kịch bản để đánh giá chức năng khả năng đăng nhập (bị hạn chế) để sử dụng trong các mô hình hỗn hợp tuyến tính. Tôi cần nó để tính toán khả năng của một mô hình với một số tham số cố định với các giá trị tùy ý. Có lẽ tập lệnh này cũng hữu ích đối với một số bạn!

Tôi sử dụng lmer() từ lme4 và logLik() để kiểm tra xem tập lệnh của tôi có hoạt động bình thường không. Và như nó có vẻ, nó không! Vì nền tảng giáo dục của tôi không thực sự quan tâm đến cấp độ thống kê này, tôi hơi lạc mất khi tìm ra sai lầm.

Tiếp theo, bạn sẽ tìm thấy một ví dụ kịch bản ngắn bằng cách sử dụng sleepstudy dữ liệu:

# * * * * * * * * * * * * * * * * * * * * * * * * 
    # * 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) 

Có thể có một số LMM-chuyên gia ở đây những người có thể giúp đỡ? Trong mọi trường hợp, các khuyến nghị của bạn được đánh giá cao!

chỉnh sửa: BTW, hàm likelihood cho LMMS có thể được tìm thấy trong Harville (1977), (hy vọng) truy cập thông qua liên kết này: Maximum likelihood approaches to variance component estimation and to related problems

Kính trọng, Simon

Answer 1

Giải pháp (như xây đá cẩm 2013) đã cài đặt phiên bản phát triển của lme4 và sử dụng đối số devFunOnly.

Phiên bản phát triển đó, cùng với khả năng này, có sẵn với lme4 on CRAN kể từ ngày 14 tháng 3 năm 2014 và reference guide cung cấp giải thích khen lời nhận xét của tác giả gói (Ben Bolker).