-
Notifications
You must be signed in to change notification settings - Fork 0
/
predictions and lower and upper limits of 95% for all data.R
146 lines (124 loc) · 4.27 KB
/
predictions and lower and upper limits of 95% for all data.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
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
library(readxl)
library(sp)
library(geoR)
library(ggplot2)
library(INLA)
setwd("C:/Users/ACER/Desktop/cupper data/1401-06-31/SamiroumGeochemistry-6352")
df <- read_xlsx("6352.xlsx",sheet = "Sheet1")
db <- data.frame(x=df$X_UTM,y=df$Y_UTM,cu=df$Cu)
##make data base
x1 <- c(min(db$x),min(db$y))
x2 <- c(min(db$x),max(db$y))
x3 <- c(max(db$x),max(db$y))
x4 <- c(max(db$x),min(db$y))
x5 <- c(min(db$x),min(db$y))
borders <- rbind(x1,x2,x3,x4,x5)
colnames(borders) <- c("x", "y")
coords <- matrix(c(db$x,db$y), ncol = 2)
colnames(coords) <- c("x", "y")
data <- db$cu
mydf <- structure(list(coords,data,borders),.Names = c("coords","data","borders"))
##show data with border
ggplot(data.frame(cbind(mydf$coords,mydf$data)))+
geom_point(aes(x, y, color =mydf$data), size = 2) +
coord_fixed(ratio = 1) +
scale_color_gradient(low = "blue", high = "orange") +
geom_path(data = data.frame(mydf$border), aes(x, y)) +
theme_bw()
coo <- mydf$coords
summary(dist(coo))
#max.edge is the largest allowed triangle length; the lower the number the higher the resolution
#offset = c(inside the boundary triangle, outside the boundary triangle)
max.edge = diff(range(mydf$coords))/500
bound.outer = diff(range(mydf$coords))/500
mesh <- inla.mesh.2d(loc=mydf$coords, max.edge = max.edge, offset=c(max.edge, bound.outer))
plot(mesh)
points(coo, col = "red")
# mesh2 <- inla.mesh.2d(loc = mydf$coords, max.edge=c(1000, 1000), offset = c(500,500),cutoff = 300)
# plot(mesh2)
# points(coo, col = "red")
## spde <- inla.spde2.pcmatern(mesh = mesh, alpha = 2, constr = TRUE)
spde <- inla.spde2.matern(mesh = mesh, alpha = 2, constr = TRUE)
indexs <- inla.spde.make.index("s", spde$n.spde)
A <- inla.spde.make.A(mesh = mesh, loc = coo)
# dimension of the projection matrix
dim(A)
# number of observations
nrow(coo)
# number of vertices of the triangulation
mesh$n
#see the elemens of each row sum to 1
rowSums(A)
#we construct a grid called coop
bb <- bbox(borders)
x <- seq(bb[1, "min"] - 1, bb[1, "max"] + 1, length.out = 100)
y <- seq(bb[2, "min"] - 1, bb[2, "max"] + 1, length.out = 100)
coop <- as.matrix(expand.grid(x, y))
ind <- point.in.polygon(coop[, 1], coop[, 2],borders[, 1], borders[, 2])
coop <- coop[which(ind == 1), ]
plot(coop, asp = 1)
Ap <- inla.spde.make.A(mesh = mesh, loc = coop)
dim(Ap)
# stack for estimation stk.e
stk.e <- inla.stack(
tag = "est",
data = list(y = db$cu),
A = list(1, A),
effects = list(data.frame(b0 = rep(1, nrow(coo))), s = indexs)
)
# stack for prediction stk.p
stk.p <- inla.stack(
tag = "pred",
data = list(y = NA),
A = list(1, Ap),
effects = list(data.frame(b0 = rep(1, nrow(coop))), s = indexs)
)
# stk.full has stk.e and stk.p
stk.full <- inla.stack(stk.e, stk.p)
formula <- y ~ 0 + b0 + f(s, model = spde)
res <- inla(formula,
data = inla.stack.data(stk.full),
control.predictor = list(
compute = TRUE,
A = inla.stack.A(stk.full)
)
)
index <- inla.stack.index(stk.full, tag = "pred")$data
pred_mean <- res$summary.fitted.values[index, "mean"]
pred_ll <- res$summary.fitted.values[index, "0.025quant"]
pred_ul <- res$summary.fitted.values[index, "0.975quant"]
dpm <- rbind(
data.frame(
east = coop[, 1], north = coop[, 2],
value = pred_mean, variable = "pred_mean"
),
data.frame(
east = coop[, 1], north = coop[, 2],
value = pred_ll, variable = "pred_ll"
),
data.frame(
east = coop[, 1], north = coop[, 2],
value = pred_ul, variable = "pred_ul"
)
)
dpm$variable <- as.factor(dpm$variable)
ggplot(dpm) + geom_tile(aes(east, north, fill = value)) +
facet_wrap(~variable, nrow = 1) +
coord_fixed(ratio = 1) +
scale_fill_gradient(
name = "Rainfall",
low = "blue", high = "orange"
) +
theme_bw()
## setwd("C:/Users/ACER/Desktop/cupper data/1401-06-31/SamiroumGeochemistry-6352/inla")
# save(pred_mean,file="pred_mean.Rda")
# save(pred_ll,file="pred_ll.Rda")
# save(pred_ul,file="pred_ul.Rda")
# data <- data.frame(east = coop[, 1], north = coop[, 2],value = pred_mean)
#
# ggplot(data) + geom_tile(aes(east, north, fill = value)) +
# scale_fill_gradient(
# name = "Rainfall",
# low = "blue", high = "orange"
# ) +
# theme_bw()