#####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)