In nmmarquez/PointPolygon: Evaluate Point and Polygon Data Simulatneously

```{R, warning=FALSE, message=F, echo=F} rm(list=ls()) library(arm) library(rgeos) library(sp) require(tidyr) require(gridExtra) require(dplyr) require(ggplot2) library(PointPolygon) library(stringr)

range of the underlying spatial process as defined

rangeE <- .5 # range of spatial prces varies from {.3, .5, .7} covVal <- 2 # covariate effect in set {.2, .4, -.5, .2, -2} M <- 100 # number of samples in Polygons chosen from U(50, 300) seed <- 12345 # RNG

set.seed(seed)

$$
i,j \in \{1, \dots, n\}\\
Y_i \sim \text{Binomial}(N_i, p_i) \\
\text{logit}(p_i) = X_i \boldsymbol{\beta} + \eta(s_i)  \\
\boldsymbol{\eta} \sim \text{MVN}(\boldsymbol{0}, \boldsymbol{\Sigma}) \\
\Sigma_{ij} = \frac{\sigma^2_\eta}{2^{\nu-1} \Gamma(\nu)}
    (\kappa ||s_i - s_j||)^\nu K_\nu(\kappa ||s_i - s_j||) \\
\Sigma_{ij} = \text{Cov}(\eta(s_i) , \eta(s_j))
$$
$$
Y_k \sim \text{Binomial}(N_k, \hat{p}_k) \\
\text{logit}(\hat{p}_k) =  \boldsymbol{\beta}\Bigg(\int_{\forall l \in \mathcal{A}_k}X_l \Bigg) + \Bigg(\int_{\forall l \in \mathcal{A}_k}\eta(s_l) \Bigg) \\
$$

$$
Y_k \sim \text{Binomial}(N_k, \hat{p}_k) \\
\hat{p}_k =  \int_{\forall l \in \mathcal{A}_k} \text{inv.logit}(X_l \boldsymbol{\beta} + \eta(s_l)) \\
$$


```{R, warning=FALSE, message=F, echo=F}
set.seed(seed)
unitSim <- simField(
    N = 60, # how detailed the grid is creates an NxN grid
    offset = c(0.1, 0.2), # all simulations will have this to create pred mesh 
    max.edge = c(0.05,0.2), # all simulations will have this to create pred mesh
    beta0 = -2, # the intercept term
    rangeE = rangeE,
    betaList = list(list(type="random", value=covVal))) # the cov type & value


# just some R plotting 
plotList <- lapply(c("V1", "z", "theta"), function(eff){
    unitSim$spdf %>%
        gather("Effect", "Value", V0:theta) %>%
        filter(Effect==eff) %>%
        ggplot(aes(x, y, fill=Value)) +
        geom_raster() +
        coord_equal() +
        theme_void() +
        scale_fill_distiller(palette = "Spectral") +
        ggtitle(NULL) +
        guides(fill=FALSE)
})

p1 <- "/tmp/1.png"
p2 <- "/tmp/2.png"
p3 <- "/tmp/3.png"

ggsave(p1, plotList[[1]], "png", width=.4, height=.4)
ggsave(p2, plotList[[2]], "png", width=.4, height=.4)
ggsave(p3, plotList[[3]], "png", width=.4, height=.4)

General Functional Form For Probability Surface

$$ \boldsymbol{p} = \text{inv.logit}(\beta_0 + \beta_1 \boldsymbol{X} + \eta(\boldsymbol{s})) $$

Simulation 1: Randomly Distributed Covariate Effect

$=\text{inv.logit}(\beta_0 + \beta_1 \times$

$+$

$)$

```{R, warning=FALSE, message=F, echo=F} set.seed(seed) unitSim <- simField( N = 60, # how detailed the grid is creates an NxN grid offset = c(0.1, 0.2), # all simulations will have this to create pred mesh max.edge = c(0.05,0.2), # all simulations will have this to create pred mesh beta0 = -2, # the intercept term rangeE = rangeE, betaList = list(list(type="spatial", value=covVal))) # the cov type & value

just some R plotting

plotList <- lapply(c("V1", "z", "theta"), function(eff){ unitSim$spdf %>% gather("Effect", "Value", V0:theta) %>% filter(Effect==eff) %>% ggplot(aes(x, y, fill=Value)) + geom_raster() + coord_equal() + theme_void() + scale_fill_distiller(palette = "Spectral") + ggtitle(NULL) + guides(fill=FALSE) })

p1 <- "/tmp/4.png" p2 <- "/tmp/5.png" p3 <- "/tmp/6.png"

ggsave(p1, plotList[[1]], "png", width=.4, height=.4) ggsave(p2, plotList[[2]], "png", width=.4, height=.4) ggsave(p3, plotList[[3]], "png", width=.4, height=.4)

#### Simulation 2: Spatially Correlated Covariate Effect

<center>
![](/tmp/6.png) $=\text{inv.logit}(\beta_0 + \beta_1 \times$ ![](/tmp/4.png) $+$ ![](/tmp/5.png) $)$
</center>

```{R, warning=FALSE, message=F, echo=F}
set.seed(seed)
unitSim <- simField(
    N = 60, # how detailed the grid is creates an NxN grid
    offset = c(0.1, 0.2), # all simulations will have this to create pred mesh 
    max.edge = c(0.05,0.2), # all simulations will have this to create pred mesh
    beta0 = -2, # the intercept term
    rangeE = rangeE,
    betaList = list(list(type="cluster", value=covVal))) # the cov type & value


