[R] hessian in solnp
CAMARDA Carlo Giovanni
c@r|o-g|ov@nn|@c@m@rd@ @end|ng |rom |ned@|r
Thu Apr 28 00:16:05 CEST 2022
In the commented snippet below, a behavior of the solnp() I cannot fully explain: when a constraint is added, hessian matrix is obviously changing, but in a way I don't understand. I know that the model below could be estimated differently and in a much efficient way. However I wanted to present the simplest example I could think.
In particular, I encountered this issue because I aim to analytically compute variance-covariance matrix in a (non-linear) model with equality constraints. Any practical hint to find a solution to this more general issue would be also welcome.
Thanks, GC
## simulating a simple linear model
n <- 100
x <- sort(runif(n))
X <- cbind(1,x)
k <- ncol(X)
y <- X%*%c(3,4) + rnorm(n, sd=0.5)
## estimating with lm()
fitlm <- lm(y~x)
## using solnp
library(Rsolnp)
## objecting function
fn1 <- function(par){
y.hat <- X%*%par
res <- y-y.hat
RSS <- t(res)%*%res/2
return(c(RSS))
}
## using solnp without any constraint
solUR <- solnp(c(3.2, 3.5), fun = fn1)
## getting the same fitted objects as in lm()
## coefficients
cbind(solUR$pars, fitlm$coef)
## and standard errors by variance-covariance matrix computed
## with hessian internally estimated
(H.UR <- solUR$hessian)
y.UR <- X%*%solUR$pars
RSS.UR <- t(y-y.UR)%*%(y-y.UR)
sigma2.UR <- RSS.UR/(n-k)
V.UR <- c(sigma2.UR)*solve(H.UR)
se.UR <- sqrt(diag(V.UR))
cbind(se.UR,summary(fitlm)$coef[,2])
## adding equality constraint
eqn1 <- function(par) sum(par)
c <- sum(solUR$pars)
## I use the sum of the previously estimated coefficients
solR <- solnp(c(3.2, 3.5), fun = fn1, eqfun = eqn1, eqB = c)
## as expected we get the same coefficients
cbind(solR$pars,fitlm$coef)
## but the hessian is rather different
(H.R <- solR$hessian)
## which leads to completely different (much larger here)
## standard errors
y.R <- X%*%solR$pars
RSS.R <- t(y-y.R)%*%(y-y.R)
sigma2.R <- RSS.R/(n-k)
V.R <- c(sigma2.R)*solve(H.R)
se.R <- sqrt(diag(V.R))
cbind(se.R,summary(fitlm)$coef[,2])
## I noticed that differences in the hessian
## is ~constant over both dimensions
## and it depends (also) on sample size
H.R-H.UR
[[alternative HTML version deleted]]
More information about the R-help
mailing list