Я использую geoR как часть моей диссертации MRes, чтобы исследовать пространственный подход к количественной оценке пространственного распределения. Я столкнулся с проблемой, и я не уверен, как двигаться дальше.
По сути, мой набор данных поступает с нескольких платформ наблюдений, и мне нужно учитывать различия между ними при оценке параметров ковариационной функции Матерна. Я использую симуляцию Гаусса-Маркова со случайным полем (GMRF) для тестирования трех различных сценариев ios:
1 Нет эффекта сосуда (параметры Матера оцениваются с использованием данных, полученных непосредственно из GMRF)
2 Эффект судна игнорируется (образцы делятся между различными платформами наблюдения и перед оценкой получают разные коэффициенты отбора проб)
3 Эффект судна включен (используется тот же набор данных, что и в сценарии 2, но судно используется в качестве сопутствующего варьировать, используя trend <- trend.spatial (h.dat2, trend = "2nd", add.to.trend = ~ Vessel) </p>
Конечная цель моего подхода заключается в том, чтобы сценарии 1 и 3 имели одинаковое значение, если не очень похожие оценки параметров. Мне удалось сделать это с помощью likfit, где:
(likfit3 <- likfit(h.dat2,
cov.model="matern",
lik.method = "ML",
trend = trend,
kappa = ini.kappa,
fix.kappa = FALSE,
ini.cov.pars = c(0.5, min(dist(h.dat2$coords))), # Sigmas, phi
fix.nugget = FALSE))
Следующий шаг, чтобы следовать тому же подходу с использованием кригинга для создания поверхности распределения. У меня возникли проблемы с trend.d и trend.l в krige.conv (). Я понимаю, что trend.d относится к данным и может быть тем же трендом, что и в likfit ()valuation.trend.l относится к сетке местоположений, подлежащих оценке. Однако я сделал это вместо ошибки:
Ошибка в файле krige.conv (h.dat2, location = pred.grid, krige = krige.control (trend.d = trend,: location и trend.l имеют несовместимые размеры
Теперь я получаю: Ошибка в trend.d! = trend.l: несовместимые массивы
Ближайшее решение, которое я нашел, находится здесь: http://r-sig-geo.2731867.n2.nabble.com/geoR-krigeing-predict-on-other-coordinates-with-a-parametric-trend-td6792393.html
Однако пример решения я не смог найти.
Воспроизводимый код выглядит следующим образом. Интересующий фрагмент кода начинается с ### ОЦЕНКИ БЕТАСОВ БИАСА
# Packages
if (!require('RandomFields')) install.packages('RandomFields'); library('RandomFields')
if (!require('ggplot2')) install.packages('ggplot2'); library('ggplot2')
if (!require('raster')) install.packages('raster'); library('raster')
if (!require('viridis')) install.packages('viridis'); library('viridis')
if (!require('geoR')) install.packages('geoR'); library('geoR')
if (!require('dplyr')) install.packages('dplyr'); library('dplyr')
if (!require('gridExtra')) install.packages('gridExtra'); library('gridExtra')
# FUNCTIONS
# Matern covaraince
Matern.f <- function(kappa,phi,distance)
# kappa = smoothness or nu parameter
# phi = range of scale parameter
# distance = vector of distances
{((2^(1-kappa)) / (gamma(kappa))) * (sqrt(2*kappa) * (distance/phi))^kappa * besselK(sqrt(kappa*2)*(distance/phi), kappa)}
SCALE <- 10# 1/SCALE. This is the scale in comparision to the Celtic Sea in km
n <- round(3302/SCALE,0)# NUMBER OF SAMPLES
i <- 15 # seed value
# RANDOM FIELD SUMLATION PARAMATERS
p <- True_phi <- 7#/SCALE # phi = range parameter or scale in random fields
k <- True_kappa <- 0.40 # kappa = Smoothness or nu in random fields
s <- True_sigmasq <- 3 # Varaince of the simulation
True_mu <- 100 # mean
# d <- seq(0.01, sqrt(366557)/SCALE, length = 1e3) # Vector of distances
d <- seq(0.01, 300, length = 1e3) # Vector of distances
seed <- i; set.seed(seed) # Set seed value
res_x <- res_y <- seq(0, sqrt(366557)/SCALE , 0.3) # Resolution of simulation
# Prediction grid for interpolation
pred.grid = expand.grid(res_x, res_y)
ini.kappa <- k*.8 # initial kappa starting value
# PLOTTING DEFAULTS
xlab <- "Distance (km)"
ylab <- "Covaraince"
lwd <- 2 # line width in all plots
xlim <- c(0,4 * p) # max distance
ylim <- c(0,1)
# Defining simulation model
model <- RMmatern(nu = k, scale = p, var = s) + RMtrend(mean = True_mu ) # for RandomFeilds simulation
# nu = Smoothness, also known as kappa
# scale = Range , also known as phi
# RMtrend = simulation mean
# Plotting the matern function according to the input parameters
Simulations_pars <- Matern.f(kappa = k, phi = p, d = d)
# Table to populate with results of each estimation
ResultsTable <- t(data.frame( Truth = c(p, k, NA)))
colnames(ResultsTable) <- c("phi", "kappa", "sigmasq")
# Generate Simulation
simu <- RFsimulate(model, res_x, res_y, grid=TRUE, seed = i, RFoptions(spConform=FALSE)) # Generate simulation
############################################# NOTE LOG TRANSFORMATION HAPPENS HERE #############################################
simu$variable1 <- log(simu$variable1)
############################################# NOTE LOG TRANSFORMATION HAPPENS HERE #############################################
simu.rast <- raster(simu); r.spdf <- as(simu, "SpatialPixelsDataFrame") # Convert to Raster
r.df <- as.data.frame(r.spdf); colnames(r.df) <- c("z", "x", "y") # Convert to dataframe
# Sample from the GMRF at each XY location
{xy <- coordinates(simu)
pts <- base::sample(nrow(xy), min(n, nrow(xy) / 2))
data <- matrix(nrow=nrow(xy), as.vector(simu))[pts, ]
data <- cbind(xy, simu@data)[pts,]
data$Point <- 1:n
colnames(data) <- c("x","y","z", "Point")
rownames(data) <- 1:n
data$Type <- factor("Sample")
# plot(simu, data, asp = 1)
MaxDist <- max(dist(cbind(data$x, data$y))) # Used as a scaling value
# Populate xyz dataframe
x <- data$x
y <- data$y
z <- data$z
# Sample data set
xyz <- data.frame(x=x,y=y,z=z)}
```
## View the true Mátern covariance strucutre generated using the true parameters
```{r}
# Plot the curve of truth
ggplot(data = data.frame(x =Simulations_pars, y = d), aes(x = y, y = x))+ geom_line()+ theme(aspect.ratio=1) + xlim(xlim) + xlab(xlab) + ylab(ylab) + ylim(ylim)+ ggtitle("Simulation Matern covaraince functon") + labs(subtitle =paste0("Simulations pars ", "p=", p, " k=", k))
```
## View sample locations from simulation
```{r}
# Plot of simulation
Sim <- ggplot() + geom_tile(data=r.df, aes(x=x, y=y, fill=z)) + scale_fill_viridis() + coord_equal() + labs( x = "x km", y = "y km" ) + labs(title = paste("0 Original Simulation"), subtitle = paste("phi =", p, " kappa = ", k, " mu = ", #log(True_mu)
True_mu ))+ geom_point(data = xyz, aes(x, y, size = z),shape = 1, inherit.aes= FALSE) + guides(size = guide_legend(reverse=TRUE))
Sim
```
# Scenatio 1 - Estimation without vessel effect
```{r pressure, echo=FALSE}
h.dat1 <- as.geodata(xyz, coords.col = c(1,2), data.col = 3)
(likfit1 <- likfit(h.dat1,
cov.model="matern",
lik.method = "ML",
trend = "2nd",
kappa = ini.kappa,
fix.kappa = FALSE,
ini.cov.pars = c(0.1, min(dist(h.dat1$coords))), # Sigmas, phi
fix.nugget = FALSE))
Table1 <- rbind(ResultsTable, t(data.frame(No_vessel_effect = c(likfit1$phi, likfit1$kappa, likfit1$sigmasq)) ))
round(Table1,2)
# Krige without vessel effect
kcE1 = krige.conv(h.dat1, loc = pred.grid, krige = krige.control(trend.d = "2nd", trend.l = "2nd", obj.m = likfit1))
# Formating for ggplot
r1 <- SpatialPointsDataFrame(pred.grid, data.frame(predict = kcE1$predict))
gridded(r1) <- TRUE
r1 <- as(r1,'RasterLayer')
test_spdf1 <- as(r1, "SpatialPixelsDataFrame")
test_df1 <- as.data.frame(test_spdf1)
colnames(test_df1) <- c("value", "x", "y")
Final1 <- data.frame(long = test_df1$x, lat = test_df1$y, z = test_df1$value)
#ggplot
Krige1 <- ggplot(Final1, aes(x=long, y=lat, fill=z)) + geom_raster() +
scale_fill_viridis_c(option = "viridis") +
theme(aspect.ratio=1) + labs(title = paste("1 No Vessel Effect"), x = "x km", y = "y km", subtitle = paste("phi =", round(likfit1$phi,2), " kappa = ", round(likfit1$kappa,2) )) + geom_point(data = xyz, aes(x, y, size = z),shape = 1, inherit.aes= FALSE) + guides(size = guide_legend(reverse=TRUE))
```
## Introduce a vessel effect to dataset
```{r}
# INTRODUCING BIAS
# Bias application
Bias1 <- function(z)
# applies simple bias to the z value of a dataset
{bias <- z * BIAS1
return(bias)}
Bias2 <- function(z)
# applies simple bias to the z value of a dataset
{bias <- z * BIAS2
return(bias)}
# BIAS WITHIN VESSEL EFFECT
## Number of samples per vessel similar to CS
A <- round(n/100*57,0)
B <- round(n/100*24,0)
C <- round(n/100*19,0)
A+B+C==n
BIAS1 <- .9 # bias applied to vesselB. Can be a positive or negative number
BIAS2 <- .95
# E1 <- 0--BIAS1 # for labelling
# E2 <- 0--BIAS2 # for labelling
# SPLITTING ORIGINAL SAMPLES INTO TWO RANDOM PARTS
# VESSELA
VesselA <- dplyr::sample_n(xyz, A)
colnames(VesselA) <- c("x","y","z")
VesselA$Vessel <- "A"
# VESSEL B
VesselB.tmp <- dplyr::anti_join(xyz, VesselA, by = c("x", "y", "z"))
colnames(VesselB.tmp) <- c("x","y","z")
VesselB.tmp$Vessel <- "B"
# Vessel C
VesselC <- dplyr::sample_n(VesselB.tmp, C)
VesselC$Vessel <- "C"
VesselC.unbias <- VesselC
# Vessel B clean up
VesselB.unbias <- VesselB <- dplyr::anti_join(VesselB.tmp, VesselC, by = c("x", "y", "z"))
nrow(VesselA) == A
nrow(VesselB) == B
nrow(VesselC) == C
VesselB$z <- Bias1(VesselB$z) # INTRODUCED BIAS
VesselC$z <- Bias2(VesselC$z)
# Combine into one dataset
Data.broken <- rbind(VesselA, VesselB, VesselC)
#Data.broken$log_z <- Data.broken$z
Data.broken <- cbind(Data.broken, p, k, True_mu, seed) # data set for export if needed
# Data.broken$x <- Data.broken$x/MaxDist
# Data.broken$y <- Data.broken$y/MaxDist
```
## Scenaro 3 - Estimation with Vessel Effect Included
```{r}
### ESTIMATING BETAS OF BIAS
trend <- trend.spatial(h.dat2, trend = "2nd", add.to.trend = ~ Vessel)
(likfit3 <- likfit(h.dat2,
cov.model="matern",
lik.method = "ML",
trend = trend,
kappa = ini.kappa,
fix.kappa = FALSE,
ini.cov.pars = c(0.5, min(dist(h.dat2$coords))), # Sigmas, phi
fix.nugget = FALSE))
# trend.d with predictions set to match vessel A
trend.d <- trend
trend.d[,7:8] <- 0
# trend.l using prediction grid locations
Loc <- pred.grid
names(Loc) <- c("x","y")
Loc$z <- NA
Loc$Vessel <- NA
total <- merge(Loc, Data.broken[,1:4], by=c("x","y"), all = TRUE)
total$z.x <- NULL
total$Vessel.x <- NULL
colnames(total) <- c("x", "y", "z", "Vessel")
for (i in 1:nrow(total)){
print(i)
ifelse( is.na(total$z[i]),
{total$z[i] <- 0
total$Vessel[i] <- "D"}
,next)}
# Make geodata
d.dat2 <- as.geodata(total, coords.col = c(1,2), data.col = 3, covar.col = 4)
# Make trend.l
trend.l <- trend.spatial(d.dat2, trend = "2nd", add = ~ Vessel)
trend.l[,7:9] <- 0
# Krige with vessel effect included
kcE3 = krige.conv(h.dat2, # Data
locations = pred.grid,# prediction grid
krige = krige.control(
trend.d = trend.d, # covaraites at data locations
trend.l = trend.l[,1:8], # covaraites at predicted locations
type.krige = "OK",
obj.m = likfit3)) # Model parameters
# Formating for ggplot
r3 <- SpatialPointsDataFrame(pred.grid, data.frame(predict = exp(kcE3$predict)))
gridded(r3) <- TRUE
r3 <- as(r3,'RasterLayer')
test_spdf3 <- as(r3, "SpatialPixelsDataFrame")
test_df3 <- as.data.frame(test_spdf3)
colnames(test_df3) <- c("value", "x", "y")
Final3 <- data.frame(long = test_df3$x, lat = test_df3$y, z = test_df3$value)
# ggplot
Krige3 <- ggplot(Final3, aes(x=long, y=lat, fill=z)) + geom_raster() +
scale_fill_viridis_c(option = "viridis") +
theme(aspect.ratio=1) + labs(title = paste("3 Vessel Effect Included"), x = "x km", y = "y km", subtitle = paste("phi =", round(likfit3$phi,2), " kappa = ", round(likfit3$kappa,2) ))+ geom_point(data = xyz, aes(x, y, size = z),shape = 1, inherit.aes= FALSE) + guides(size = guide_legend(reverse=TRUE))
```
## Compare Scenario 1 No Vessel Effect and Scenario 3 Vessel Effect Accounted For
```{r}
grid.arrange(Krige1, Krige3)
```