【问题标题】：Unable to fit model using “lrm.fit”无法使用“lrm.fit”拟合模型
【发布时间】：2019-01-13 08:53:38
【问题描述】：

我已经对各种结果进行了一系列逻辑回归，并且一直使用rms 包计算 Nagelkerke R^2，直到现在没有任何问题......

当我尝试在具有肥胖 (=1) 或“正常”体重 (=0) 的二元结果的回归中使用它时，我收到一条消息，指出我的公式中有错误，无法拟合模型使用lrm.fit。我尝试将我的结果编码为数字和一个因素，但得到了同样的错误。我还尝试使用数据的一个子集，不包括既不是肥胖也不是“正常”体重的病例。那也行不通。

我在下面粘贴了我的公式。（我无法包含数据，因为我无权共享它。）我还粘贴了使用 glm 函数获得的回归输出。

我做错了什么？我确信已经发布了解决方案，但在 Stack Overflow 或其他论坛上找不到解决方案。（我的主管要求我在结果中包含伪 R^2。我还包括平均边际效应。）

提前感谢您的帮助！

代码

require(rms)
logit.obnorm.base <- lrm(obnorm13 ~ wgt9 + gender + birthweight + breastfed + 
pcgweight + smoked + pcgeducation + income + welfare, data = gui)
print(logit.obnorm.base)

错误信息

Error in lrm(obnorm13 ~ wgt9 + gender + birthweight + breastfed + pcgweight +  : 
  Unable to fit model using “lrm.fit”

使用`glm` 进行回归时的回归输出

glm(formula = obnorm13 ~ wgt9 + gender + birthweight + breastfed + 
    pcgweight + smoked + pcgeducation + income + welfare, family = binomial(link = "logit"), 
    data = gui, weights = bdwg01, na.action = na.omit)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.4402  -0.1509  -0.0876  -0.0618   6.2779  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)   -5.73687    0.73321  -7.824 5.10e-15 ***
wgt92          3.25940    0.20237  16.106  < 2e-16 ***
wgt93          6.38153    0.28246  22.593  < 2e-16 ***
gender         0.26701    0.17444   1.531 0.125854    
birthweight    0.14304    0.15234   0.939 0.347751    
breastfed      0.24493    0.19173   1.277 0.201430    
pcgweight2     0.35537    0.22073   1.610 0.107401    
pcgweight3     1.11191    0.21671   5.131 2.89e-07 ***
smoked         0.50183    0.19956   2.515 0.011913 *  
pcgeducation2 -0.71452    0.21294  -3.356 0.000792 ***
pcgeducation3 -0.81295    0.27853  -2.919 0.003515 ** 
pcgeducation4 -0.88217    0.32623  -2.704 0.006848 ** 
income        -0.01348    0.07976  -0.169 0.865744    
welfare        0.11014    0.06298   1.749 0.080290 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2410.1  on 4758  degrees of freedom
Residual deviance: 1037.4  on 4745  degrees of freedom
  (2536 observations deleted due to missingness)
AIC: 1048.4

Number of Fisher Scoring iterations: 7

`dput(head(mydata))` 的输出

