線形混合モデル（lme4）のikelihood関数の評価

Question

現在、線形混合モデルで使用する（制限された）対数尤度関数を評価するためのスクリプトを作成しています。私はいくつかのパラメータを任意の値に固定してモデルの可能性を計算する必要があります。 多分、このスクリプトは皆さんにとっても役に立ちます！

を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-専門家があるのですか？いずれにしても、あなたのお勧めは大いに評価されます！

編集：ところで、LMMSの尤度関数はHarville（1977）で見つけることができる、（うまくいけば）このリンクを介してアクセス：月のよう Maximum likelihood approaches to variance component estimation and to related problems

よろしく、 サイモン

Answer 1

溶液（ 2013）は、lme4の開発バージョンをインストールし、devFunOnly引数を使用することでした。

この開発版では、lme4 on CRANと1414年3月14日から利用できます。reference guideは、パッケージ作成者（Ben Bolker）が元の質問に対してコメントしたことを説明しています。