从逻辑回归中手动计算 logLik
Calculating logLik by hand from a logistic regression
我 运行 一个混合模型逻辑回归调整我的模型与遗传关系矩阵使用称为 GMMAT
的 R 包(功能:glmmkin()
)。
我的模型输出包括(取自用户手册):
theta
:离散参数估计[1]和方差分量参数估计[2]
coefficients
:固定效应参数估计值(包括截距)。
linear.predictors
:线性预测变量。
fitted.values
:原始尺度上的拟合平均值。
Y
:长度等于最终工作向量样本大小的向量。
P
:维度等于样本大小的投影矩阵。
residuals
:原始尺度上的残差。未通过色散参数重新缩放。
cov
:固定效应的协方差矩阵(包括截距)。
converged
:收敛的逻辑指标。
我正在尝试获取对数似然以计算解释的方差。我的第一直觉是拆开 logLik.glm
函数来计算这个 "by hand",但我一直在尝试计算 AIC。我使用了 here.
的答案
我用逻辑回归 运行 和 stats::glm()
进行了合理性检查,其中 model1$aic
是 4013.232 但使用我找到的 Stack Overflow 答案,我得到了 30613.03。
我的问题是 -- 有谁知道如何使用我在上面列出的 R 中的输出手动计算逻辑回归的对数似然?
这里没有统计见解,只是我从 glm.fit
中看到的解决方案。这仅在您在拟合模型时未指定权重时有效(或者如果您指定了权重,则需要在模型对象中包含这些权重)
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)
}
例如...
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
我 运行 一个混合模型逻辑回归调整我的模型与遗传关系矩阵使用称为 GMMAT
的 R 包(功能:glmmkin()
)。
我的模型输出包括(取自用户手册):
theta
:离散参数估计[1]和方差分量参数估计[2]coefficients
:固定效应参数估计值(包括截距)。linear.predictors
:线性预测变量。fitted.values
:原始尺度上的拟合平均值。Y
:长度等于最终工作向量样本大小的向量。P
:维度等于样本大小的投影矩阵。residuals
:原始尺度上的残差。未通过色散参数重新缩放。cov
:固定效应的协方差矩阵(包括截距)。converged
:收敛的逻辑指标。
我正在尝试获取对数似然以计算解释的方差。我的第一直觉是拆开 logLik.glm
函数来计算这个 "by hand",但我一直在尝试计算 AIC。我使用了 here.
我用逻辑回归 运行 和 stats::glm()
进行了合理性检查,其中 model1$aic
是 4013.232 但使用我找到的 Stack Overflow 答案,我得到了 30613.03。
我的问题是 -- 有谁知道如何使用我在上面列出的 R 中的输出手动计算逻辑回归的对数似然?
这里没有统计见解,只是我从 glm.fit
中看到的解决方案。这仅在您在拟合模型时未指定权重时有效(或者如果您指定了权重,则需要在模型对象中包含这些权重)
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)
}
例如...
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