3

我被要求为我在 R 中运行的回归提供伪 R2 值。我正在使用 vglm 命令对 4 类患者生活质量数运行回归,给定患者肺动脉应变、脉压和最大肺动脉压。

但是,vglm 只提供 Residual Deviation,不提供 Null Deviance。我从以前的帖子中查找的伪 R2 拟合优度测量都需要零偏差。有没有办法从 vglm 获取或计算零偏差? https://www.stat.auckland.ac.nz/~yee/VGAM/doc/VGAMrefcard.pdf

谢谢,肖娜

> summary(fit)

Call:
vglm(formula = resp_var ~ Strain + PP + Pmax, family = multinomial, data = DF)

Pearson Residuals:
                   Min       1Q    Median        3Q    Max
log(mu[,1]/mu[,4]) -3.8849 -0.45457  0.237302  0.456967 4.9936
log(mu[,2]/mu[,4]) -2.7173 -0.45528 -0.255479  0.601062 2.1737
log(mu[,3]/mu[,4]) -2.7175 -0.15174 -0.096334 -0.041271 6.6040

Coefficients:
               Estimate Std. Error    z value
(Intercept):1 -18.37926   1609.958 -0.0114160
(Intercept):2 -22.19514   1609.958 -0.0137862
(Intercept):3 -24.72028   1609.958 -0.0153546
Strain:1      286.66750   3989.303  0.0718590
Strain:2      285.95551   3989.303  0.0716806
Strain:3      284.37742   3989.303  0.0712850
PP:1           -0.87220     67.951 -0.0128356
PP:2           -0.77133     67.951 -0.0113512
PP:3           -0.93263     67.951 -0.0137250
Pmax:1          0.15845     28.841  0.0054941
Pmax:2          0.17548     28.841  0.0060843
Pmax:3          0.28930     28.841  0.0100307

Number of linear predictors:  3 

Names of linear predictors: 
log(mu[,1]/mu[,4]), log(mu[,2]/mu[,4]), log(mu[,3]/mu[,4])

Dispersion Parameter for multinomial family:   1

Residual deviance: 102.8416 on 222 degrees of freedom

Log-likelihood: -51.42081 on 222 degrees of freedom

Number of iterations: 19 
4

1 回答 1

2

零偏差是零(仅截取)模型的残余偏差,例如resp_var ~ 1。所以要计算 McFadden 的伪 R 平方:

# generate sample data
set.seed(101)
resp_var = t(rmultinom(100, 1, rep(0.25, 4)))
DF = data.frame(Strain = rnorm(100), PP = rnorm(100), Pmax = rnorm(100))

# fit models
m1 <- vglm(formula = resp_var ~ Strain + PP + Pmax, family = multinomial, data = DF)
null <- vglm(formula = resp_var ~ 1, family = multinomial, data = DF)

# McFadden's pseudo-R^2
print(pseudo_R2 <- 1 - deviance(m1) / deviance(null))

# [1] 0.05300655

例如,这同意mlogit

library(mlogit)
DF$cat_resp = factor(apply(resp_var, 1, function(.) which(. == 1) ) )
mDf <- mlogit.data(DF, choice = "cat_resp", shape = "wide")
summary(mlogit(cat_resp ~ 1|Strain + PP + Pmax, data = mDf))

# Call:
# mlogit(formula = cat_resp ~ 1 | Strain + PP + Pmax, data = mDf, 
#    method = "nr", print.level = 0)
# ...
# ...
# Log-Likelihood: -130.03
# McFadden R^2:  0.053007 
# Likelihood ratio test : chisq = 14.557 (p.value = 0.10386)
于 2013-09-26T23:00:47.427 回答