Hessian is not positive definite when using UTM coordinates
Opened this issue · 0 comments
elizabethng commented
Here is a reproducible example of an issue I had when fiting a spatio-temporal index standardization model to UTM-georeferenced data. If I make my spde mesh using UTMs, the Hessian is not positive definite. Scaling the location data seems to fix the problem. The example below uses the .cpp file from the Week 7 Lecture.
Running with UTM = TRUE
and RESCALE = FALSE
will give the error message
Warning message: In TMBhelper::Optimize(obj = Obj, newtonsteps = 1, getsd = TRUE, : Hessian is not positive definite, so standard errors are not available
Otherwise, running with UTM = FALSE
or with RESCALE = TRUE
(which re-scales the coordinates) will run fine.
setwd("~/GitHub/2018_FSH556/Week 7 -- spatiotemporal models/Lecture")
## Libraries and functions
library( TMB )
library( INLA )
library( RandomFields )
library( RANN )
## Simulate Data
Sim_Fn = function( UTM = FALSE, logmean=1, Scale=0.2, SD_omega=1, SD_epsilon=1, n_per_year=100, n_years=10 ){
loc_xy = cbind( "x"=runif(n_per_year), "y"=runif(n_per_year))
if(UTM == TRUE) loc_xy = 10^6*loc_xy # UTM-like scale
RF_omega = RMgauss(var=SD_omega^2, scale=Scale)
RF_epsilon = RMgauss(var=SD_epsilon^2, scale=Scale)
log_d_rt = Epsilon_rt = matrix(NA, ncol=n_years, nrow=n_per_year)
Omega_r = RFsimulate(model=RF_omega, x=loc_xy[,'x'], y=loc_xy[,'y'])@data[,1]
for(t in 1:n_years){
Epsilon_rt[,t] = RFsimulate(model=RF_epsilon, x=loc_xy[,'x'], y=loc_xy[,'y'])@data[,1]
log_d_rt[,t] = logmean + Epsilon_rt[,t] + Omega_r
}
b_t = colSums( exp(log_d_rt) )
# Simulate samples for each site and year
DF = expand.grid("s_i"=1:n_per_year, "t_i"=1:10)
DF = cbind( DF, "log_d_rt"=log_d_rt[ as.matrix(DF[,c('s_i','t_i')]) ] )
DF = cbind( DF, "c_i"=rpois(n_per_year*n_years, lambda=exp(DF[,'log_d_rt'])) )
# Return stuff
Return = list("DF"=DF, "Epsilon_rt"=Epsilon_rt, "Omega_r"=Omega_r, "loc_xy"=loc_xy, "b_t"=b_t, "n_per_year"=n_per_year, "n_years"=n_years)
return( Return )
}
UTM = TRUE # Simulate data with UTM-like coordinates
RESCALE = FALSE # Don't fix the problem
SimList = Sim_Fn(UTM = UTM)
## Fit the Model
# Area for each location
loc_extrapolation = expand.grid( "x"=seq(0,1,length=1e3), "y"=seq(0,1,length=1e3) )
if(UTM == TRUE) loc_extrapolation = loc_extrapolation*10^6
NN_extrapolation = nn2( data=SimList$loc_xy, query=loc_extrapolation, k=1 )
a_s = table(factor(NN_extrapolation$nn.idx,levels=1:nrow(SimList$loc_xy))) / nrow(loc_extrapolation)
# Scale down
if(UTM == TRUE & RESCALE == TRUE) SimList$loc_xy = SimList$loc_xy*10^-6
# Make triangulated mesh
mesh = inla.mesh.create( SimList$loc_xy )
spde = inla.spde2.matern(mesh)
# Compile
Version = "spatial_index_model"
compile( paste0(Version,".cpp") )
dyn.load( dynlib(Version) )
# Make inputs
Data = list("n_s"=SimList$n_per_year, "n_t"=SimList$n_years, "a_s"=a_s, "c_i"=SimList$DF$c_i, "s_i"=SimList$DF$s_i-1, "t_i"=SimList$DF$t_i-1, "M0"=spde$param.inla$M0, "M1"=spde$param.inla$M1, "M2"=spde$param.inla$M2)
Params = list("beta0"=0, "ln_tau_O"=log(1), "ln_tau_E"=log(1), "ln_kappa"=1, "omega_s"=rep(0,mesh$n), "epsilon_st"=matrix(0,nrow=mesh$n,ncol=Data$n_t))
Random = c("omega_s", "epsilon_st")
# Build and run
Obj = MakeADFun( data=Data, parameters=Params, random=Random, DLL = paste(Version) )
Opt = TMBhelper::Optimize( obj=Obj, newtonsteps=1, getsd=TRUE, bias.correct=TRUE )