5

私は現在、線形混合モデルで使用するための（制限された）対数尤度関数を評価するためのスクリプトを書いています。いくつかのパラメーターが任意の値に固定されているモデルの尤度を計算するために必要です。たぶん、このスクリプトはあなたの何人かにも役立つでしょう！

lmer()fromlme4とを使用しlogLik()て、スクリプトが正常に機能するかどうかを確認します。そして、そう思われるように、そうではありません！私の学歴はこのレベルの統計にあまり関心がなかったので、間違いを見つけるのに少し迷いました。

次に、sleepstudy-dataを使用した短いサンプルスクリプトを示します。

  # * * * * * * * * * * * * * * * * * * * * * * * * 
  # * 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）にあり、（うまくいけば）このリンクからアクセスできます： 分散成分推定および関連する問題への最尤法

よろしく、サイモン

4

1 に答える 1

2

解決策（2013年3月現在）は、の開発バージョンをインストールし、引数lme4を利用することでした。devFunOnly

この開発バージョンは、この機能とともに、2014年3月14日からCRANのlme4で利用可能であり、リファレンスガイドには、パッケージ作成者（Ben Bolker）の元の質問に対するコメントを補足する説明が記載されています。

于 2014-05-01T16:21:00.483 に答える