Оценка функции правдоподобия в линейных смешанных моделях (lme4)

В настоящее время я пишу сценарий для оценки (ограниченной) функции логарифмического правдоподобия для использования в линейных смешанных моделях. Мне это нужно для расчета вероятности модели с некоторыми параметрами, фиксированными на произвольные значения. Возможно, этот скрипт будет полезен и для некоторых из вас!

Я использую lmer() из lme4 и logLik(), чтобы проверить, работает ли мой скрипт должным образом. И как кажется, это не так! Поскольку мое образование не очень-то касалось такого уровня статистики, я немного растерялся, обнаружив ошибку.

Далее вы найдете краткий пример скрипта, использующего данные sleepstudy:

  # * * * * * * * * * * * * * * * * * * * * * * * * 
  # * 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-специалисты, которые могут помочь? В любом случае ваши рекомендации очень ценны!

редактировать: Кстати, функцию правдоподобия для LMM можно найти в Harville (1977), (надеюсь) доступную по этой ссылке: Максимальная вероятность подходы к оценке компонента дисперсии и к связанным проблемам

С уважением, Саймон