时间戳

首页 > 解决方案 > R：与 R 中的 glm 和 mle2 包不同的结果

r - R：与 R 中的 glm 和 mle2 包不同的结果

问题描述

所以我想使用 GLM 找到估计参数并将其与 mle2 包进行比较。这是我的 GLM 代码

d <- read.delim("http://dnett.github.io/S510/Disease.txt")

d$disease=factor(d$disease)
d$ses=factor(d$ses)
d$sector=factor(d$sector)
str(d)
glm2 <- glm(disease~ses+sector, family=binomial(link=logit), data=d)
summary(glm2)

我的 mle2() 代码

y<-as.numeric(as.character(d$disease))
x1<-as.numeric(as.character(d$age))
x2<-as.numeric(as.character(d$sector))
x3<-as.numeric(as.character(d$ses))

library(bbmle)
nlldbin=function(A,B,C,D){
  eta<-A+B*(x3==2)+C*(x3==3)+D*(x2==2)
  p<-1/(1+exp(-eta))
  joint.pdf= (p^y)*((1-p)^(1-y))
  -sum(joint.pdf, log=TRUE ,na.rm=TRUE)
}
st <- list(A=0.0001,B=0.0001,C=0.0001,D=0.0001)
est_mle2<-mle2(start=st,nlldbin,hessian=TRUE)
summary(est_mle2)

但结果却截然不同。请帮我解决这个问题，谢谢！

> summary(est_mle2)
Maximum likelihood estimation

Call:
mle2(minuslogl = nlldbin, start = st, hessian.opts = TRUE)

Coefficients:
     Estimate  Std. Error z value  Pr(z)
A    -20.4999   5775.1484 -0.0035 0.9972
B     -5.2499 120578.9515  0.0000 1.0000
C     -7.9999 722637.2670  0.0000 1.0000
D     -2.2499  39746.6639 -0.0001 1.0000

> summary(glm2)

Call:
glm(formula = disease ~ ses + sector, family = binomial(link = logit), 
    data = d) 

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.52001    0.33514  -4.535 5.75e-06 ***
ses2        -0.08525    0.41744  -0.204 0.838177    
ses3         0.16086    0.39261   0.410 0.682019    
sector2      1.28098    0.34140   3.752 0.000175 ***

标签： rglmestimationmle

解决方案

这条线

-sum(joint.pdf, log=TRUE ,na.rm=TRUE)

是错的。sum没有特殊的log论点；您正在做的是将值TRUE（转换为 1）添加到 pdf 中。

你想要的是

-sum(log(joint.pdf), na.rm=TRUE)

但由于数字原因，这也不是很好，因为 pdf 可能会下溢。更好的写作方式是

logpdf <- y*log(p) + (1-y)*log(1-p)
-sum(logpdf, na.rm=TRUE)

推荐阅读