plotList <- lapply(c("V1", "z", "theta"), function(eff){
    unitSim$spdf %>%
        gather("Effect", "Value", V0:theta) %>%
        filter(Effect==eff) %>%
        ggplot(aes(x, y, fill=Value)) +
        geom_raster() +
        coord_equal() +
        theme_void() +
        scale_fill_distiller(palette = "Spectral") +
        ggtitle(NULL) +
        guides(fill=FALSE)
})

p1 <- "/tmp/7.png"
p2 <- "/tmp/8.png"
p3 <- "/tmp/9.png"

ggsave(p1, plotList[[1]], "png", width=.4, height=.4)
ggsave(p2, plotList[[2]], "png", width=.4, height=.4)
ggsave(p3, plotList[[3]], "png", width=.4, height=.4)

Simulation 3: Clustered Covariate Effect

$=\text{inv.logit}(\beta_0 + \beta_1 \times$

$+$

$)$

pointSampleDF <- samplePoints(unitSim, N=250, M=50)
polyFrame <- dividePolygon(unitSim$bound, 3)
polyFrame$id <- 1:9
polyDF <- fortify(polyFrame, region="id") %>%
    mutate(id=as.numeric(id)) %>%
    left_join(polyFrame@data, by="id")
boundDF <- fortify(dividePolygon(unitSim$bound, 1))

ggplot() +
    geom_path(
        aes(long,lat, group=group, fill=NULL, z=NULL),
        color="black",
        size=.8,
        data=boundDF) +
    theme_void() + 
    geom_raster(
        aes(x=x, y=y, fill=is.na(obs)),
        data=pointSampleDF %>% right_join(unitSim$spdf)) +
    scale_fill_manual(
        breaks = c(TRUE, FALSE),
        values=c("black", "white")) +
    guides(fill=FALSE)

polyDF %>%
    mutate(val="1") %>%
    ggplot(aes(long,lat,group=group, fill=val)) + 
    geom_polygon() +
    geom_path(color="black", size=.1) +
    scale_fill_manual(breaks="1", values="grey50") +
    theme_void() +
    guides(fill=FALSE)

ggplot() +
    theme_void() +
    geom_raster(
        aes(x=x, y=y, fill=is.na(obs)),
        data=pointSampleDF %>%
            right_join(unitSim$spdf) %>%
            mutate(val=as.character(is.na(obs) + 2))) +
    scale_fill_manual(
        breaks = c("1", "2", "3"),
        values=c("black", "white", "grey50")) +
    geom_path(
        aes(long,lat, group=group, fill=NULL, z=NULL),
        color="black",
        size=.3,
        data=polyDF) +
    # geom_polygon(
    #     aes(long,lat,group=group,fill="grey50"),
    #     data = polyDF %>%
    #         mutate(val=factor("3", as.character(1:3)))) +
    guides(fill=FALSE)

fits <- readRDS("./range=0.5,cov=2,covtype=cluster,M=200,seed=3.Rds")
sampling <- "rwidth_3 poly"
predList <- list(
            riemann = fits$pred$riemann[[sampling]],
            utazi = fits$pred$utazi[[sampling]],
            point = fits$pred$point$point
        )
field <- fits$sim
DF <- field$spdf@data
allDF <- DF
allDF$Type <- "True"
allDF$sd <- NA
for(i in 1:length(predList)){
    DFpred <- DF
    DFpred$Type <- names(predList)[i]
    DFpred$theta <- predList[[i]]$mu
    DFpred$sd <- predList[[i]]$sd
    allDF <- dplyr::bind_rows(allDF, DFpred)
}

new <- field
new$spdf@data <- allDF %>%
    mutate(Type=str_to_title(Type))

ggplot(new$spdf@data, aes_string("x", "y", fill="theta")) +
    geom_raster(alpha=.7) +
    coord_equal() +
    theme_void() +
    scale_fill_distiller(palette = "Spectral") +
    facet_wrap(~Type) +
    labs(fill="Probability")


new$spdf@data <- allDF[allDF$Type != "True",]
ggplot(new$spdf@data, aes_string("x", "y", fill="sd")) +
    geom_raster(alpha=.7) +
    coord_equal() +
    theme_void() +
    scale_fill_distiller(palette = "Spectral") +
    facet_wrap(~Type) +
    labs(fill="SD")