Różnica błędów resztkowych standardowych między optym a glm

16

Staram się odtworzyć z optimwynikami prostej regresji liniowej zaopatrzonej glmlub nawet nlsfunkcje R.
Oszacowania parametrów są takie same, ale oszacowanie wariancji rezydualnej i błędy standardowe innych parametrów nie są takie same, szczególnie gdy wielkość próby jest niska. Przypuszczam, że jest to spowodowane różnicami w sposobie obliczania resztkowego błędu standardowego między podejściami maksymalnego prawdopodobieństwa i najmniejszych kwadratów (dzielenie przez n lub n-k + 1 patrz poniżej w przykładzie).
Z moich odczytów w sieci rozumiem, że optymalizacja nie jest prostym zadaniem, ale zastanawiałem się, czy można w prosty sposób odtworzyć standardowe oszacowania błędów glmpodczas używania optim.

Symuluj mały zestaw danych

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) 

Oszacuj z optym

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

Porównanie z glm i 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

Mogę odtworzyć różne szacunkowe wartości błędu standardowego w następujący sposób:

> # 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
źródło

Odpowiedzi:

9

Problem polega na tym, że pochodzą ze standardowych błędów

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

Aby rozwinąć więcej w żądaniach usεr11852 , prawdopodobieństwo dziennika jest następujące

l(β,σ)=-n2)log(2)π)-nlogσ-12)σ2)(y-Xβ)(y-Xβ)

Xn

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

σ

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 

Możemy zrobić to samo z rozkładem QR, co lmrobi

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

Tak, aby odpowiedzieć

Z moich odczytów w sieci rozumiem, że optymalizacja nie jest prostym zadaniem, ale zastanawiałem się, czy można w prosty sposób odtworzyć standardowe oszacowania błędów glmpodczas używania optim.

następnie musisz przeskalować standardowe błędy w użytym przykładzie Gaussa.

Benjamin Christoffersen
źródło
1
+1. Nie jestem w 100% przekonany, że masz to w pełni poprawne, ale zdecydowanie jest to właściwy kierunek. Czy możesz wyjaśnić, dlaczego oczekujesz tego czynnika?
usεr11852 mówi: Przywróć Monic
Czy teraz jest bardziej jasne?
Benjamin Christoffersen
1
Tak. Dobra odpowiedź! (Już go głosowałem)
usεr11852 mówi: Przywróć Monic w dniu
1

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

papgeo
źródło
1
Dzięki za odpowiedź. Zdaję sobie sprawę, że moje pytanie nie było wystarczająco jasne (właśnie je edytowałem). Chcę nie tylko odtworzyć resztkowe obliczenia błędów standardowych, ale także błędy standardowe parametrów ...
Gilles
@Gilles Nie wiem, jak odtworzyć standardowe błędy. Różnice wynikają z: 1. glm używa macierzy informacji Fishera, podczas gdy optymalizuje hessian, i 2. glm uważa to za problem 2 parametrów (znajdź b0 i b1), podczas gdy optymalizuje problem 3 parametrów (b0, b1 i sigma2) . Nie jestem pewien, czy te różnice można pokonać.
papgeo