Ich versuche, mit optim
den Ergebnissen einer einfachen linearen Regression mit zu reproduzierenglm
oder sogar nls
R-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 glm
Verwendung 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
quelle
optim
sqrt(4.717216^2*4/2) = 6.671151
quelle