Hessian is not positive definite when using UTM coordinates

[elizabethng]

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 )