Presentations will be 10 minutes plus 5 minutes questions.
Course schedule
No class on 5th and 26th November
Week 1: Introduction to Spatial Data Science with R
Week 2-3: Areal data
Week 4: Dashboards
Week 5: Geostatistical data
Week 6: Student presentations project 1
Week 7: Geostatistical data
Weeek 8-10: Point patterns
Week 11: Student presentations project 2
Week 12-13: Point patterns
Week 14: Interactive web applications
Week 15: Student presentations project 3
Evaluation
General guidelines
The evaluation of STAT 300 consists of 3 projects to be done individually:
Project 1: Areal data analysis (33%). Due date October 4 by 5pm
Project 2: Geostatistical data analysis (33%). Due date November 8 by 5pm
Project 3: Point pattern data analysis (34%). Due date December 6 by 5pm
For each of the projects, students will look for an appropriate dataset and analyze it using the statistical techniques learned during the course. They will need to submit a written report and reproducible R code, and give an oral presentation about their work.
Students will need to submit their work via Blackboard by the due date. Late coursework will not be accepted, unless prior arrangements have been made.
Files to be submitted
Students will need to submit via Blackboard a zip file called LastnameFirstname.zip containing the following:
written report in PDF format (LastnameFirstname.pdf),
folder data with the data,
folder code with reproducible R code (code.R). The code should load data that is in the data folder, work without errors, and include comments,
presentation to be shown in the oral presentation (PDF, HTML or PowerPoint) (LastnameFirstnamePresentation.pdf).
Written report
The written report will include a title, abstract, introduction, data, methods, results and interpretation, a brief discussion section, and references. Make sure you wholly and accurately interpret your results, tables and figures in the context of the study. The report does not have to include R code. The R code will be submitted separately. Tables and figures captions need to be self-contained. For the report, please use a maximum of 10 pages with 12 point font, one column, and 2cm margins, including tables and figures. Reports that exceed the specified number of pages will not be assessed. You will be judged on the quality and not on the quantity of material you submit.
Oral presentation
In the oral presentation you can use slides, animations and/or dashboards to show your work. You should show good communication and presentation skills. Presentations will be 10 minutes plus 5 minutes questions.
Data
The data can be from your own research (and that you have permission to openly share) or data that you can find in open repositories. The website https://rspatialdata.github.io/ can be helpful to find data. Examples of data that can be used for each of the projects are the following.
An example of geostatistical data would be air quality data in the US measured at specific monitoring stations. These data can be found at https://www.epa.gov/outdoor-air-quality-data.
Project 3: Point pattern data
Point pattern datasets can be found in the spatstat package. For example, data(betacells) contains data of cells in the retina.
Evaluation
The project will be marked out of 100. There is no single correct analysis for this project, so you will not be scored based on how close you get to one particular answer. Instead, the marks will be allocated according to the following principles.
70-100 A project that could be presented with little or no revision. Analysis should be soundly done so that conclusions are well supported statistically. Interpretation should be reasonably mature. The project should demonstrate a clear overview of the work without getting lost in details and be free of minor statistical errors.
60-69 A project that could be presented after a round of revision, but without having to re-do much of the actual analysis. Some substantive flaws in the analysis or presentation (or more minor flaws in both), but sound. A good grasp of the statistics and context, so that interpretation is reasonable.
50-59 Major re-working required before the project could be presented, but containing some sound statistics demonstrating understanding of statistical methods and their application. Good presentation and organization.
40-49 Major flaws in analysis and presentation, but demonstrating some understanding of statistical methods and a reasonable attempt to present the results.
0-39 Flawed analysis demonstrating little or no understanding of statistical methods, and/or incomprehensible or very poorly organized presentation.
Questions?
Course schedule
Week 1: Introduction to Spatial Data Science with R
Week 2: Areal data: modeling
Week 3: Areal data: neighborhood matrices, spatial autocorrelation. R packages for making maps and downloading open spatial data
Week 4: Interactive dashboards. Saudi National Day
Week 5: Geostatistical data: spatial interpolation methods
Week 6: Student presentations project 1
Week 7: Geostatistical data: modeling. Spatially misaligned data
Week 8: Midsemester break. Point patterns
Week 9: Point patterns: simulation, complete spatial randomness
Week 10: Point patterns: K-function
Week 11: Student presentations project 2
Week 12: Point patterns: intensity
Week 13: Point patterns: modeling
Week 14: Interactive web applications
Week 15: Student presentations project 3
Questions?
Spatial Data Science with R
Paula Moraga, Ph.D.
Associate Professor of Statistics
King Abdullah University of Science and Technology (KAUST), Saudi Arabia
Fixed effects quantify the effects of the covariates on the disease risk
Random effects represent residual variation not explained by the covariates. \(u_i\) spatial effect to account for spatial dependence between relative risks. \(v_i\) unstructured effect to account for independent noise.
Covariates based on characteristics known to affect disease transmission (temperature, precipitation, vegetation, elevation, land cover, population, etc.). Random effects model residual variation not explained by covariates
Point patterns
Assume point pattern \(\{s_i: i=1, \ldots, n\}\) has been generated as a realization of a point process \(Z = \{Z(s): s \in \mathbb{R}^2\}\)
A point process model can be used to estimate the intensity of events, identify patterns in the distribution of the observed locations, and learn about the correlation between the locations and spatial covariates
State whether each of the following is an example of areal data, geostatistical data or a point pattern.
Geographical locations of the home addresses of a sample of cardiovascular patients in Italy.
Yesterday’s mid-day temperature at the centres of fifty cities in Germany.
Carbon dioxide levels measured at twenty locations near a road.
The positions of hair follicles on a human head.
Crop yields from each of thirteen adjacent rows of wheat in a field.
Positions of snooker balls on a snooker table after one shot has been played.
The average house price in each postcode region in Paris.
Spatial data and visualization with R
Spatial data in R
Spatial data can be represented using vector and raster data
Vector data displays points, lines and polygons, and associated information
Examples: locations of monitoring stations, road networks, municipalities
Raster data are regular grids with cells of equal size that are used to store values of spatially continuous phenomena
Examples: elevation, temperature, air pollution values
R packages: sf (vector data) and terra (raster and vector data)
sf to work with vector data
Vector data are often represented using a data format called shapefile.
unprojected or geographic: Latitude and Longitude for referencing location on the ellipsoid Earth
(decimal degrees (DD) or degrees, minutes, and seconds (DMS))
projected: Easting and Northing for referencing location on 2-dimensional representation of Earth
Common projection: Universal Transverse Mercator (UTM)
Location is given by the zone number (60 zones), hemisphere (north or south), and Easting and Northing coordinates in the zone in meters
Coordinate Reference Systems (CRS)
Most common CRSs can be specified with their EPSG (European Petroleum Survey Group) codes
EPSG 4326 refers to the geographic CRS (latitude and longitude)
Interactive documents and dashboards for communication
Reproducible documents with R Markdown
R Markdown can be used to turn our analysis into fully reproducible documents that can be shared with others. Output formats include HTML, PDF or Word. An R Markdown file is written with Markdown syntax with embedded R code, and can include narrative text, tables and visualizations
Quarto can be used to create reproducible documents using multiple languages such as R, Python and Julia. A Quarto document can be rendered as formats like PDF and Word. A Quarto document has extension .qmd and is formed of a YAML header, Markdown text, and R code chunks
Interactive dashboards with flexdashboard
flexdashboard uses R Markdown to publish a group of related data visualizations as a dashboard
Integrated nested Laplace approximation (INLA) is a computational approach to perform approximate Bayesian inference in latent Gaussian models such as generalized linear mixed models and spatial and spatio-temporal models
INLA uses a combination of analytical approximations and numerical integration to obtain approximated posterior distributions of parameters and is very fast compared to MCMC. The combination of INLA and the stochastic partial differential equation (SPDE) approaches permits to analyze point data
INLA and SPDE can be applied with the R package R-INLA
Observations belong to an exponential family with mean \(\mu_i=g^{-1}(\eta_i)\) (e.g. Gaussian, Poisson, Binomial, …) \[y_i|\boldsymbol{x},\boldsymbol{\theta} \sim \pi(y_i|x_i,\boldsymbol{\theta}),\ i=1, \ldots, n\]
\(\alpha\) intercept, \(\{\beta_k\}\)’s quantify the linear effects of covariates \(\{z_{ki}\}\) on response, \(\{f^{(j)}(\cdot)\}\)’s set of random effects defined in terms of some covariates \(\{u_{ji}\}\) (e.g. iid, rw1, ar1, CAR, …)
This approximation can be integrated numerically with respect to \(\boldsymbol{\theta}\)\[\tilde \pi(x_i|\boldsymbol{y})=\sum_k \tilde \pi(x_i|\boldsymbol{\theta_k},\boldsymbol{y})\tilde \pi(\boldsymbol{\theta_k}|\boldsymbol{y})\times \Delta_k,\
\tilde \pi(\theta_j|\boldsymbol{y})=\sum_l \tilde \pi(\boldsymbol{\theta_l^*}|\boldsymbol{y})\times \Delta_l^*\]
\(\Delta_k\) area weight corresponding to \(\boldsymbol{\theta_k}\) \(\Delta_l^*\) area weight corresponding to \(\boldsymbol{\theta_l^*}\)
INLA
The approximated posterior distributions \(\tilde \pi (x_i|\boldsymbol{y})\) can be post-processed to compute quantities of interest like posterior expectations and quantiles
Write the linear predictor as a formula object in R
formula <- y ~ x1 + x2 +f(id, model ="iid")
Random effects specified with f(). First argument id is an index vector that specifies the element of the random effect that applies to each observation, second argument model name
Call the function inla()
res <-inla(formula, family ="gaussian", data = d)
Inspect the res object wich contains the fitted model. Posteriors can be post-processed using a set of functions provided by R-INLA
summary(res)res$summary.fixed
Areal data
Areal data
Standardized Mortality Ratio (SMR) is often used to estimate disease risk
\[
SMR_i = \frac{Y_i}{E_i} = \frac{\mbox{number observed cases in area } i}{\mbox{number expected cases in area } i}\]
\(SMR_i = 1\) same number observed as expected \(SMR_i > 1\) more observed than exp. (high risk) \(SMR_i < 1\) less observed than exp. (low risk)
The expected number of cases in area \(i\) are the number of cases one would expect if population in area \(i\) behaved the way the standard pop. behaves
Standard population is considered as the whole population (all areas)
Population is stratified by several factors (e.g., age and gender)
\[E_i = \sum_{j=1}^m r_j^{(std)} n_j^{(i)}\]
\(r_j^{(std)}=\frac{\mbox{number cases}}{\mbox{population}}\) rate in stratum \(j\) in the standard population
\(n_j^{(i)}\): population in stratum \(j\) of area \(i\)
Limitations
SMRs easy to calculate but may be misleading and unreliable in areas with small populations or rare diseases. Models enable to incorporate covariates & borrow information from neighboring areas to obtain smoothed relative risks
Areal models
Model to estimate disease relative risk \(\theta_i\) in areas \(i=1,\ldots,n\)
Neighbors based on contiguity assume that neighbors of a given area are other areas that share a common boundary. Neighbors can be of type Queen if a single shared boundary point meets the contiguity condition, or Rook if more than one shared point is required to meet the contiguity condition
nb <-poly2nb(map)id <-20# neighbors of area id 20map$neighbors <-"other"; map$neighbors[id] <-"area"; map$neighbors[nb[[id]]] <-"neighbors"plot(st_geometry(map), col =c("gray30", "gray", "white")[as.factor(map$neighbors)])
Neighbors based on k nearest neighbors
coo <-st_centroid(map); nb <-knn2nb(knearneigh(coo, k =3)) # k number nearest neighborsid <-20# neighbors of area id 20map$neighbors <-"other"; map$neighbors[id] <-"area"; map$neighbors[nb[[id]]] <-"neighbors"plot(st_geometry(map), col =c("gray30", "gray", "white")[as.factor(map$neighbors)])
Neighbors based on distance
nb <-dnearneigh(x =st_centroid(map), d1 =0, d2 =1)id <-20# neighbors of area id 20map$neighbors <-"other"; map$neighbors[id] <-"area"; map$neighbors[nb[[id]]] <-"neighbors"plot(st_geometry(map), col =c("gray30", "gray", "white")[as.factor(map$neighbors)])
Neighbors of order k based on contiguity
Neighbors of order k based on contiguity
nb <-poly2nb(map, queen =TRUE)nblags <- spdep::nblag(neighbours = nb, maxlag =2)
par(mfrow =c(1, 2))nb1 <- nblags[[1]]; nb2 <- nblags[[2]] # Neighbors of first and second orderid <-20# neighbors of area id 20map$neighbors <-"other"; map$neighbors[id] <-"area"; map$neighbors[nb1[[id]]] <-"neighbors"plot(st_geometry(map), col =c("gray30", "gray", "white")[as.factor(map$neighbors)])map$neighbors <-"other"; map$neighbors[id] <-"area"; map$neighbors[nb2[[id]]] <-"neighbors"plot(st_geometry(map), col =c("gray30", "gray", "white")[as.factor(map$neighbors)])
Neighborhood matrices
A neighborhood matrix \(W\) defines a neighborhood structure over the study region, and its elements can be viewed as weights. \(w_{ij}\) spatially connects areas \(i\) and \(j\) in some fashion. More weight is associated with areas closer to \(i\) than farther away from \(i\)
Neighborhood matrices
If neighbors are based on contiguity, we can construct a binary spatial matrix with \(w_{ij} = 1\) if regions \(i\) and \(j\) share a common boundary, and \(w_{ij} = 0\) otherwise. Customarily, \(w_{ii}\) is set to 0 for \(i = 1, \ldots, n\)
\[
\begin{matrix}
& A & B & C & D & E & \text{Sum} \\
A & 0 & 1 & 1 & 1 & 0 & 3 \\
B & 1 & 0 & 1 & 1 & 1 & 4 \\
C & 1 & 1 & 0 & 1 & 0 & 3 \\
D & 1 & 1 & 1 & 0 & 1 & 4 \\
E & 0 & 1 & 0 & 1 & 0 & 2
\end{matrix}
\]
Neighborhood matrices
Other definitions could be \(w_{ij} = 1\) for all \(i\) and \(j\) within a specified distance, \(w_{ij} = 1\) if \(j\) is one of the \(k\) nearest neighbors of \(i\), or \(w_{ij}\) = inverse distance between areas
We may also want to adjust for the total number of neighbors in each area and use a standardized matrix with entries \(w_{std,i,j} = w_{ij}/\sum_{j=1}^n w_{ij}\)
\[
\begin{matrix}
& A & B & C & D & E & \text{Sum} \\
A & 0 & 1 & 1 & 1 & 0 & 3 \\
B & 1 & 0 & 1 & 1 & 1 & 4 \\
C & 1 & 1 & 0 & 1 & 0 & 3 \\
D & 1 & 1 & 1 & 0 & 1 & 4 \\
E & 0 & 1 & 0 & 1 & 0 & 2
\end{matrix}
\]
\[
\begin{matrix}
& A & B & C & D & E & \text{Sum} \\
A & 0 & 0.333 & 0.333 & 0.333 & 0 & 1 \\
B & 0.25 & 0 & 0.25 & 0.25 & 0.25 & 1 \\
C & 0.333 & 0.333 & 0 & 0.333 & 0 & 1 \\
D & 0.25 & 0.25 & 0.25 & 0 & 0.25 & 1 \\
E & 0 & 0.5 & 0 & 0.5 & 0 & 1
\end{matrix}
\]
Spatial weights matrix based on a binary neighbor list
nb <-poly2nb(map, queen =TRUE)nb[1:2]
[[1]]
[1] 2 3
[[2]]
[1] 1 3 4
# style B binary, W row standardized (1/n.neigh)nbw <- spdep::nb2listw(nb, style ="W")nbw$weights[1:2]
Spatial autocorrelation is used to describe the extent to which a variable is correlated with itself through space
Positive spatial autocorrelation occurs when observations with similar values are closer together (i.e., clustered). Negative spatial autocorrelation occurs when observations with dissimilar values are closer together (i.e., dispersed)
Global Moran’s \(I\)
Moran’s \(I\) is used to assess spatial autocorrelation in areal data. It summarizes the degree to which similar observations tend to occur near each other over the study region
Moran’s \(I\)\(>\)\(E[I] = -1/(n-1)\) indicate positive spatial autocorrelation or clustering, neighboring regions tend to have similar values
Moran’s \(I\)\(<\)\(E[I]\) indicate negative spatial autocorrelation or dispersion, regions that are close to one another tend to have different values
Moran’s \(I\) around \(E[I]\) indicate randomness, absence of spatial pattern
# Neighborsnb <- spdep::poly2nb(map, queen =TRUE) # queen shares point or bordernbw <-nb2listw(nb, style ="W")# Global Moran's Igmoran <-moran.test(map$vble, nbw, alternative ="two.sided")gmoran
Moran I test under randomisation
data: map$vble
weights: nbw
Moran I statistic standard deviate = 23.35, p-value < 2.2e-16
alternative hypothesis: two.sided
sample estimates:
Moran I statistic Expectation Variance
0.6266753872 -0.0019801980 0.0007248686
Moran’s \(I\)\(>\)\(E[I] = -1/(n-1)\) indicate positive spatial autocorrelation or clustering, neighboring regions tend to have similar values
Assessing significance
\(H_0: I = E[I]\) no spatial autocorrelation \(H_1: I \neq E[I]\) spatial autocorrelation
When the number of regions is sufficiently large, \(I\) has normal distribution and we can assess whether a pattern deviates significantly from a random pattern by comparing the z-score \(z = \frac{I-E(I)}{Var(I)^{1/2}}\) to the standard normal distribution
An alternative approach to assess significance is Monte Carlo randomization. Create random patterns by reassigning the observed values among the areas and calculate Moran’s \(I\) for each pattern, obtaining a randomization distribution for Moran’s \(I\). If the observed Moran’s \(I\) lies in the tails of this distribution, the assumption of no spatial autocorrelation is rejected
Assessing significance
\(H_0: I = E[I]\) no spatial autocorrelation \(H_1: I \neq E[I]\) positive or negative spatial autocorrelation
Assessing significance with Normal(0, 1)
State the null and alternative hypotheses: \(H_0: I = E[I]\) no spatial autocorrelation \(H_1: I \neq E[I]\) spatial autocorrelation
Choose the significance level \(\alpha\) we are willing to tolerate, which represents the maximum value for the probability of incorrectly rejecting the null hypothesis when it is true (usually \(\alpha = 0.05\)).
Calculate the test statistic: \[z = \frac{I-E(I)}{Var(I)^{1/2}}.\]
Find the p-value for the observed data by comparing the z-score to the standard normal distribution or via Monte Carlo randomization. The p-value is the probability of obtaining a test statistic as extreme as or more extreme than the one observed test statistic in the direction of the alternative hypothesis, assuming the null hypothesis is true.
Make one of these two decisions and state a conclusion:
If p-value \(< \alpha\), we reject the null hypothesis. We conclude data provide evidence for the alternative hypothesis.
If p-value \(\geq \alpha\), we fail to reject the null hypothesis. The data do not provide evidence for the alternative hypothesis.
Assessing significance with Normal(0, 1)
\(H_0: I = E[I]\) no spatial autocorrelation \(H_1: I \neq E[I]\) positive or negative spatial autocorrelation
# Neighbors# queen shares point or bordernb <- spdep::poly2nb(map, queen =TRUE)nbw <-nb2listw(nb, style ="W")# Global Moran's Igmoran <-moran.test(map$vble, nbw, alternative ="two.sided")gmoran
Moran I test under randomisation
data: map$vble
weights: nbw
Moran I statistic standard deviate = 23.35, p-value < 2.2e-16
alternative hypothesis: two.sided
sample estimates:
Moran I statistic Expectation Variance
0.6266753872 -0.0019801980 0.0007248686
gmoran[["estimate"]][["Moran I statistic"]]
[1] 0.6266754
gmoran[["statistic"]] # z-score
Moran I statistic standard deviate
23.3498
gmoran[["p.value"]] # p-value
[1] 1.384667e-120
p-value \(<\) significance level 0.05. Then, we reject the null hypothesis and conclude there is evidence for spatial autocorrelation
Assessing significance with Monte Carlo
\(H_0: I = E[I]\) no spatial autocorrelation \(H_1: I \neq E[I]\) positive or negative spatial autocorrelation
Create random patterns by reassigning values among areas and calculate Moran’s \(I\) for each pattern. p-value is proportion of values as extreme or more extreme than the statistic observed in the direction of alternative hypothesis. p-value \(<\) 0.05 indicating data presents spatial autocorrelation
Monte-Carlo simulation of Moran I
data: map$vble
weights: nbw
number of simulations + 1: 1000
statistic = 0.62668, observed rank = 1000, p-value < 2.2e-16
alternative hypothesis: two.sided
# Moran's I for simulated patternshist(gmoranMC$res)# Red line Moran's I for real dataabline(v = gmoranMC$statistic, col ="red")
Assessing significance
\(H_0: I \leq E[I]\) negative spatial autocorrelation or no spatial autocorrelation \(H_1: I > E[I]\) positive spatial autocorrelation
Assessing significance with Normal(0, 1)
\(H_0: I \leq E[I]\) negative spatial autocorrelation or no spatial autocorrelation \(H_1: I > E[I]\) positive spatial autocorrelation
# Neighbors# queen shares point or bordernb <- spdep::poly2nb(map, queen =TRUE)nbw <-nb2listw(nb, style ="W")# Global Moran's Igmoran <-moran.test(map$vble, nbw, alternative ="greater")gmoran
Moran I test under randomisation
data: map$vble
weights: nbw
Moran I statistic standard deviate = 23.35, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.6266753872 -0.0019801980 0.0007248686
gmoran[["estimate"]][["Moran I statistic"]]
[1] 0.6266754
gmoran[["statistic"]] # z-score
Moran I statistic standard deviate
23.3498
gmoran[["p.value"]] # p-value
[1] 6.923334e-121
p-value \(<\) significance level 0.05. Then, we reject the null hypothesis and conclude there is evidence for positive spatial autocorrelation
Assessing significance with Monte Carlo
\(H_0: I \leq E[I]\) negative spatial autocorrelation or no spatial autocorrelation \(H_1: I > E[I]\) positive spatial autocorrelation
Create random patterns by reassigning values among areas and calculate Moran’s \(I\) for each pattern. p-value is proportion of values as extreme or more extreme than the statistic observed in the direction of alternative hypothesis. p-value \(<\) 0.05 indicating data presents positive spatial autocorrelation
(gmoranMC <-moran.mc(map$vble, nbw, nsim =999))
Monte-Carlo simulation of Moran I
data: map$vble
weights: nbw
number of simulations + 1: 1000
statistic = 0.62668, observed rank = 1000, p-value = 0.001
alternative hypothesis: greater
# Moran's I for simulated patternshist(gmoranMC$res)# Red line Moran's I for real dataabline(v = gmoranMC$statistic, col ="red")
Moran’s \(I\) scatterplot
Moran’s \(I\) scatterplot can be used to visualize the spatial autocorrelation in the data. It displays observations of each area against its spatially lagged values
The spatially lagged value for a given area is calculated as a weighted average of the neighboring values for that area
\[(Y_i,\ \sum_{j=1}^n w_{ij}Y_j)\]
moran.plot(map$vble, nbw)
Exercise
Reproduce the Moran’s \(I\) scatterplot Spatial lagged variable for area \(i = \sum_{j=1}^n w_{ij}Y_j\) (sum(neighmat[i, ]*map$vble))
Spatial lagged variable can be obtained with spdep::lag.listw() or apply(neighmat %*% map$vble, 1, sum)
nb <-poly2nb(map, queen =TRUE)nbw <-nb2listw(nb, style ="W") # style = "W" is row standardized# spatial lagged variablevblelag <-lag.listw(nbw, map$vble) # spatial weights of areas and valuesres <-lm(vblelag ~ map$vble)plot(map$vble, vblelag, xlab ="vble", ylab ="Spatial lag of vble")abline(res$coefficients[1], res$coefficients[2])abline(h =mean(vblelag), lty =2); abline(v =mean(map$vble), lty =2)
Local Moran’s \(I\)
Global Moran’s \(I\) assesses spatial autocorrelation for the whole study region. We can also provide a local measure of similarity between each area’s value and those of nearby areas
Local Indicators of Spatial Association (LISA) provide indicate the extent of significant spatial clustering of similar values around each observation
A high value for \(I_i\) suggests that the area is surrounded by areas with similar values. Such an area is part of a cluster of high observations, low observations, or moderate observations.
A low value for \(I_i\) indicates that the area is surrounded by areas with dissimilar values. Such an area is an outlier indicating that the observation of area \(i\) is different from most or all of the observations of its neighbors.
Local Moran’s \(I\)
\(H_0\): no or negative spatial autocorrelation, \(H_1\): positive spatial autocorrelation
lmoran <-localmoran(map$vble, nbw, alternative ="greater")head(lmoran)
We identify the clusters of each type by using the information from the Moran’s \(I\) scatterplot showing the scaled values against its spatially lagged values
mp <-moran.plot(as.vector(scale(map$vble)), nbw)
Clusters
We identify cluster types by using the quadrants of the scaled values (mp$x) and their spatially lagged values (mp$wx), and p-values obtained with the local Moran’s \(I\) for each area (map$lmp). Areas with significant local Moran’s \(I\) are:
High-high if both the value and its spatially lagged value are positive
Low-low if both the value and its spatially lagged value are negative
High-low if the the value is positive and the spatially lagged value negative
Low-high is the value is negative and the spatially lagged value positive
Identify cluster types of household income in Ohio
# https://jakubnowosad.com/spData/reference/columbus.html?utm_source=chatgpt.com# x and y coordinates (in arbitrary digitizing units, not polygon coordinates)library(spData); library(sf); library(mapview)map <-st_read(system.file("shapes/columbus.gpkg", package ="spData"), quiet =TRUE)# Assign CRS so mapview can plot it. https://epsg.io/# It will not be able to transform to geographical CRS 4326 since coordinates are in arbitrary unitsst_crs(map) <-"EPSG:3631"map$vble <- map$INCmapview(map, zcol ="vble")
\(Y_i\) number people positive, \(n_i\) number people tested, \(P(\boldsymbol{x}_i)\) prevalence at \(\boldsymbol{x}_i\)
Covariates based on characteristics known to affect disease transmission (temperature, precipitation, vegetation, elevation, land cover, population density, etc.)
Random effects model residual variation not explained by covariates
Inference using INLA and SPDE
Integrated nested Laplace approximations (INLA) is a computational approach to perform approximate Bayesian inference in latent Gaussian models
In the Stochastic partial differential equation (SPDE) approach, the continuously indexed Gaussian random field \(S\) is represented as a discretely indexed Gaussian Markov random field (GMRF) by means of a finite basis function defined on a triangulation of the study region
\(\psi_g(\cdot)\) piecewise polynomial basis functions on each triangle \(\{S_g \}\) zero-mean Gaussian distributed \(G\) number of vertices in triangulation
Projection matrix
\(S(\boldsymbol{x})\) weighted average of the GMRF values at the vertices of the triangle containing the point. Weights = barycentric coordinates
\(T_1, T_2, T_3\) areas subtriangles formed by \(\boldsymbol{x}\) and vertices. \(T\) area whole triangle
Projection matrix
\[S(\boldsymbol{x}_i) \approx \sum_{g=1}^G A_{ig} S_g\]\(A\) projection matrix that maps GMRF from observations to triangulation nodes
Row \(i\) corresponding to observation \(i\) has possibly three non-zero values at columns that represent the vertices of the triangle containing the point
(non-zero values are equal to the barycentric coordinates)
Assume point pattern \(\{s_i: i=1, \ldots, n\}\) has been generated as a realization of a point process \(Z = \{Z(s): s \in \mathbb{R}^2\}\)
A point process model can be used to estimate the intensity of events, identify patterns in the distribution of the observed locations, and learn about the correlation between the locations and spatial covariates
Intensity of the process
Intensity of the process: mean number of events per unit area at location \(s\)\[\lambda(s) = lim_{|ds|\rightarrow 0} \frac{E[N(ds)]}{|ds|}\]
library(spatstat)plot(X)den <-density(x = X, sigma =10)plot(den, main ="Intensity")contour(den, add =TRUE)
Complete Spatial Randomness
Point processes are stochastic models that describe the locations of events of interest and possibly some additional information such as marks that inform about different types of events (e.g., cases and controls)
Point processes provide models for point patterns, with complete spatial randomness (CSR) being the simplest theoretical model
CSR helps differentiate between regular and clustered patterns
Homogeneous Poisson process (CSR)
Simple model that assumes events are equally likely to occur at any location within the study area, independent of the locations of other events
the number of events in any region \(A\) follows a Poisson distribution with mean \(\lambda |A|\), where \(\lambda\) is a constant value denoting the intensity and \(|A|\) is the area of \(A\)
given that there are \(n\) events inside \(A\), the locations of these events are independent and uniformly distributed in \(A\)
Inhomogeneous Poisson process
In homogeneous Poisson processes, the intensity is constant (\(\lambda(s)= \lambda,\ \ \forall s\)), whereas in inhomogeneous Poisson processes, the intensity varies in space.
the number of events in any region \(A\) follows a Poisson distribution with mean \(\mu(A) = \int_A \lambda(s)ds\)
given that there are \(n\) events inside \(A\), the locations of these events are independent with probability density function proportional to the intensity \(\lambda(\cdot)\)
Log-Gaussian Cox process (LGCP)
Log-Gaussian Cox processes (LGCP) are typically used to model phenomena that are environmentally driven
A LGCP is a Poisson process with a varying intensity which is itself a stochastic process of the form \[\Lambda(s) = exp(\eta(s)),\] where \(\eta = \{\eta(s): s \in \mathbb{R}^2\}\) is spatial Gaussian process
Conditional on \(\Lambda(s)\), the point process is a Poisson process with intensity \(\Lambda(s)\) 1. the number of events in any region \(A\) is distributed as \(Poisson(\int_A \Lambda(s)ds)\) 2. the locations of events are an independent random sample with probability density proportional to \(\Lambda(s)\)
Fitting a log-Gaussian Cox process model
Assume point pattern \(\{s_i: i=1, \ldots, n\}\) has been generated as a realization of a LGCP with intensity \(\Lambda(s)= exp(\eta(s))\)
A LGCP can be fitted by discretizing the study region into a grid with \(n_1 \times n_2\) cells \(\{s_{ij}: i=1,\ldots,n_1, j=1,\ldots,n_2\}\). \(|s_{ij}|\) area of cell \(s_{ij}\)
R Markdown can be used to turn our analysis into fully reproducible documents that can be shared with others. Output formats include HTML, PDF or Word. An R Markdown file is written with Markdown syntax with embedded R code, and can include narrative text, tables and visualizations
A Shiny app can be built by creating a directory (called, for example, appdir) that contains an R file (called, for example, app.R) with three components:
ui user interface object which controls layout and appearance of the app
server() function with instructions to build objects displayed in the ui
call to shinyApp() that creates the Shiny app from the ui/server pair
# define user interface objectui <-fluidPage( )# define server() functionserver <-function(input, output){ }# call to shinyApp() which returns the Shiny appshinyApp(ui = ui, server = server)
Save app.R inside directory appdir. Launch the app:
shiny::runApp("appdir_path")
Directory appdir can also contain data or other R scripts needed by the app. We can also write two separate files ui.R and server.R. This permits an easier management of the code for large apps
Inputs
Outputs
Plots, tables, texts, images
Inputs, outputs and reactivity
Inputs: we can interact with the app by modifying their values
Outputs: objects we want to show in the app
ui <-fluidPage(*Input(inputId = myinput, label = mylabel, ...)*Output(outputId = myoutput, ...))server <-function(input, output){ output$myoutput <- render*({# code to build the output.# If it uses an input value (input$myinput),# the output will be rebuilt whenever the input value changes })}
We can include inputs by writing an *Input() function in the ui (e.g., textInput(), dateRangeInput(), fileInput())
Outputs can be created by writing an *Output() function in the ui, and the R code to build the output inside a render*() function in server() (e.g., textOutput() and renderText(), tableOutput() and renderTable())
fluidPage() creates a display that automatically adjusts the app to the dimensions of the browser window
Inputs: we can interact with the app by modifying their values
Outputs: objects we want to show in the app
ui <-fluidPage(*Input(inputId = myinput, label = mylabel, ...)*Output(outputId = myoutput, ...))server <-function(input, output){ output$myoutput <- render*({# code to build the output.# If it uses an input value (input$myinput),# the output will be rebuilt whenever the input value changes })}
We can include inputs by writing an *Input() function in the ui (e.g., textInput(), dateRangeInput(), fileInput())
Outputs can be created by writing an *Output() function in the ui, and the R code to build the output inside a render*() function in server() (e.g., textOutput() and renderText(), tableOutput() and renderTable())
fluidPage() creates a display that automatically adjusts the app to the dimensions of the browser window
Build a Shiny app containing the following inputs and outputs:
Slider input with label "Select dates:", minimum date "2020-01-01", maximum date "2020-12-31", and initial value "2020-01-31".
Slider input to select a value between 0 and 100 where the interval between each selectable value on the slider is 5. Add animation so when the user presses play the input widget scrolls through automatically.
Output with a text showing the input values selected.
Exercise
library(shiny)ui <-fluidPage(sliderInput("date", "Select date:",min =as.Date("2020-01-01"),max =as.Date("2020-12-31"),value =as.Date("2020-01-31")),sliderInput("number", "Select a number:",min =0, max =100, value =0,step =5, animate =TRUE),textOutput("inputvaluesselected"))server <-function(input, output, session) { output$inputvaluesselected <-renderText({paste("Date", input$date, ",","number", input$number, ".") })}shinyApp(ui, server)
# Function to calculate nth number in Fibonacci sequencefib <-function(n){ifelse(n<3, 1, fib(n-1)+fib(n-2))}ui <-fluidPage(numericInput("n", "Enter number", value =1),textOutput("nthvalue"),textOutput("nthvalueInv") )server <-function(input, output) {# time consuming calculation that we only want to do once# do calculation in a reactive expression currentFib <-reactive({fib(as.numeric(input$n))}) output$nthvalue <-renderText({currentFib()}) output$nthvalueInv <-renderText({1/currentFib()})}shinyApp(ui, server)
Observers
Observers, like reactive expressions, can take reactive values and reactive expressions. However, they do not return any values. Instead of returning values they have side effects
ui <-fluidPage(sliderInput("controller", "Controller", 0, 20, 10),textInput("inText", "Input text"),textInput("inText2", "Input text 2"))server <-function(input, output, session) {observe({# We use input$controller multiple times so we save it as x x <- input$controller# changes the value of input$inTextupdateTextInput(session, "inText", value =paste("New text", x))# changes the value and label of input$inText2updateTextInput(session, "inText2",label =paste("New label", x),value =paste("New text", x)) })}shinyApp(ui, server)
Sometimes we may want to wait for a specific action to be taken from the user, like clicking an actionButton() or an actionLink(), before calculating an expression or taking an action
An event is a reactive value or expression that is used to trigger other calculations in this way
observeEvent() and eventReactive() can be used for event handling:
observeEvent() to perform an action in response to an event
eventReactive() to create a calculated value that only updates in response to an event
ui <-fluidPage(numericInput("number", "Value", 5),actionButton("button", "Show"),tableOutput("table"))server <-function(input, output) {# observeEvent() and eventReactive()# take a dependency on input$button, but# not on any of the stuff inside the functions# Take an action every time button is pressed # Here, print a message to the consoleobserveEvent(input$button, {cat("Showing",input$number,"rows\n")})# Calculate df every time button is pressed df <-eventReactive(input$button, {head(cars, input$number)}) output$table <-renderTable({df()})}shinyApp(ui, server)
fileInput() returns a data frame with four columns: name, size, type, datapath
Uploads
input$upload initialized to NULL on page load. Need req(input$upload) to make sure code waits until first file is uploaded
accept argument limits the possible inputs. But it is only a suggestion to browser not always enforced. Good practice validate (tools::file_ext())
ui <-fluidPage(fileInput("upload", "Select file", accept =c(".csv", ".tsv")),numericInput("n", "Rows", value =5, min =1, step =1),tableOutput("head"))server <-function(input, output, session) { data <-reactive({req(input$upload) ext <- tools::file_ext(input$upload$name)switch(ext, # file extension excluding leading dotcsv = vroom::vroom(input$upload$datapath, delim =","),tsv = vroom::vroom(input$upload$datapath, delim ="\t"),validate("Invalid file; Please upload a .csv or .tsv file")) }) output$head <-renderTable({head(data(), input$n)})}shinyApp(ui, server)
Downloads
In ui, downloadButton(id) or downloadLink(id) so user can click to download a file. In server(), downloadHandler() with two fn arguments
ui <-fluidPage(selectInput("dataset", "Pick a dataset", ls("package:datasets")),tableOutput("preview"),downloadButton("download", "Download .tsv"))server <-function(input, output, session) { data <-reactive({ out <-get(input$dataset, "package:datasets")if(!is.data.frame(out)){validate(paste0("'", input$dataset, "' is not a data frame"))} out }) output$preview <-renderTable({head(data())}) output$download <-downloadHandler(# returns a file name shown in the download dialog boxfilename =function(){paste0(input$dataset, ".tsv")},# argument file is the path to save the file. Saves the file# in a place Shiny knows about, so it can send it to the usercontent =function(file){vroom::vroom_write(data(), file)})}shinyApp(ui, server)
library(leaflet)ui <-fluidPage(leafletOutput("map"),actionButton("button", "New locations"))server <-function(input, output, session){# when button is pressed, generate 40 new locations points <-eventReactive(input$button, {cbind(rnorm(40)*2+13, rnorm(40)+48)}, ignoreNULL =FALSE) output$map <-renderLeaflet({leaflet() %>%addTiles() %>%addMarkers(data =points()) })}shinyApp(ui, server)
Set ignoreNULL = FALSE to initially perform the calculation and let the user recalculate (reference)
Leaflet
This works, but reactive inputs and expressions that affect the renderLeaflet() expression will cause the entire map to be redrawn from scratch and reset the map position and zoom level
Use leafletProxy() instead of leaflet() to have more control over the map, such as changing the color of a single polygon or adding a marker at the point of a click without redrawing the entire map
Typically, use leaflet() to create the static aspects of the map, and leafletProxy() to modify the dynamic elements of a map that is already running in the page
Leaflet
library(leaflet)ui <-fluidPage(leafletOutput("map"),actionButton("button", "New locations"))server <-function(input, output, session){# when button is pressed, generate 40 new locations points <-eventReactive(input$button, {cbind(rnorm(40)*2+13, rnorm(40)+48)}, ignoreNULL =FALSE)# leaflet() to include aspects that do not change in the map output$map <-renderLeaflet({leaflet() %>%addTiles() %>%setView(lng =15, lat =48, zoom =6) })# leafletproxy() to include aspects that change in the mapobserve({leafletProxy("map") %>%clearMarkers() %>%addMarkers(data =points()) })}shinyApp(ui, server)
Shiny practical 3
Interactive dashboards with flexdashboard and Shiny
Interactive dashboards with flexdashboard and Shiny
# leaflet() to include aspects that do not change in the mapoutput$map <-renderLeaflet({leaflet() %>%addTiles() %>%setView(lng =0, lat =30, zoom =2)})# leafletproxy() to include aspects that change in the mapobserve({if (nrow(mapFiltered()) ==0) {return(NULL)}leafletProxy("map", session, data =mapFiltered()) %>%clearShapes() %>%clearControls() %>%addPolygons(fillColor =~pal(PM2.5), color ="white",fillOpacity =0.7, label =~labels,highlight =highlightOptions(color ="black",bringToFront =TRUE)) %>% leaflet::addLegend(pal = pal, values =~PM2.5, opacity =0.7, title ="PM2.5")})leafletOutput("map")
Interactive dashboards with shinyflexdashboard
Package shinydashboard provides another way to create dashboards with Shiny
Before deployApp(), use setwd() to set the working directory to the folder where app.R is. Or use the argument appDir of deployApp() to set the folder where app.R is
Ribeiro Amaral, A., et al. (2022). Spatio-temporal modeling of infectious diseases by integrating compartment and point process models. SERRA, 37, 1519-1533
Moraga, P. and Baker, L. (2022). rspatialdata: a collection of data sources and tutorials on downloading and visualising spatial data using R. F1000Research, 11:77
Moraga, P., et al. (2019). epiflows: an R package for risk assessment of travel-related spread of disease. F1000Research, 7:1374
Moraga, P. (2018). Small Area Disease Risk Estimation and Visualization Using R. The R Journal, 10(1):495-506
Moraga, P. (2017). SpatialEpiApp: A Shiny Web Application for the analysis of Spatial and Spatio-Temporal Disease Data. Spatial and Spatio-temporal Epidemiology, 23:47-57
Moraga, P., et al. (2017). A geostatistical model for combined analysis of point-level and area-level data using INLA and SPDE. Spatial Statistics, 21:27-41
Moraga, P. and Kulldorff, M. (2016). Detection of spatial variations in temporal trends with a quadratic function. Statistical Methods for Medical Research, 25(4):1422-1437
Hagan, J. E., Moraga, P., et al. (2016). Spatio-temporal determinants of urban leptospirosis transmission: Four-year prospective cohort study of slum residents in Brazil. PLOS Neglected Tropical Diseases, 10(1): e0004275
Moraga, P., et al. (2015). Modelling the distribution and transmission intensity of lymphatic filariasis in sub-Saharan Africa prior to scaling up interventions: integrated use of geostatistical and mathematical modelling. Parasites & Vectors, 8:560
