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]]
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
