change the value of zero.variance in spatial error model case lead INLA crashes

[shubingtang]

I'm trying to simulate a spatially autocorrelated random data, and I draw both co-ordinates of 40 points,by using Voronoi tessellation to construct units and get rook type spatial weight matrix w（row-standardisation ). #generating data N=40 set.seed(1) beta=c(1,3) rho=0.7 alpha=1 x1=rnorm(N,40,2) x2=rnorm(N,36,5) Irhow.t=Diagonal(nrow(w))-rho*w Irhow2.t=crossprod(Irhow.t) sigma=1 sigma2=sigma*solve(Irhow2.t) epsilon=mvrnorm(1,rep(0,N),sigma2) y=solve(Irhow.t)%*%(x1*beta[1]+x2*beta[2]+alpha)+epsilon

I don't konw what's the zero.variance use for ,so I try to change its value,the simulation either crashes (INLA crashes in early iteration xxx)or its result not working well，even I start point that very close to true values.

Thank you for any help in advance.

[becarioprecario]

The zero.variance stuff is used to set the variance of the Gaussian error to zero. In this way, we remove the Gaussian error.

El 24 abr 2019, a las 9:04, shubingtang notifications@github.com<mailto:notifications@github.com> escribió:

I'm trying to simulate a spatially autocorrelated random data, and I draw both co-ordinates of 40 points,by using Voronoi tessellation to construct units and get rook type spatial weight matrix w（row-standardisation ).

generating data N=40 set.seed(1) beta=c(1,3) rho=0.7 alpha=1 x1=rnorm(N,40,2) x2=rnorm(N,36,5) Irhow.t=Diagonal(nrow(w))-rhow Irhow2.t=crossprod(Irhow.t) sigma=1 sigma2=sigmasolve(Irhow2.t) epsilon=mvrnorm(1,rep(0,N),sigma2) y=solve(Irhow.t)%%(x1beta[1]+x2*beta[2]+alpha)+epsilon

I don't konw what's the zero.variance use for ,so I try to change its value,the simulation either crashes (INLA crashes in early iteration xxx)or its result not working well，even I start point that very close to true values.

Thank you for any help in advance.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHubhttps://github.com/becarioprecario/INLAMCMC_examples/issues/1, or mute the threadhttps://github.com/notifications/unsubscribe-auth/ABYD6WQPKSE4XODTJJF3EFTPSABADANCNFSM4HIA6I3A.

[shubingtang]

Thank you very much,I already understand zero.variance is used for additional Gaussian error.So I keep the inital value as 25.But still "Error in inla.inlaprogram.has.crashed() ",and I try to rerun my code ,the INLA error number 1 happened at differnt iteration every time.

[shubingtang]

I turn on verbose=True,and I get this"GMRFLib version 3.0-0-snapshot, has recived error no [21] Reason : This should not happen Message : Condition `density->Pinv' is not TRUE Function : GMRFLib_density_Pinv File : density.c Line : 1036 RCSId : file: density.c hgid: b84cb08f11c9 date: Thu Jul 12 13:38:25 2018 +0300"

where did I go wrong?

[becarioprecario]

You need to provide a fully reproducible example. Otherwise, it is imposible to say what is going on.

[shubingtang]

You need to provide a fully reproducible example. Otherwise, it is imposible to say what is going on.

Sorry,there is my code: `#construct spatial weight matrix library(sp) library(deldir) library(spdep) voronoipolygons = function(layer) { require(deldir) crds = layer@coords z = deldir(crds[,1], crds[,2]) w = tile.list(z) polys = vector(mode='list', length=length(w)) require(sp) for (i in seq(along=polys)) { pcrds = cbind(w[[i]]$x, w[[i]]$y) pcrds = rbind(pcrds, pcrds[1,]) polys[[i]] = Polygons(list(Polygon(pcrds)), ID=as.character(i)) } SP = SpatialPolygons(polys) voronoi = SpatialPolygonsDataFrame(SP, data=data.frame(dummy = seq(length(SP)), row.names=sapply(slot(SP, 'polygons'), function(x) slot(x, 'ID')))) }

create some points

set.seed(1) N=100 x.coord=runif(N,0,200) y.coord=runif(N,0,200) points=data.frame(Longitude=x.coord,Latitude =y.coord) coordinates(points) <- ~ Longitude + Latitude

get Voronoi polygons

vp=voronoipolygons(points) temp=poly2nb(vp,queen = FALSE) w.lw=nb2listw(temp,style = 'W') w=listw2mat(w.lw)

############# library(INLA) library(INLABMA) library(MASS)

generating data

beta=c(1,3) rho=0.7 alpha=1 x1=rnorm(N,40,2) x2=rnorm(N,17,5) Irhow1=Diagonal(nrow(w))-rhow sigma=0.09 epsilon=rnorm(N,0,sigma) epsilon=solve(Irhow1)%%epsilon y=solve(Irhow1)%%(x1beta[1]+x2*beta[2]+alpha)+epsilon y=as.vector(y) data=data.frame(x1=x1,x2=x2,y=y,idx=1:N)

Formula

formula=y~x1+x2

zero.variance = list(prec=list(initial = 25, fixed=TRUE))

Fit model for a given rho

fit.inla <- function(data, rho) { res=slm.inla(formula,d=data,W=w,rho=rho,family="gaussian", impacts = FALSE,verbose = FALSE, control.family = list(hyper = zero.variance))

return(list(mlik = res$mlik[1,1], model = res)) }

Proposal distribution x -> y

density

dq.rho<- function(x, y, sigma = .15, log =TRUE) { res <- dnorm(y, mean = x, sd = sigma, log = log) }

random

rq.rho <- function(x, sigma = .15) { rnorm(1, mean = x, sd = sigma) }

Prior for rho Uniform [-0.6, 1] (based on the eigenvalues of W)

prior.rho <- function(x, log = TRUE) { res <- dunif(x, -0.6, 1, log = log) }

Data as 'd'

d=data

Run INLA within MCMC

inlamh.res <- INLAMH(d, fit.inla, 0.5, rq.rho, dq.rho, prior.rho, n.sim = 200, n.burnin = 20, n.thin = 1, verbose =FALSE)`

[becarioprecario]

I get:

Irhow1=Diagonal(nrow(w))-rhow Error in eval(ei, envir) : objeto 'rhow' no encontrado

Please, make sure that the code that you provide runs smoothly.

[shubingtang]

I just copied my code and didn't notice all '*' are disappeared in this page,I am so sorry about that.

slm.zip

I am so sorry again and thank you for your time.

[becarioprecario]

Thanks. ‘rho’ should be between -0.6 and 1, and the proposal distribution does not take care of that. You may want to change this or, alternatively, set sing a to something smaller in rq and dq. 0.05 seems to be fine to produce small ’steps’ in the values of rho.