> dput(head(gui))
structure(list(zid01 = c(1000, 2000, 3000, 4000, 5000, 6000), 
    bdwg01 = c(0.222555072843037, 1.95115048656779, 1.08399948558217, 
    0.697544449190721, 0.910021907765603, 0.278315557267728), 
    gender = c(0, 0, 0, 0, 0, 0), wgt9 = structure(c(1L, 1L, 
    1L, 1L, 1L, 1L), .Label = c("1", "2", "3"), class = "factor"), 
    apcb01 = c(3.4, 3.1, 3.4, 3.3, 3.6, 2.8), smoked = c(1, 1, 
    0, 0, 1, 1), breastfed = c(0, 0, 0, 0, 1, 0), pcgweight = structure(c(2L, 
    1L, 1L, 1L, 2L, 2L), .Label = c("1", "2", "3"), class = "factor"), 
    pcgeducation = structure(c(4L, 1L, 4L, 4L, 4L, 2L), .Label = c("1", 
    "2", "3", "4"), class = "factor"), welfare = c(2, 7, 2, 3, 
    2, 3), adsd59b = c(5, 2, 4, 3, 5, 2), aded08a = c(51, 31, 
    35, 44, 52, 30), aded08b = c(11, 10, 4, 10, 11, 8), aded08c = c(15, 
    7, 8, 11, 14, 8), aded08d = c(10, 7, 7, 8, 8, 7), aded08e = c(12, 
    6, 9, 11, 13, 5), aded08f = c(11, 5, 9, 8, 10, 5), aded08g = c(8, 
    4, 7, 8, 10, 6), aded10a = c(0, 8, 4, 2, 1, 2), aded10b = c(1, 
    4, 1, 1, 0, 0), aded10c = c(4, 10, 4, 3, 2, 4), aded10d = c(1, 
    3, 1, 0, 2, 0), aded10e = c(9, 7, 9, 10, 5, 9), aded10f = c(6, 
    25, 10, 6, 5, 6), aded11a = c(1.8, 3, 1.2, 1.2, 2.8, 2), 
    aded11b = c(1.4, 5, 1.2, 1.4, 1.2, 1.8), aded11c = c(4.8, 
    5, 4.8, 5, 5, 4.6), aded11d = c(3.8, 3.4, 4, 4.2, 4.2, 3.6
    ), obese13 = c(0, 0, 0, 0, 0, 0), obnorm13 = c(0, 0, 0, 0, 
    0, 0), xswgt13 = c(0, 0, 0, 0, 0, 0), up = c("0", "0", "0", 
    "0", "0", "0"), down = c(NA_character_, NA_character_, NA_character_, 
    NA_character_, NA_character_, NA_character_), obese9 = c(0, 
    0, 0, 0, 0, 0), obnorm9 = c("0", "0", "0", "0", "0", "0"), 
    xswgt9 = c(0, 0, 0, 0, 0, 0), wgt9.1 = c("1", "1", "1", "1", 
    "1", "1"), income = c(5, 2, 4, 3, 5, 2), birthweight = c(3.4, 
    3.1, 3.4, 3.3, 3.6, 2.8), PHtotal = c(51, 31, 35, 44, 52, 
    30), PHbehaviour = c(11, 10, 4, 10, 11, 8), PHintel_school = c(15, 
    7, 8, 11, 14, 8), PHphys_app = c(10, 7, 7, 8, 8, 7), PHfreedom_anx = c(12, 
    6, 9, 11, 13, 5), PHpopularity = c(11, 5, 9, 8, 10, 5), PHhappiness = c(8, 
    4, 7, 8, 10, 6), SDQEmotional = c(0, 8, 4, 2, 1, 2), SDQConduct = c(1, 
    4, 1, 1, 0, 0), SDQHyperactivity = c(4, 10, 4, 3, 2, 4), 
    SDQPeer_Rel = c(1, 3, 1, 0, 2, 0), SDQProsocial = c(9, 7, 
    9, 10, 5, 9), SDQTotalDiff = c(6, 25, 10, 6, 5, 6), EASShyness = c(1.8, 
    3, 1.2, 1.2, 2.8, 2), EASEmotionality = c(1.4, 5, 1.2, 1.4, 
    1.2, 1.8), EASActivity = c(4.8, 5, 4.8, 5, 5, 4.6), EASSociability = c(3.8, 
    3.4, 4, 4.2, 4.2, 3.6), factor.wgt9 = structure(c(1L, 1L, 
    1L, 1L, 1L, 1L), .Label = c("1", "2", "3"), class = "factor"), 
    factor.pcged = structure(c(4L, 1L, 4L, 4L, 4L, 2L), .Label = c("1", 
    "2", "3", "4"), class = "factor"), factor.pcgweight = structure(c(2L, 
    1L, 1L, 1L, 2L, 2L), .Label = c("1", "2", "3"), class = "factor"), 
    orig.wgt9 = c("1", "1", "1", "1", "1", "1"), orig.pcgweight = c(2, 
    1, 1, 1, 2, 2), orig.pcgeducation = c(4, 1, 4, 4, 4, 2), 
    PHintelschool = c(15, 7, 8, 11, 14, 8)), .Names = c("zid01", 
"bdwg01", "gender", "wgt9", "apcb01", "smoked", "breastfed", 
"pcgweight", "pcgeducation", "welfare", "adsd59b", "aded08a", 
"aded08b", "aded08c", "aded08d", "aded08e", "aded08f", "aded08g", 
"aded10a", "aded10b", "aded10c", "aded10d", "aded10e", "aded10f", 
"aded11a", "aded11b", "aded11c", "aded11d", "obese13", "obnorm13", 
"xswgt13", "up", "down", "obese9", "obnorm9", "xswgt9", "wgt9.1", 
"income", "birthweight", "PHtotal", "PHbehaviour", "PHintel_school", 
"PHphys_app", "PHfreedom_anx", "PHpopularity", "PHhappiness", 
"SDQEmotional", "SDQConduct", "SDQHyperactivity", "SDQPeer_Rel", 
"SDQProsocial", "SDQTotalDiff", "EASShyness", "EASEmotionality", 
"EASActivity", "EASSociability", "factor.wgt9", "factor.pcged", 
"factor.pcgweight", "orig.wgt9", "orig.pcgweight", "orig.pcgeducation", 
"PHintelschool"), row.names = c(NA, 6L), class = "data.frame")

【问题讨论】：

'gui' 是我的数据集，它在其他任何地方都没有出现任何问题。我需要做什么来修复它？（我对 R 还是很陌生——我是一名硕士生。）
如果您只需要那个伪 R 方块，您可以使用库 DescTools 中的 PseudoR2。
谢谢@tmfmnk。我会试一试的。
谢谢@jaySF。我已经编辑了这个问题以包含它。很遗憾。我无法包含实际数据集的样本，因为我无权共享它。
@tmfmnk - 非常感谢您的帮助！这完全没有任何问题。还要感谢 jaySF。我是 R 新手，非常感谢像你这样回复这些 cmets 并帮助像我这样的新手的人。

标签： r logistic-regression model-fitting rms

【解决方案1】：

解决方案

tmfmnk 建议的计算 R^2 的替代方法没有任何问题。

我把它贴在这里以防它帮助别人。谢谢大家的回复。

# Install the DescTools package and open the library   
install.packages("DescTools")
    library(DescTools)

# Run the baseline logistic regression for obese / no excess weight at 13, including weights 
summary(logit.obnorm.base <- glm(obnorm13 ~ wgt9 + gender + birthweight + breastfed + 
pcgweight + smoked + pcgeducation + income + welfare, data = gui, weights = bdwg01,
family = binomial(link = "logit"), na.action=na.omit))

# Calculate Nagelkerke R^2
PseudoR2(logit.obnorm.base, which = "Nagelkerke")
>Nagelkerke 
 0.7309776

【讨论】：

【解决方案2】：

我遇到了类似的问题，我找到了解决方案here。问题可能是 MLE 没有在多个步骤中收敛。我在参数中添加了maxit=1000，它起作用了。

【讨论】：

代码

错误信息

使用glm 进行回归时的回归输出

dput(head(mydata)) 的输出

解决方案

使用`glm` 进行回归时的回归输出

`dput(head(mydata))` 的输出