#1.Load data & packages
#Load libraries
library(INLA)
library(spatialreg)
library(parallel)

#Load data
data(boston)
#library(maptools)
library(sf)
library(tidyverse)

#readShapePoly(system.file("etc/shapes/boston_tracts.shp",
#  package="spdep")[1], ID="poltract",
#  proj4string=CRS(paste("+proj=longlat +datum=NAD27 +no_defs +ellps=clrk66",
#  "+nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat")))
boston_nb <- spdep::poly2nb(boston.tr)
censored <- boston.tr$CMEDV == 50

#Set censored areas to NA for spatial prediction
boston.c <- boston.tr
boston.c$CMEDV[censored]<-NA
lw <- spdep::nb2listw(boston_nb, style="W")

#2 Create an index
boston.c$idx <- 1:n

#3 duplicate data 
boston1012.c <- rbind(boston.c, boston.c)
boston1012.tr <- rbind(boston.tr, boston.tr)

# so now I have 1012 data points, 2 for each area. 

4# W matrix 
#(the index is already created index in step 2)
# W matrix 
#(the index is already created index in step 2)
#Define some indices used in the models
#n <- nrow(boston.c)
#boston.c$idx <- 1:n

#Adjcency matrix
#W<-nb2mat(boston.soi)
W <- as(as_dgRMatrix_listw(lw), "CsparseMatrix")
#5 Create mmatrix 
#Model matrix for SLM models
f1 <- log(CMEDV) ~ CRIM 

#+ ZN + INDUS + CHAS + I(NOX^2)+ I(RM^2) +  AGE + log(DIS) + log(RAD) + TAX + PTRATIO + B + log(LSTAT)

 
mmatrix <- model.matrix(f1, boston1012.tr)#NOTE: We use boston.tr because NA's in boston.c


#5 priors
#DEFINE PRIORS TO BE USED WITH R-INLA

#Zero-variance for error term
zero.variance = list(prec = list(initial = 15, fixed = TRUE))

#Compute eigenvalues for SLM model (as in Havard's code)

re.idx <- which(abs(Im(e)) < 1e-6)
rho.max <- 1 / max(Re(e[re.idx]))
rho.min <- 1 / min(Re(e[re.idx]))
rho <- mean(c(rho.min, rho.max))


#
#Variance-covariance matrix for beta coeffients' prior
#
betaprec <- 0.0001
#Standard regression model
Q.beta <- Diagonal(n = ncol(mmatrix), x = 1)
Q.beta <- betaprec * Q.beta 

#This defines the range for the spatial autocorrelation parameters
# the spatial adjacenty matrix W, the associated
#matrix of covariates X and the precision matrix Q.beta for the prior
#on the coefficients of the covariates 
#
args.slm <- list(
  rho.min = rho.min,
  rho.max = rho.max,
  W = W,
  X = matrix(0, nrow(mmatrix), 0),
  Q.beta = matrix(1, 0, 0)
)

#Definition of priors for precision of the error effect in the slm latent
#effect and the spatial autocorrelation parameter (in the (0,1) interval).
#
hyper.slm <- list(
  prec = list(
    prior = "loggamma", param = c(0.01, 0.01)),
    rho = list(initial=0, prior = "logitbeta", param = c(1, 1))
)


#6 run model with boston1012.c dataset with 1012 records and 506 areas (idx)

#6 run model with Boston2 dataset with 1012 records and 506 areas (idx)

#SLM model
slmm1<-inla(log(CMEDV) ~ -1 + 
   f(idx, model="slm", 
      args.slm=list(rho.min = rho.min, rho.max = rho.max, W=W, X=mmatrix, 
         Q.beta=Q.beta), 
      hyper=hyper.slm),
   data=as.data.frame(boston1012.c), family="gaussian",
   control.family = list(hyper=zero.variance),
   control.compute=list(dic=TRUE, cpo=TRUE)
)

END

Throws off an error
Error in INLA::f(idx, model = "slm", args.slm = list(rho.min = rho.min, : all(dim(args.slm$W) == c(slm.n, slm.n)) is not TRUE *** inla.core.safe: inla.program has crashed: rerun to get better initial values. try=1/1 Error in INLA::f(idx, model = "slm", args.slm = list(rho.min = rho.min, : all(dim(args.slm$W) == c(slm.n, slm.n)) is not TRUE
Error in inla.core.safe(formula = formula, family = family, contrasts = contrasts, : *** Failed to get good enough initial values. Maybe it is due to something else. Error in inla.core.safe(formula = formula, family = family, contrasts = contrasts, : *** Failed to get good enough initial values. Maybe it is due to someth



ICAR model 
As a side point, the above hierarchical data set up works for an ICAR framework within INLA
form <- log(CMEDV) ~ 1 + CRIM

inla_model <- inla(
  update(form, . ~ . + f(idx, model = "besag", graph = W)), # Update formula with additional terms (to_check)
  family = "gaussian", # Specify the family for the model
  data = boston1012, # Provide the dataset
  control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE),
  control.predictor = list(compute = TRUE)
)

summary(inla_model)

inla_model$summary.random$idx$mean


boston_506$inla_model <- inla_model$summary.random$idx$mean # to_check: need to make joins on area code rather than cbind/ appending a column

boston_506_drop <- st_drop_geometry(boston_506)

boston_506_drop %>% head(10) %>% display() # Display for user