Différence d'erreur standard résiduelle entre optim et glm

16

J'essaie de reproduire avec optimles résultats d'une simple régression linéaire équipée glmou même de nlsfonctions R.
Les estimations des paramètres sont les mêmes, mais l'estimation de la variance résiduelle et les erreurs-types des autres paramètres ne sont pas les mêmes, en particulier lorsque la taille de l'échantillon est faible. Je suppose que cela est dû à des différences dans la façon dont l'erreur standard résiduelle est calculée entre les approches du maximum de vraisemblance et des moindres carrés (en divisant par n ou par n-k + 1, voir ci-dessous dans l'exemple).
Je comprends de mes lectures sur le Web que l'optimisation n'est pas une tâche simple, mais je me demandais s'il serait possible de reproduire de manière simple les estimations d'erreur standard à partir de l' glmutilisation optim.

Simuler un petit ensemble de données

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) 

Estimer avec 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

Comparaison avec glm et 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

Je peux reproduire les différentes estimations d'erreur standard résiduelle comme ceci:

> # 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
Gilles
la source

Réponses:

9

Le problème est que les erreurs standard proviennent de

σ^2(XX)1

σ^2summary.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    ...

(β0,β1)σ^2(β0,β1,σ)σn/(n3+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

Pour en savoir plus sur les demandes usεr11852 , la log-vraisemblance est

l(β,σ)=n2log(2π)nlogσ12σ2(yXβ)(yXβ)

Xn

-ββl(β,σ)=1σ2XX

σ

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 

On peut faire la même chose avec une décomposition QR comme le lmfait

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

Donc pour répondre

Je comprends de mes lectures sur le Web que l'optimisation n'est pas une tâche simple, mais je me demandais s'il serait possible de reproduire de manière simple les estimations d'erreur standard à partir de l' glmutilisation optim.

vous devez ensuite mettre à l'échelle les erreurs standard dans l'exemple gaussien que vous utilisez.

Benjamin Christoffersen
la source
1
+1. Je ne suis pas sûr à 100% que vous l'ayez entièrement compris, mais c'est certainement dans la bonne direction. Pouvez-vous expliquer pourquoi vous vous attendez à ce facteur?
usεr11852 dit Reinstate Monic
Est-ce plus clair maintenant?
Benjamin Christoffersen
1
Oui. Bonne réponse! (Je l'ai déjà voté)
usεr11852 dit Reinstate Monic
1

optimnn-k+1nn-k+1: sqrt(4.717216^2*4/2) = 6.671151

papgeo
la source
1
Merci pour votre réponse. Je me rends compte que ma question n'était pas assez claire (je l'ai maintenant éditée). Je ne veux pas seulement reproduire le calcul de l'erreur standard résiduelle mais aussi les paramètres des erreurs standard ...
Gilles
@Gilles Je ne sais pas reproduire les erreurs standard. Les différences sont dues à: 1. glm utilise la matrice d'informations de Fisher, tout en optimisant la toile de jute, et 2. glm considère cela comme un problème à 2 paramètres (trouver b0 et b1), tandis qu'optim un problème à 3 paramètres (b0, b1 et sigma2) . Je ne sais pas si ces différences peuvent être comblées.
papgeo