J'essaie de reproduire avec optim
les résultats d'une simple régression linéaire équipée glm
ou même de nls
fonctions 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' glm
utilisation 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
la source
optim
sqrt(4.717216^2*4/2) = 6.671151
la source