-
Notifications
You must be signed in to change notification settings - Fork 2
/
Run models unstructured bias covariate.R
123 lines (85 loc) · 4.43 KB
/
Run models unstructured bias covariate.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
## Function to run unstructured only models with simulated data
unstructured_model_cov <- function(unstructured_data, dat1, biasfield, dim = dim, plotting=FALSE, biascov){
#packages
library(INLA)
library(reshape2)
library(deldir)
library(rgeos)
library(fields)
#preparation - mesh construction - use the loc.domain argument
mesh <- inla.mesh.2d(loc.domain = biasfield[,c(1,2)],max.edge=c(20,40),cutoff=2, offset = c(5,20))
#plot the mesh to see what it looks like
if(plotting == TRUE){plot(mesh)}
##set the spde representation to be the mesh just created
spde <- inla.spde2.matern(mesh)
#make A matrix for unstructured data
unstructured_data_A <- inla.spde.make.A(mesh = mesh, loc = as.matrix(unstructured_data[,1:2]))
#make integration stack for unstructured data
#get dimensions
max_x <- max(biasfield$x)
max_y <- max(biasfield$y)
loc.d <- t(matrix(c(0,0,max_x,0,max_x,max_y,0,max_y,0,0), 2))
#make dual mesh
dd <- deldir::deldir(mesh$loc[, 1], mesh$loc[, 2])
tiles <- deldir::tile.list(dd)
#make domain into spatial polygon
domainSP <- SpatialPolygons(list(Polygons(list(Polygon(loc.d)), '0')))
#intersection between domain and dual mesh
poly.gpc <- as(domainSP@polygons[[1]]@Polygons[[1]]@coords, "gpc.poly")
# w now contains area of voronoi polygons
w <- sapply(tiles, function(p) rgeos::area.poly(rgeos::intersect(as(cbind(p$x,p$y), "gpc.poly"), poly.gpc)))
#check some have 0 weight
table(w>0)
nv <- mesh$n
n <- nrow(unstructured_data)
#change data to include 0s for nodes and 1s for presences
y.pp <- rep(0:1, c(nv, n))
#add expectation vector (area for integration points/nodes and 0 for presences)
e.pp <- c(w, rep(0, n))
#diagonal matrix for integration point A matrix
imat <- Diagonal(nv, rep(1, nv))
#combine integration point A matrix with data A matrix
A.pp <- rbind(imat, unstructured_data_A)
#get covariate for integration points
covariate = dat1$gridcov[Reduce('cbind', nearest.pixel(
mesh$loc[,1], mesh$loc[,2],
im(dat1$gridcov)))]
biascovariate = biascov[Reduce('cbind', nearest.pixel(mesh$loc[,1], mesh$loc[,2], im(biascov)))]
# Create data stack
stk_unstructured_data <- inla.stack(data=list(y=y.pp, e = e.pp),
effects=list(list(data.frame(interceptB=rep(1,nv+n), env = c(covariate, unstructured_data$env), bias = c(biascovariate, unstructured_data$bias))), list(Bnodes=1:spde$n.spde)),
A=list(1,A.pp),
tag="unstructured_data")
source("Create prediction stack.R")
join.stack <- create_prediction_stack(stk_unstructured_data, c(10,10), biasfield = biasfield, dat1 = dat1, mesh, spde)
formulaN = y ~ -1 + interceptB + env + bias + f(Bnodes, model = spde)
result <- inla(formulaN,family="poisson",
data=inla.stack.data(join.stack),
control.predictor=list(A=inla.stack.A(join.stack), compute=TRUE),
control.family = list(link = "log"),
E = inla.stack.data(join.stack)$e,
control.compute = list(cpo=TRUE, waic = TRUE, dic = TRUE)
)
result$summary.fixed
##project the mesh onto the initial simulated grid
proj1<- inla.mesh.projector(mesh,ylim=c(1,max_y),xlim=c(1,max_x),dims=c(max_x,max_y))
##pull out the mean of the random field
xmean1 <- inla.mesh.project(proj1, result$summary.random$Bnodes$mean)
##plot the estimated random field
# plot with the original
# some of the commands below were giving warnings as not graphical parameters - I have fixed what I can
# scales and col.region did nothing on my version
if(plotting == TRUE){png("unstructured_model with bias covariate.png", height = 1000, width = 2500, pointsize = 30)
par(mfrow=c(1,3))
image.plot(1:dim[1],1:dim[2],xmean1, col=tim.colors(),xlab='', ylab='',main="mean of r.f",asp=1)
#plot truth
image.plot(list(x=dat1$Lam$xcol*100, y=dat1$Lam$yrow*100, z=t(dat1$rf.s)), main='Truth', asp=1) # make sure scale = same
points(unstructured_data[,1:2], pch=16)
##plot the standard deviation of random field
xsd1 <- inla.mesh.project(proj1, result$summary.random$Bnodes$sd)
#library(fields)
image.plot(1:dim[1],1:dim[2],xsd1, col=tim.colors(),xlab='', ylab='', main="sd of r.f",asp=1)
dev.off()}
#return from function
return(list(join.stack = join.stack, result = result))
}