Diferença residual de erro padrão entre optim e glm

16

Eu tento reproduzir com optimos resultados de uma regressão linear simples ajustada com glmou até nlsfunções R.
As estimativas de parâmetros são as mesmas, mas a estimativa de variação residual e os erros padrão dos outros parâmetros não são os mesmos, principalmente quando o tamanho da amostra é baixo. Suponho que isso ocorra devido a diferenças na maneira como o erro padrão residual é calculado entre as abordagens de Máxima Verossimilhança e Menos Quadrados (dividindo por n ou por n-k + 1, veja abaixo no exemplo).
Compreendo pelas minhas leituras na Web que a otimização não é uma tarefa simples, mas fiquei pensando se seria possível reproduzir de maneira simples o erro padrão estimado glmdurante o uso optim.

Simule um pequeno conjunto de dados

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) 

Estimar com otimização

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

Comparação com glm e 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

Posso reproduzir as diferentes estimativas de erro padrão residual como este:

> # 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
fonte

Respostas:

9

O problema é que os erros padrão vêm 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 0,β1)σ^2(β0 0,β1,σ)σ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

Para elaborar mais solicitações de usεr11852 , a probabilidade de log é

eu(β,σ)=-n2registro(2π)-nregistroσ-12σ2(y-Xβ)(y-Xβ)

Xn

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

Agora podemos conectar o MLE ou o estimador não-avaliado de σ

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 

Podemos fazer o mesmo com uma decomposição QR como lmfaz

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

Então, para responder

Compreendo pelas minhas leituras na Web que a otimização não é uma tarefa simples, mas fiquei pensando se seria possível reproduzir de maneira simples o erro padrão estimado glmdurante o uso optim.

então você precisa aumentar os erros padrão no exemplo gaussiano que você usa.

Benjamin Christoffersen
fonte
1
+1. Não estou 100% certo de que você acertou totalmente, mas isso definitivamente está na direção correta. Você pode explicar por que espera esse fator?
usεr11852 diz Reinstate Monic
Está mais claro agora?
Benjamin Christoffersen
1
Sim. Boa resposta! (Eu já votei)
usεr11852 diz Reinstate Monic
1

optimnnk+1nnk+1sqrt(4.717216^2*4/2) = 6.671151

papgeo
fonte
1
Obrigado pela sua resposta. Percebo que minha pergunta não foi clara o suficiente (agora a editei). Eu não só deseja reproduzir o cálculo residual padrão de erro, mas também os parâmetros erros padrão ...
Gilles
@ Gilles Não sei como reproduzir os erros padrão. As diferenças são devido a: 1. glm usa a matriz de informações de Fisher, enquanto otimiza o hessian, e 2. glm considera isso um problema de 2 parâmetros (encontre b0 e b1), enquanto otimiza um problema de 3 parâmetros (b0, b1 e sigma2) . Não tenho certeza se essas diferenças podem ser superadas.
Papgeo 13/0818