Reststandardfehlerdifferenz zwischen optim und glm

16

Ich versuche, mit optimden Ergebnissen einer einfachen linearen Regression mit zu reproduzierenglm oder sogar nlsR-Funktionen ausgestattet ist.
Die Parameterschätzungen sind die gleichen, aber die Restvarianzschätzung und die Standardfehler der anderen Parameter sind nicht die gleichen, insbesondere wenn die Stichprobengröße niedrig ist. Ich nehme an, dass dies auf Unterschiede in der Art und Weise zurückzuführen ist, wie der verbleibende Standardfehler zwischen Maximum Likelihood und Least Square berechnet wird (dividiert durch n oder durch n-k + 1, siehe unten im Beispiel).
Ich verstehe aus meinen Lesungen im Internet, dass Optimierung keine einfache Aufgabe ist, aber ich habe mich gefragt, ob es möglich ist, die Standardfehlerschätzungen aus der glmVerwendung auf einfache Weise zu reproduzieren optim.

Simulieren Sie einen kleinen Datensatz

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) 

Schätzung mit 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

Vergleich mit glm und 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

Ich kann die verschiedenen Standardfehlerschätzungen wie folgt reproduzieren:

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

Antworten:

9

Das Problem ist, dass die Standardfehler von stammen

σ^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

Um mehr als usεr11852- Anforderungen auszuarbeiten , ist die Protokollwahrscheinlichkeit

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

Xn

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

Jetzt können wir entweder den MLE oder den Schätzer von Anhaltspunkt einstecken (siehe Abbildung)σ

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 

Wir können mit einem QR - Zerlegung das gleiche tun wie lmtut

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

Also zu antworten

Ich verstehe aus meinen Lesungen im Internet, dass Optimierung keine einfache Aufgabe ist, aber ich habe mich gefragt, ob es möglich ist, die Standardfehlerschätzungen aus der glmVerwendung auf einfache Weise zu reproduzieren optim.

Dann müssen Sie die Standardfehler in dem von Ihnen verwendeten Gaußschen Beispiel hochskalieren.

Benjamin Christoffersen
quelle
1
+1. Ich bin nicht zu 100% der Meinung, dass Sie es richtig verstanden haben, aber dies ist definitiv in die richtige Richtung. Können Sie erklären, warum Sie diesen Faktor erwarten?
usεr11852 sagt Reinstate Monic
Ist es jetzt klarer?
Benjamin Christoffersen
1
Ja. Gute Antwort! (Ich habe es bereits hochgestuft)
usεr11852 sagt Reinstate Monic
1

optimnn-k+1. Machen Sie also die Teilung durch rückgängign und dividieren durch n-k+1: sqrt(4.717216^2*4/2) = 6.671151

papgeo
quelle
1
Danke für deine Antwort. Mir ist klar, dass meine Frage nicht klar genug war (ich habe sie jetzt bearbeitet). Ich möchte nicht nur die Berechnung der verbleibenden Standardfehler, sondern auch die Parameter der Standardfehler reproduzieren ...
Gilles
@ Gilles Ich weiß nicht, wie ich die Standardfehler reproduzieren soll. Die Unterschiede sind folgende: 1. glm verwendet die Fisher-Informationsmatrix, während das Hessische optimiert wird, und 2. glm betrachtet dies als ein 2-Parameter-Problem (finde b0 und b1), während ein 3-Parameter-Problem optimiert wird (b0, b1 und sigma2). . Ich bin nicht sicher, ob diese Unterschiede überbrückt werden können.
Papgeo