geoR krige.conv () Полином 2-го порядка, тренд.d и тренд. - PullRequest
0 голосов
/ 01 марта 2020

Я использую 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)

```
...