0

For a project I'm doing a need to create a LGM from points on a [0, 1]^2 grid. I'd really appreciate any help you all can give me

> set.seed(2000)
> # Define number of points and locations
> par(mfrow=c(1,1))#Ensure 1 graph on screen at a time 
> n <- 10
> s <- matrix(runif(n*2, 0, 1), ncol=2) 
> plot(s)
> 
> # Define mean and covariance matrix for multivariate normal distribution
> mu <- rep(0, n)
> 
> # Create the variance matrix
> r <- 12
> dist_mat <- as.matrix(dist(s))
> rho_ij_mat <- exp(-r*dist_mat)
> Sigma <- rho_ij_mat # Covariance matrix
> 
> #Cholesky
> Cholesky <- chol(Sigma)
> 
> # Generate random values from multivariate normal distribution
> a <- rmvnorm(n, mean=mu, sigma=Sigma)
> plot(a)

> b <- rpois(n, exp(a))
> plot(b) #Our observation
> 
> formula = b ~ 1 + f(a, model = "iid") 
> 
> 
> # Fit the model using INLA
> model <- inla(formula, 
+               family = "poisson", # likelihood family
+               data = data.frame(a, b), # data frame with response and 
+                                        #spatial location variables
+               control.compute = list(dic = TRUE, waic = TRUE)) # compute DIC and WAIC
Error in inla.core(formula = formula, family = family, contrasts = contrasts,  : 
  length(covariate[[r]]) == NPredictor is not TRUE

 *** inla.core.safe:  inla.program has crashed: rerun to get better initial values. try=1/2 
Error in inla.core(formula = formula, family = family, contrasts = contrasts,  : 
  length(covariate[[r]]) == NPredictor is not TRUE
 *** inla.core.safe:  inla.program has crashed: rerun to get better initial values. try=2/2 
Error in inla.core(formula = formula, family = family, contrasts = contrasts,  : 
  length(covariate[[r]]) == NPredictor is not TRUE
Error in inla.core.safe(formula = formula, family = family, contrasts = contrasts,  : 
  *** Fail to get good enough initial values. Maybe it is due to something else.

I've tried converting "a" to a GMRF with limited success as f(a, model = "iid") is the piece of code causing this issue. The only way the formula works is in the form "formula = b ~ 1"

Phil
  • 7,287
  • 3
  • 36
  • 66
user21343597
  • 1
  • 1

1 Answers1

0

Not quite sure what you are trying to estimate. If you want a random effects model tying the values of your multivariate normal to the observations from the Poisson, then you need to flatten the matrices. The data frame input isn't appropriately shaped for that kind of problem. You should have all the y responses in 1 column, and have each covariate in a named column.

user1848065
  • 1
  • 1