INLA graph

82 views
Skip to first unread message

Pietro Maietta

unread,
Jan 18, 2024, 12:45:26 PM1/18/24
to R-inla discussion group
good evening, I have a question regarding the congruity matrix, starting from a shape file of southern Brazil, a double congruity is processed with the function nb2INLA("nbgINLA.graph",pr.nb), why, I modified the isolated points and removed all the NAs? then I also have difficulties in the formula, I made some changes to the dataframe, and once in a while when viewing in debug it returns some values and then returns the usual error message for missing values.the code is missing the final part which I haven't inserted yet since I can't move forward.
Thanks for the attention
plot_zoom_png.pngplot_zoom_png.png
#####5.8Spatial Car Accident Insurance Analysis
NbMean<-function(shp,vari){
  library(spdep)
  shpnb<-poly2nb(shp)
  shpnb.mat<-nb2mat(shpnb,style="B",zero.policy=TRUE)#matrice adiacente
  selNA<-which(is.na(shp@data[,vari]))
  NAnb<-shpnb.mat[selNA,]
  shp@data[selNA,vari]<-apply(NAnb,1,FUN=function(x)mean(shp@data[which(x==1),vari],na.rm=TRUE))
return(shp)
}

#####Modello INLA (Integrated Nested Laplace Approximations).

library(maptools)
library(INLA)
library(fmesher)
library(Rgraphviz)
library(spdep)
library(igraph)

#source("http://www.math.ntnu.no/inla/givemeINLA.R")
#update(inla.version())
shape<-readShapeSpatial("sul+sp_shape")
#text(coordinates(shape),label=shape@data$City,cex=.3)#aggiunge il nome

pr.nb<-poly2nb(shape)#lista ngb per confinanti
is.list(pr.nb)#output è vero
(pr.nb[[233]])#neighbors regionale isolato
summary(pr.nb)
class(pr.nb)

plot(shape)
plot(pr.nb,coordinates(shape),add=TRUE,col="blue")
# Crea la matrice di pesi vicini
# Aggiungi collegamenti fittizi
pr.nb[[233]] <- sample(1:length(pr.nb), 1)

pr.listw<-nb2listw(pr.nb,style="W")#weighted ngb list
length(pr.listw);names(pr.listw)
pr.listw$weights[[1]]
pr.listw$weights[[233]]
pos<-c(1,pr.nb[[233]])
map4<-shape[pos,]
plot(map4)

colnames(shape@data)
shape@data$SIN_LUX<-rowSums(shape@data[,13:16],na.rm = F)
shape@data$SIN_POP<-rowSums(shape@data[,7:10],na.rm = F)
shape@data$POP<-log(shape@data$CityDens10)
pos<-which(!is.na(shape@data$PopExpo))#posizioni non mancanti
pos.ms<-which(is.na(shape@data$PopExpo))#posizioni mancanti
expo<-log(shape@data$PopExpo[pos]+1)
pop<-shape@data$POP[pos]
reg<-lm(expo~pop)
summary(reg)
pred<-as.vector(predict(reg,data.frame(pop=shape@data$POP[pos.ms])))
shape@data$PopExpo[pos.ms]<-sapply((exp(pred)-1),function(x)max(x,0))
pos<-which(!is.na(shape@data$LuxExpo))#posizioni non mancanti
pos.ms<-which(is.na(shape@data$LuxExpo))#posizioni mancanti
expo<-log(shape@data$LuxExpo[pos]+1)
pop<-shape@data$POP[pos]
reg<-lm(expo~pop)
summary(reg)
pred<-as.vector(predict(reg,data.frame(pop=shape@data$POP[pos.ms])))
shape@data$LuxExpo[pos.ms]<-sapply((exp(pred)-1),function(x)max(x,0))
shape<-NbMean(shape,"E_POP")
shape@data$E_POP<-shape@data$PopExpo*sum(shape@data$SIN_POP,na.rm = TRUE)/sum(shape@data$PopExpo)
shape@data$E_LUX<-shape@data$LuxExpo*sum(shape@data$SIN_LUX,na.rm = TRUE)/sum(shape@data$LuxExpo)
shape@data$struct<-rep(1:dim(shape@data)[1])
shape@data$unstruct<-rep(1:dim(shape@data)[1])
#nb=poly2nb(shape)
nb2INLA("nbgINLA.graph",pr.nb)
#nb2INLA("nbgINLA.graph",poly2nb(shape))
g<- inla.read.graph(filename = "nbgINLA.graph",size.only =F )
plot(g)
summary(g)
#inla(formula,family="gaussian",data=data.frame(),....)
##meshbuilder()
data<-shape@data
str(shape@data)
summary(pr.nb)
#install.packages("INLA",repos = c(getOption("repos"), INLA="https://inla.r-inla-download.org/R/testing"),dep=TRUE)
f.pop<-shape@data$SIN_POP~shape@data$HDIcity00+f(unstruct,model="iid")+f(struct,model="besag",graph = "g",scale.model = TRUE)

f.lux<-shape@data$SIN_LUX~shape@data$HDIcity00+f(unstruct,model="iid")+f(struct,model="besag",graph="g")
#f(name,model,...) generica sintassi della funzione R-INLA

m.pop<-inla(formula = f.pop,family="poisson",data=data,E=shape@data$E_POP,
            control.compute = list(dic=TRUE,cpo=TRUE),
            verbose = T,
            debug=T,safe = F,
            control.predictor = list(compute=TRUE,link=1))
summary(m.pop)
mod<- inla(f.pop,family="poisson",data=shape@data,E=E_POP,verbose = T)
nbgINLA.graph
sul+sp_shape.dbf
sul+sp_shape.shx
sul+sp_shape.shp

Elias T. Krainski

unread,
Jan 21, 2024, 1:00:44 AM1/21/24
to R-inla discussion group
Hi,
In the formulae where you have f(..., graph = "g") it should be f(..., graph = g). I'm attaching the code with a few changes as suggestions.

Regards,
Elias

--
You received this message because you are subscribed to the Google Groups "R-inla discussion group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/2cb9841b-2099-4934-a8e5-2d53927d727bn%40googlegroups.com.
code.R
Reply all
Reply to author
Forward
0 new messages