optim和glm之间的剩余标准误差差异

机器算法验证 r 最大似然 优化
2022-01-25 18:06:07

我尝试使用拟合甚至R 函数optim的简单线性回归的结果来重现。 参数估计值相同，但残差估计值和其他参数的标准误不同，特别是在样本量较小时。我想这是由于最大似然法和最小二乘法之间计算剩余标准误差的方式存在差异（除以 n 或 n-k+1 见下文示例）。 我从网上的阅读中了解到优化不是一项简单的任务，但我想知道是否有可能以简单的方式重现使用.glmnls

glmoptim

模拟一个小数据集

set.seed(1)
n = 4 # very small sample size !
b0 <- 5
b1 <- 2
sigma <- 5
x <- runif(n, 1, 100)
y =  b0 + b1*x + rnorm(n, 0, sigma)

用 optim 估计

negLL <- function(beta, y, x) {
    b0 <- beta[1]
    b1 <- beta[2]
    sigma <- beta[3]
    yhat <- b0 + b1*x
    likelihood <- dnorm(y, yhat, sigma)
    return(-sum(log(likelihood)))
}

res <- optim(starting.values, negLL, y = y, x = x, hessian=TRUE)
estimates <- res$par     # Parameters estimates
se <- sqrt(diag(solve(res$hessian))) # Standard errors of the estimates
cbind(estimates,se)


    > cbind(estimates,se)
      estimates         se
b0     9.016513 5.70999880
b1     1.931119 0.09731153
sigma  4.717216 1.66753138

与 glm 和 nls 的比较

> m <- glm(y ~ x)
> summary(m)$coefficients
            Estimate Std. Error   t value    Pr(>|t|)
(Intercept) 9.016113  8.0759837  1.116411 0.380380963
x           1.931130  0.1376334 14.030973 0.005041162
> sqrt(summary(m)$dispersion) # residuals standard error
[1] 6.671833
> 
> summary(nls( y ~ b0 + b1*x, start=list(b0 = 5, b1= 2)))

Formula: y ~ b0 + b1 * x

Parameters:
   Estimate Std. Error t value Pr(>|t|)   
b0   9.0161     8.0760   1.116  0.38038   
b1   1.9311     0.1376  14.031  0.00504 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.672 on 2 degrees of freedom

我可以像这样重现不同的残差标准误差估计：

> # optim / Maximum Likelihood estimate
> sqrt(sum(resid(m)^2)/n)
[1] 4.717698
> 
> # Least squares estimate (glm and nls estimates)
> k <- 3 # number of parameters
> sqrt(sum(resid(m)^2)/(n-k+1))
[1] 6.671833
2个回答

问题是标准错误来自

σ^2(X⊤X)−1

在哪里σ^2是无偏估计量，而不是 MLE。看summary.lm

summary.lm
#R function (object, correlation = FALSE, symbolic.cor = FALSE, 
#R     ...) 
#R {
#R    z <- object
#R    p <- z$rank
#R    rdf <- z$df.residual
#R    ...
#R    Qr <- qr.lm(object) 
#R    ... 
#R    r <- z$residuals
#R    f <- z$fitted.values
#R    w <- z$weights
#R    if (is.null(w)) {
#R         mss <- if (attr(z$terms, "intercept")) 
#R             sum((f - mean(f))^2)
#R         else sum(f^2)
#R         rss <- sum(r^2)
#R    }
#R    ...
#R    resvar <- rss/rdf
#R    ...
#R    R <- chol2inv(Qr$qr[p1, p1, drop = FALSE])
#R    se <- sqrt(diag(R) * resvar)
#R    ...

这是逆观察到的 Fisher信息(β0,β1)有条件的σ^2. 现在，您计算的反向观察到的 Fisher 信息是针对三元组的(β0,β1,σ). 即，您使用的 MLEσ而不是无偏估计量。因此，我认为标准误差应该因因素而异n/(n−3+1)或类似的东西。这是这种情况

set.seed(1)
n = 4 # very small sample size !
b0 <- 5
b1 <- 2
sigma <- 5
x <- runif(n, 1, 100)
y =  b0 + b1*x + rnorm(n, 0, sigma) 

negLL <- function(beta, y, x) {
  b0 <- beta[1]
  b1 <- beta[2]
  sigma <- beta[3]
  yhat <- b0 + b1*x
  return(-sum(dnorm(y, yhat, sigma, log = TRUE)))
}

res <- optim(c(0, 0, 1), negLL, y = y, x = x, hessian=TRUE)
estimates <- res$par     # Parameters estimates
(se <- sqrt(diag(solve(res$hessian))))
#R [1] 5.690 0.097 1.653
k <- 3
se * sqrt(n / (n-k+1))
#R [1] 8.047 0.137 2.338

为了详细说明usεr11852的要求，对数似然是

l(β→,σ)=−n2log⁡(2π)−nlog⁡σ−12σ2(y→−Xβ→)⊤(y→−Xβ→)

在哪里X是设计矩阵和n是观察次数。因此，观察到的信息矩阵是

−∇β→∇β→⊤l(β→,σ)=1σ2X⊤X

现在我们可以插入 MLE 或σ如下图

m <- lm(y ~ x)
X <- cbind(1, x)
sqrt(sum(resid(m)^2)/n       * diag(solve(crossprod(X))))
#R                     x 
#R 5.71058285 0.09732149
k <- 3
sqrt(sum(resid(m)^2)/(n-k+1) * diag(solve(crossprod(X))))
#R                   x 
#R 8.0759837 0.1376334

我们可以对 QR 分解做同样的lm事情

obj <- qr(X)
sqrt(sum(resid(m)^2)/(n-k+1) * diag(chol2inv(obj$qr)))
#R [1] 8.0759837 0.1376334

所以来回答

我从网上的阅读中了解到优化不是一项简单的任务，但我想知道是否有可能以简单的方式重现glm使用optim.

那么您需要在您使用的高斯示例中扩大标准误差。

如果我理解得很好，解决方案很简单：optim通过将残差平方和除以最大化可能性n. 你想要的是将平方和除以n−k+1. 所以撤消除法n并除以n−k+1：sqrt(4.717216^2*4/2) = 6.671151

其它你可能感兴趣的问题
上一篇使用批量标准化时是否需要缩放数据 [0,1]？ 下一篇梯度提升是否适合事件率低（如 1%）的数据？