I'm using the solnp()
function in the R package Rsolnp to solve a nonlinear regression with constraints. It works well, converges with no problem. I want to use the Hessian matrix to calculate standard errors of the four parameter estimates, but the Hessian is not 4 by 4 as I had expected, but 5 by 5. I looked around on SO and didn't see anyone else with an unexpected Hessian size. All the examples I found with the Hessians printed showed them to be the expected size of p by p (e.g., 2x2, 3x3, and 4x4).
How can I get standard errors for my 4 parameters from this 5 by 5 Hessian?
df <- data.frame(
Recruit.N = c(78.4, 79.8, 106, 57.4, 81.7, 94.4, 74.1, 42, 61.6, 47.7, 61.8,
28.1, 32.3, 19, 23.4, 20.1, 27),
Stock.5 = c(66.6, 90.3, 138.5, 79.8, 77.3, 78.4, 79.8, 106, 57.4, 81.7, 94.4,
74.1, 42, 61.6, 47.7, 61.8, 28.1),
Stock.6 = c(25.2, 66.6, 90.3, 138.5, 79.8, 77.3, 78.4, 79.8, 106, 57.4, 81.7,
94.4, 74.1, 42, 61.6, 47.7, 61.8),
Stock.7 = c(23.8, 25.2, 66.6, 90.3, 138.5, 79.8, 77.3, 78.4, 79.8, 106, 57.4,
81.7, 94.4, 74.1, 42, 61.6, 47.7)
)
lossfcn <- function(parz, mydat) {
alpha <- parz[[1]]
beta <- parz[[2]]
p5 <- parz[[3]]
p6 <- parz[[4]]
p7 <- 1 - p5 - p6
S <- with(mydat, p5*Stock.5 + p6*Stock.6 + p7*Stock.7)
Obs <- mydat$Recruit.N
Pred <- alpha * S * exp(-beta*S)
Resid <- log(Obs) - log(Pred)
sigma <- sqrt(mean(Resid^2))
LL <- dlnorm(Obs, meanlog=log(Pred), sdlog=sigma, log=TRUE)
-sum(LL)
}
inequal <- function(parz, mydat) {
parz[3] + parz[4]
}
library(Rsolnp)
solnp(pars=c(1, 0.008, 1/3, 1/3), fun=lossfcn, mydat=df,
ineqfun=inequal, ineqLB=0, ineqUB=1,
LB=c(0, 0, 0, 0), UB=c(1000, 1000, 1, 1), control=list(trace=0))
$pars
[1] 6.731317e-01 1.888572e-10 8.141363e-01 1.858631e-01
$convergence
[1] 0
$values
[1] 79.87150 75.50927 75.50927 75.50927
$lagrange
[,1]
[1,] -2.028222
$hessian
[,1] [,2] [,3] [,4] [,5]
[1,] 0.3350868 -3.359077e-01 17.84919 -0.4306057 -0.3382811
[2,] -0.3359077 1.993956e+02 -10161.63351 -7.0844295 -2.2749785
[3,] 17.8491854 -1.016163e+04 548099.69224 -85.9544831 -224.0362766
[4,] -0.4306057 -7.084429e+00 -85.95448 25.1086694 5.8817704
[5,] -0.3382811 -2.274979e+00 -224.03628 5.8817704 4.1978178
$ineqx0
[1] 0.9999995
$nfuneval
[1] 142
$outer.iter
[1] 3
$elapsed
Time difference of 0.03016496 secs
$vscale
[1] 1 1 1 1 1 1