require(sp) require(spdep) require(INLA) require(raster) require(ggplot2) require(SpatialEpi) require(maptools) require(spatialEco) rm(list = ls()) cpn<-read.csv("cpn.csv") full<-read.csv("grid.csv") mpios<-readShapePoly("mpios.shp") ful1<-merge(cpn, full, by.x="id", by.y="NEAR_FID", all.y=T) coordinates(ful1) = ~ X_DD+Y_DD mpios$label<-seq(1, nrow(mpios)) a.data <- over(ful1, mpios[,c("NAME_0", "NAME_1", "label")]) ful1$Country<-NULL ful1$ColPro2<-NULL ful1$label<-as.numeric(as.character(a.data$label)) ful1$country<-as.numeric(as.factor(a.data$NAME_0)) ful1$depto<-as.numeric(as.factor(a.data$NAME_1)) ful1$CP<-as.numeric(as.factor(ful1$ColPro)) rast <- raster(ncol = 482, nrow = 658) extent(rast) <- extent(ful1) cras<-rasterize(ful1, rast, "coca", fun = sum) nras<-rasterize(ful1, rast, "NEAR_DIST", fun = median) pras<-rasterize(ful1, rast, "CP", fun = mean) oras<-rasterize(ful1, rast, "country", fun = mean) dras<-rasterize(ful1, rast, "depto", fun = mean) lras<-rasterize(ful1, rast, "label", fun = mean) cras1<-aggregate(cras, 3, max) nras1<-aggregate(nras, 3, median) pras1<-aggregate(pras, 3, min) oras1<-aggregate(oras, 3, min) dras1<-aggregate(dras, 3, min) lras1<-aggregate(lras, 3, min) coca<-as.data.frame(cras1, xy=T) near<-as.data.frame(nras1, xy=T) proj<-as.data.frame(pras1, xy=T) coun<-as.data.frame(oras1, xy=T) dept<-as.data.frame(dras1, xy=T) labe<-as.data.frame(lras1, xy=T) rm(nras1, pras1, oras1, dras1, lras1) reco<-coca reco$near<-near$layer reco$project<-as.factor(as.integer(proj$layer)) reco$country<-as.factor(as.integer(coun$layer)) reco$depto<-as.factor(as.integer(dept$layer)) reco$label<-as.integer(labe$layer) colnames(reco)[3]<-"coca" reco<-na.omit(reco) c.data <- over(mpios, ful1) coordinates(c.data)<-coordinates(mpios) nb.map<-dnearneigh(c.data, d1 = 0, d2 = 2) nb2INLA("map.graph",nb.map) reco$dum<-rep(1, dim(reco)[1]) reco$d2p<-log10((reco$near/1000)+1) rm(coca, near, proj, coun, dept, labe) g = inla.read.graph("map.graph") f0 = coca~f(label, model="besag", graph = g, param=c(1,0.01)) mod0 = inla(f0, data=reco, family="binomial", Ntrials=dum, control.inla = list(strategy = "gaussian", int.strategy="ccd"), control.predictor=list(compute=TRUE), control.compute=list(waic=TRUE)) h.mod0 = inla.hyperpar(mod0) f1 = coca~f(label, model="besag", graph = g, param=c(1,0.01))+d2p mod1 = inla(f1, data=reco, family="binomial", Ntrials=dum, control.inla = list(strategy = "gaussian", int.strategy="ccd"), control.predictor=list(compute=TRUE), control.compute=list(waic=TRUE)) h.mod1 = inla.hyperpar(mod1) f2 = coca~f(label ,model="besag", graph = g, param=c(1,0.01))+f(project, model="iid", param=c(0.001,0.001))+d2p mod2 = inla(f2, data=reco, family="binomial", Ntrials=dum, control.inla = list(strategy = "gaussian", int.strategy="ccd"), control.predictor=list(compute=TRUE), control.compute=list(waic=TRUE)) h.mod2 = inla.hyperpar(mod2) f3 = coca~f(label, model="besag", graph = g, param=c(1,0.01))+f(country, model="iid", param=c(0.001,0.001))+d2p mod3 = inla(f3, data=reco, family="binomial", Ntrials=dum, control.inla = list(strategy = "gaussian", int.strategy="ccd"), control.predictor=list(compute=TRUE), control.compute=list(waic=TRUE)) #h.mod3 = inla.hyperpar(mod3) f4 = coca~f(label, model="besag", graph = g, param=c(1,0.01))+f(depto, model="iid", param=c(0.001,0.001))+d2p mod4 = inla(f4, data=reco, family="binomial", Ntrials=dum, control.inla = list(strategy = "gaussian", int.strategy="ccd"), control.predictor=list(compute=TRUE), control.compute=list(waic=TRUE)) h.mod4 = inla.hyperpar(mod4) f5 = coca~f(label, model="bym", graph = g, param=c(1,0.01)) mod5 = inla(f5, data=reco, family="binomial", Ntrials=dum, control.inla = list(strategy = "gaussian", int.strategy="ccd"), control.predictor=list(compute=TRUE), control.compute=list(waic=TRUE)) #h.mod5 = inla.hyperpar(mod5) f6 = coca~f(label, model="bym", graph = g, param=c(1,0.01))+d2p mod6 = inla(f6, data=reco, family="binomial", Ntrials=dum, control.inla = list(strategy = "gaussian", int.strategy="ccd"), control.predictor=list(compute=TRUE), control.compute=list(waic=TRUE)) #h.mod6 = inla.hyperpar(mod6) sink("spatial_models3_hires.txt") print(summary(h.mod0)) print(summary(h.mod1)) print(summary(h.mod2)) print(summary(mod3)) print(summary(h.mod4)) print(summary(mod5)) print(summary(mod6)) sink() save.image("d_models3.Rdata")