2016-09-21 11 views
2

GMMAT(関数:glmmkin())というRパッケージを使用して、私のモデルを遺伝関係行列で調整した混合モデルロジスティック回帰を実行しました。ロジスティック回帰から手でlogLikを計算

モデルから私の出力は、(マニュアルから取られた)を含む:

  • theta:分散パラメータ推定値を[1]および分散成分のパラメータ推定値[2]
  • coefficients:固定効果パラメータ見積もり(傍受を含む)。
  • linear.predictors:線形予測子。
  • fitted.values:元の縮尺に適合した平均値。
  • Y:最終作業ベクトルのサンプルサイズに等しい長さのベクトル。
  • P:サンプルサイズに等しい寸法の投影行列。
  • residuals:元の縮尺での残差。分散パラメータによって再スケーリングされません。
  • cov:固定効果(迎撃を含む)のための共分散行列。
  • converged:収束のための論理指標。

分散を計算するために対数尤度を求めようとしています。私の最初の本能は、これを "手で"計算するためにlogLik.glm関数を引き離すことでした。私はAICを計算しようとしていました。私は答えをhereから使いました。

ロジスティック回帰をstats::glm()で実行して健全性チェックを行いましたが、model1$aicは4013.232ですが、スタックオーバーフローの回答を使用して30613.03を取得しました。

私の質問は - 私は上記のRでリストアップした出力を使ってロジスティック回帰からの対数尤度を手で計算する方法を知っていますか?

答えて

1

ここでは統計的な洞察はありません。ただ、私が見ているのはglm.fitです。この...

model <- glm(vs ~ mpg, mtcars, family = binomial(logit)) 
logLik(model) 
# 'log Lik.' -12.76667 (df=2) 

sparse_model <- model[c("theta", "coefficients", "linear.predictors", "fitted.values", "y", "P", "residuals", "cov", "converged")] 
get_logLik(sparse_model) 
#[1] -12.76667 
+0

感謝をモデルにフィットしながら、あなたが重みを指定していない(または、あなたがした場合、あなたはモデルオブジェクトでこれらの重みを含める必要があります)例えば

get_logLik <- function(s_model, family = binomial(logit)) { n <- length(s_model$y) wt <- rep(1, n) # or s_model$prior_weights if field exists deviance <- sum(family$dev.resids(s_model$y, s_model$fitted.values, wt)) mod_rank <- sum(!is.na(s_model$coefficients)) # or s_model$rank if field exists aic <- family$aic(s_model$y, rep(1, n), s_model$fitted.values, wt, deviance) + 2 * mod_rank log_lik <- mod_rank - aic/2 return(log_lik) } 

場合のみの作品はあなた、クリス!それは魅力のように働いた:) – mkv8

関連する問題