0.1 Patrones puntuales
Point Pattern Analysis responde a preguntas sobre por qué las cosas aparecen donde aparecen. Las cosas podrían ser árboles, casos de enfermedades, crímenes, rayos - cualquier cosa con una ubicación puntual.
0.1.0.1 Proceso generador del patron puntual
0.1.0.2 Inferencias sobre el proceso
0.1.1 Quadrant Test (Conteo en “cuadrantes”)
Los humanos tendemos a ver patrones en arreglos aleatorios, así que necesitamos pruebas estadísticas. La prueba de recuento de cuadrantes fue uno de los primeros métodos estadísticos espaciales desarrollados. Se puede utilizar para comprobar si los puntos son completamente aleatorios desde el punto de vista espacial, es decir, si son uniformemente aleatorios en toda el área de interés.
La ventana de observación se divide en baldosas, y se cuenta el número de puntos de datos de cada baldosa, tal y como se describe en el recuento cuádruple. Los cuadrantes son rectangulares por defecto, o pueden ser regiones de forma arbitraria especificadas por el argumento tess. También se calcula el número esperado de puntos en cada cuadrante, determinado por el CSR (Complete Spatial Randomness)
#install.packages("spatstat")
library(spatstat)
ls.str()## ppxy : List of 5
## $ window :List of 5
## $ n : int 300
## $ x : num [1:300] -2.03 7.83 -7.57 -7.41 8.81 ...
## $ y : num [1:300] 2.69 -0.942 -1.888 2.72 -2.902 ...
## $ markformat: chr "none"# Set coordinates and window
# ppxy <- ppp(x = x, y = y, window = disc(radius))
# Test the point pattern
qt <- quadrat.test(ppxy)## Warning: Some expected counts are small; chi^2 approximation may be
## inaccurate# Inspect the results
plot(qt)
print(qt)##
## Chi-squared test of CSR using quadrat counts
## Pearson X2 statistic
##
## data: ppxy
## X2 = 20.617, df = 24, p-value = 0.6776
## alternative hypothesis: two.sided
##
## Quadrats: 25 tiles (irregular windows)El valor p de la prueba del cuadrante es mucho mayor que 0,05, por lo que no se puede rechazar la hipótesis de que los puntos son completamente aleatorios desde el punto de vista espacial.
0.1.2 Desventajas
- El tamaño de la partición afecta los resultados
Alternativas, nearest Neighbors (distancia al vecino más cercano) y Ripley’s K Function (Individuos un una bola de radio D)
0.1.3 Creación de un patrón de puntos uniforme con spatstat
Un proceso de punto Poisson crea eventos de acuerdo a una distribución de Poisson con un parámetro de intensidad que especifica los eventos esperados por unidad de área. El número total de eventos generados es un único número de una distribución de Poisson, por lo que múltiples realizaciones del mismo proceso pueden tener fácilmente diferentes números de eventos.
En el ejercicio anterior se utilizó un conjunto de 300 eventos dispersos uniformemente dentro de un círculo. Si repite la generación de los eventos de nuevo, todavía tendrá 300 de ellos, pero en diferentes lugares. El conjunto de datos de exactamente 300 puntos es de un proceso de punto de Poisson condicionado a que el total sea de 300.
El paquete spatstat puede generar procesos espaciales de Poisson con la función rpoispp() dada una intensidad y una ventana, que no están condicionadas al total.
Así como las funciones del generador de números aleatorios en R comienzan con una “r”, la mayoría de las funciones de patrón de puntos aleatorios en spatstat comienzan con una “r”.
La función area() de spatstat calculará el área de una ventana como un disco.
0.1.4 Simulación de patrones agrupados e inhibitorios
El paquete spatstat también tiene funciones para generar patrones de puntos a partir de otros modelos de proceso. Estos generalmente se dividen en dos clases: procesos agrupados, donde los puntos ocurren juntos más que bajo un proceso Poisson uniforme, y procesos regulares (también conocidos como inhibidores) donde los puntos están más separados que bajo un proceso Poisson de intensidad uniforme. Algunos modelos de proceso pueden generar patrones en un continuo, desde agrupados, pasando por uniformes, hasta regulares, dependiendo de sus parámetros.
La función quadrat.test() puede probar contra hipótesis agrupadas o alternativas regulares. Por defecto prueba contra cualquiera de ellos, pero esto se puede cambiar con el parámetro alternativo para crear una prueba unilateral.
Un proceso Thomas es un patrón agrupado donde un número de puntos “padre”, distribuidos uniformemente, crean un número de puntos “hijo” en su vecindario. El niño se señala a sí mismo para formar el patrón. Este es un patrón de puntos atractivo, y tiene sentido para modelar cosas como árboles, ya que los nuevos árboles crecerán cerca del árbol original. Los patrones de puntos aleatorios de Thomas pueden ser generados usando rThomas(). Esto toma tres números que determinan la intensidad y agrupación de los puntos, y un objeto ventana.
Por el contrario, los puntos de un proceso de Strauss causan una disminución en la probabilidad de encontrar otro punto cercano. Los parámetros de un proceso Strauss pueden ser tales que se trata de un proceso “hard-core”, en el que no puede haber dos puntos más cercanos que un umbral establecido. La creación de puntos a partir de este proceso implica algunos algoritmos de simulación inteligentes. Este es un patrón de puntos repulsivo, y tiene sentido para modelar cosas como animales territoriales, ya que los otros animales de esa especie evitarán el territorio de un animal dado. Se pueden generar patrones aleatorios de puntos Strauss usando rStrauss(). Esto toma tres números que determinan la intensidad y el “territorio” de los puntos, y un objeto ventana. Los puntos generados por un proceso Strauss a veces se denominan espaciados regularmente.
# Create a disc of radius 10
disc10 <- disc(10)
# Generate clustered points from a Thomas process
set.seed(123)
p_cluster <- rThomas(kappa = 0.35, scale = 1, mu = 3, win = disc10)
plot(p_cluster)
# Run a quadrat test
quadrat.test(p_cluster, alternative = "clustered")## Warning: Some expected counts are small; chi^2 approximation may be
## inaccurate##
## Chi-squared test of CSR using quadrat counts
## Pearson X2 statistic
##
## data: p_cluster
## X2 = 53.387, df = 24, p-value = 0.0005142
## alternative hypothesis: clustered
##
## Quadrats: 25 tiles (irregular windows)# Regular points from a Strauss process
set.seed(123)
p_regular <- rStrauss(beta = 2.9, gamma = 0.025, R = .5, W = disc10)## Warning: Simulation will be performed in the containing rectangle and
## clipped to the original window.plot(p_regular)
# Run a quadrat test
quadrat.test(p_regular, alternative = "regular")## Warning: Some expected counts are small; chi^2 approximation may be
## inaccurate##
## Chi-squared test of CSR using quadrat counts
## Pearson X2 statistic
##
## data: p_regular
## X2 = 8.57, df = 24, p-value = 0.001619
## alternative hypothesis: regular
##
## Quadrats: 25 tiles (irregular windows)0.2 Patrones Puntuales Bivariados
Aunque la aleatoriedad espacial es fácil de entender y probar, en la realidad es poco frecuente. Piense por ejemplo en la ubicación de las personas en el espacio, los humanos tendemos a aglomerarnos en ciudades 
Por lo cual muchos de los mapas de calor (HeatMaps) que se presentas sobre distintas variables no son más que representaciones de la densidad de las ciudades.
Un mapa de frecuencia de accidentalidad vial puede solo representar donde las personas cruzan la calle. Suponga que tiene las ubicaciones de enfermedades de las personas, es posible concluir que los casos se deban a la fuente de agua y no solamente a que la densidad de personas en esa zona es mayor?
En este caso necesitariamos, no solo la ubicación de las personas enfermas, sino de las personas sanas, en este caso tenemos un caso de patrón puntual bivariado, también podría ser multivariado.
0.2.1 Crimen en PReston
El conjunto de datos que utilizaremos para este ejemplo consiste en crímenes en un radio de 4 km del centro de Preston, una ciudad del noroeste de Inglaterra. Queremos buscar focos de crímenes violentos en el área.
Se ha construido un objeto ppp llamado preston_crime. Se trata de un proceso de puntos marcados, en el que cada punto se marca como un “delito violento” o como un “delito no violento”. Las marcas de cada punto se pueden recuperar utilizando la función marks(). La ventana es un círculo de 4 km centrado en el centro de la ciudad.
# Load the spatstat package
library(spatstat)
# Get some summary information on the dataset
summary(preston_crime)
# Get a table of marks
table(marks(preston_crime))
# Define a function to create a map
preston_map <- function(cols = c("green", "red"), cex = c(1, 1), pch = c(1, 1)) {
plotRGB(preston_osm) # from the raster package
plot(preston_crime, cols = cols, pch = pch, cex = cex, add = TRUE, show.window = TRUE)
}
# Draw the map with colors, sizes and plot characters
preston_map(
cols = c("black", "red"),
cex = c(0.5, 1),
pch = c(19, 19)
)
0.2.2 Estimación de la proporción de delitos violentos
Un método para calcular una superficie de intensidad suave a partir de un conjunto de puntos es usar el suavizado del núcleo. Imagine reemplazar cada punto con un punto de tinta sobre papel absorbente. Cada gota de tinta individual se extiende en un parche con un centro oscuro, y múltiples gotas se suman y hacen que el papel sea aún más oscuro. Con la cantidad correcta de tinta en cada gota, y con un papel de la absorbencia correcta, puede crear una impresión justa de la densidad de los puntos originales. En la jerga de suavizado del kernel, esto significa calcular un ancho de banda y usar una función particular del kernel.
Para obtener un mapa sin problemas de la proporción de crímenes violentos, podemos estimar la intensidad de la superficie de los crímenes violentos y no violentos, y tomar la proporción. Para hacer esto con la función density() en spatstat, tenemos que dividir los puntos de acuerdo a los dos valores de las marcas y luego calcular la relación entre la superficie del crimen violento y el total. La función tiene valores por defecto razonables para la función del núcleo y el ancho de banda para garantizar algo que parezca al menos plausible.
# preston_crime has been pre-defined
preston_crime
# Use the split function to show the two point patterns
crime_splits <- split(preston_crime)
# Plot the split crime
plot(crime_splits)
# Compute the densities of both sets of points
crime_densities <- density(crime_splits)
# Calc the violent density divided by the sum of both
frac_violent_crime_density <- crime_densities[[2]] /
(crime_densities[[1]] + crime_densities[[2]])
# Plot the density of the fraction of violent crime
plot(frac_violent_crime_density)
0.2.3 Selección de ancho de banda
Podemos obtener una medida más basada en principios de la proporción de delitos violentos utilizando un modelo de segregación espacial. El paquete spatialkernel implementa la teoría de la segregación espacial.
El primer paso es calcular el ancho de banda óptimo para el alisamiento del núcleo bajo el modelo de segregación. Un ancho de banda pequeño resultaría en una densidad que es casi nula, con picos en las ubicaciones de los eventos. Un gran ancho de banda aplanaría cualquier estructura en los eventos, resultando en una gran “mancha” a lo largo de toda la ventana. En algún punto entre estos extremos hay un ancho de banda que representa mejor una densidad subyacente para el proceso.
spseg() escaneará sobre un rango de anchos de banda y calculará una estadística de prueba usando un método de validación cruzada. El ancho de banda que maximiza esta estadística de prueba es el que se debe utilizar. El valor devuelto por spseg() en este caso es una lista, con elementos h y cv dando los valores de la estadística sobre los valores de entrada h. El paquete spatialkernel proporciona una función plotcv para mostrar cómo varía el valor de la prueba. El elemento hcv tiene el valor del mejor ancho de banda.
# Scan from 500m to 1000m in steps of 50m
bw_choice <- spseg(
preston_crime,
h = seq(500, 1000, by = 50),
opt = 1)
# Plot the results and highlight the best bandwidth
plotcv(bw_choice); abline(v = bw_choice$hcv, lty = 2, col = "red")
# Print the best bandwidth
print(bw_choice$hcv)
0.2.4 Probabilidades de segregación
El segundo paso es calcular las probabilidades de crímenes violentos y no violentos como una superficie plana, así como los valores p para una prueba de segregación por puntos. Esto se hace llamando a spseg() con opt = 3 y un parámetro de ancho de banda fijo h.
Normalmente se ejecutaría este proceso para al menos 100 simulaciones, pero eso tardará demasiado tiempo en ejecutarse aquí. En su lugar, ejecute sólo 10 simulaciones. A continuación, puede utilizar un segmento de objeto precargado que es la salida de una ejecución de simulación de 1000 que tardó unos 20 minutos en completarse.
# Set the correct bandwidth and run for 10 simulations only
seg10 <- spseg(
pts = preston_crime,
h = 800,
opt = 3,
ntest = 10,
proc = FALSE)
# Plot the segregation map for violent crime
plotmc(seg10, "Violent crime")
# Plot seg, the result of running 1000 simulations
plotmc(seg, "Violent crime")
0.2.5 Mapeo de la segregación
Con un mapa base y algunas funciones de imagen y contorno podemos mostrar tanto las probabilidades como las pruebas de significación sobre el área con más control que la función plotmc().
El objeto seg es una lista con varios componentes. Las coordenadas X e Y de la cuadrícula se almacenan en los elementos $gridx y $gridy. Las probabilidades de cada clase de datos (delitos violentos o no violentos) están en un elemento de la matriz \(p\) con una columna para cada clase. El valor p de la prueba de significación se encuentra en un elemento de matriz similar llamado $stpvalue. La reordenación de las columnas de estas matrices en una cuadrícula de valores se puede hacer con la función matrix() de R. Desde allí se pueden construir objetos de lista con un vector $x de coordenadas X, $y de coordenadas Y y $z como matriz. A continuación, puede alimentar a image() o a contour() para su visualización.
Este proceso puede parecer complejo, pero recuerde que con R siempre puede escribir funciones para realizar tareas complejas y las que puede repetir con frecuencia. Por ejemplo, para ayudar con el mapeo en este ejercicio, usted creará una función que construye un mapa a partir de cuatro elementos diferentes.
Se carga el objeto seg de 1000 simulaciones, así como los puntos preston_crime y la imagen del mapa preston_osm.
# Inspect the structure of the spatial segregation object
str(seg)
# Get the number of columns in the data so we can rearrange to a grid
ncol <- length(seg$gridx)
# Rearrange the probability column into a grid
prob_violent <- list(x = seg$gridx,
y = seg$gridy,
z = matrix(seg$p[, "Violent crime"],
ncol = ncol))
image(prob_violent)
# Rearrange the p-values, but choose a p-value threshold
p_value <- list(x = seg$gridx,
y = seg$gridy,
z = matrix(seg$stpvalue[, "Violent crime"] < 0.05,
ncol = ncol))
image(p_value)
# Create a mapping function
segmap <- function(prob_list, pv_list, low, high){
# background map
plotRGB(preston_osm)
# p-value areas
image(pv_list,
col = c("#00000000", "#FF808080"), add = TRUE)
# probability contours
contour(prob_list,
levels = c(low, high),
col = c("#206020", "red"),
labels = c("Low", "High"),
add = TRUE)
# boundary window
plot(Window(preston_crime), add = TRUE)
}
# Map the probability and p-value
segmap(prob_violent, p_value, 0.05, 0.15)
0.2.6 Space-Time Data
http://plotkml.r-forge.r-project.org/bigfoot.html
El sasquatch, o “bigfoot”, es una gran criatura parecida a un simio que vive en los bosques de Norteamérica. La Organización de Investigadores de Campo de Pie Grande mantiene una base de datos de avistamientos y permite su uso para la enseñanza y la investigación. Se ha creado un subconjunto limpio de datos en el noroeste de EE.UU. como el objeto ppp sasq y se carga para que usted explore el patrón espacio-temporal de avistamientos en el área.
library(sp)
data(bigfoot)
aea.prj <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96
+x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"
library(sp)
coordinates(bigfoot) <- ~Lon+Lat
proj4string(bigfoot) <- CRS("+proj=latlon +datum=WGS84")
library(rgdal)
bigfoot.aea <- spTransform(bigfoot, CRS(aea.prj))
# Load the covariates:
data(USAWgrids)
gridded(USAWgrids) <- ~s1+s2
proj4string(USAWgrids) <- CRS(aea.prj)
# Visualize data:
data(SAGA_pal)
pnts <- list("sp.points", bigfoot.aea, pch="+", col="yellow")
spplot(USAWgrids[2], col.regions=rev(SAGA_pal[[3]]), sp.layout=pnts)library(spatstat)
# Get a quick summary of the dataset
summary(sasq)## Marked planar point pattern: 423 points
## Average intensity 2.097156e-09 points per square unit
##
## *Pattern contains duplicated points*
##
## Coordinates are given to 1 decimal place
## i.e. rounded to the nearest multiple of 0.1 units
##
## Mark variables: date, year, month
## Summary:
## date year month
## Min. :1990-05-03 Y2004 : 41 Sep : 59
## 1st Qu.:2000-04-30 Y2000 : 36 Oct : 56
## Median :2003-11-05 Y2002 : 30 Aug : 54
## Mean :2003-08-11 Y2005 : 30 Jul : 50
## 3rd Qu.:2007-11-02 Y2001 : 26 Nov : 43
## Max. :2016-04-05 Y2008 : 26 Jun : 41
## (Other):234 (Other):120
##
## Window: polygonal boundary
## single connected closed polygon with 64 vertices
## enclosing rectangle: [368187.8, 764535.6] x [4644873, 5434933] units
## Window area = 2.01702e+11 square units
## Fraction of frame area: 0.644# Plot unmarked points
plot(unmark(sasq))
# Plot the points using a circle sized by date
plot(sasq, which.marks = "date")
0.2.7 Patrón espacial de avistamientos de Pie Grande
¿El patrón de avistamientos se ve regular, homogéneo o agrupado en la región de estudio?
Estos términos se definieron en Simulación de patrones agrupados e inhibidores.
quadrat.test(sasq,alternative = "clustered")## Warning: Some expected counts are small; chi^2 approximation may be
## inaccurate##
## Chi-squared test of CSR using quadrat counts
## Pearson X2 statistic
##
## data: sasq
## X2 = 261.32, df = 22, p-value < 2.2e-16
## alternative hypothesis: clustered
##
## Quadrats: 23 tiles (irregular windows)0.2.8 Patrón temporal de avistamientos de pie grande
Una vez establecido que hay algún tipo de agrupamiento espacial, es necesario explorar el comportamiento temporal. ¿Está aumentando el número de avistamientos? ¿Disminuyendo? ¿Varía la tasa a lo largo de un año (“estacionalidad”)? ¿Cambia mucho el patrón espacial a lo largo de un año?
La función base R hist() tiene un método para fechas que permite especificar una unidad de tiempo para las pausas. Pasa una cadena al argumento de rupturas, como “días”, “semanas”, “meses”, “cuartos” o “años”.
# Show the available marks
names(marks(sasq))## [1] "date" "year" "month"# Histogram the dates of the sightings, grouped by year
hist(marks(sasq)$date, "years", freq = TRUE)
# Plot and tabulate the calendar month of all the sightings
plot(table(marks(sasq)$month))
# Split on the month mark
sasq_by_month <- split(sasq, "month", un = TRUE)
# Plot monthly maps
plot(sasq_by_month)
# Plot smoothed versions of the above maps.
plot(density(sasq_by_month))
0.3 Agrupación espacio-temporal
Para hacer una prueba de agrupación espacio-temporal con stmctest() desde el paquete splancs, primero necesita convertir partes de su objeto ppp. Las funciones en los splancs tienden a utilizar datos de matriz en lugar de marcos de datos.
# Get a matrix of event coordinates
sasq_xy <- as.matrix(coords(sasq))
# Check the matrix has two columns
dim(sasq_xy)## [1] 423 2# Get a vector of event times
sasq_t <- marks(sasq)$date
# Extract a two-column matrix from the ppp object
sasq_poly <- as.matrix(as.data.frame(Window(sasq)))
dim(sasq_poly)## [1] 64 2# Set the time limit to 1 day before and 1 day after the range of times
tlimits <- range(sasq_t) + c(-1, 1)
# Scan over 400m intervals from 100m to 20km
s <- seq(100, 20000, by = 400)
# Scan over 14 day intervals from one week to 31 weeks
tm <- seq(7, 217, by = 14)0.4 Pruebas de monteCarlo
El método de Montecarlo proporciona soluciones aproximadas a una gran variedad de problemas matemáticos posibilitando la realización de experimentos con muestreos de números pseudoaleatorios en una computadora. El método es aplicable a cualquier tipo de problema, ya sea estocástico o determinista. A diferencia de los métodos numéricos que se basan en evaluaciones en N puntos en un espacio M-dimensional para producir una solución aproximada, el método de Montecarlo tiene un error absoluto de la estimación que decrece como \({\displaystyle {\frac {1}{\sqrt {N}}}}\)
0.4.1 Ejemplo
Los programas de diseño asistido por ordenador (CAD) pueden determinar rápidamente el volumen de modelos muy complejos. Estos modelos, en general, no tienen una expresión analítica para determinar su volumen (por ejemplo, para un prisma, área de la base multiplicada por la altura), y la única solución es dividir el modelo en un conjunto de pequeños submodelos (teselación) cuyo volumen pueda determinarse (por ejemplo, dividir el modelo en miles de tetraedros). Sin embargo, esto consume muchos recursos, tanto para la teselación como para el cálculo del volumen de cada uno de los elementos. Por ello utilizan métodos de Montecarlo, más robustos y eficientes.
Como el software sí que conoce la expresión analítica de la geometría del modelo (posición de los nodos, aristas y superficies) puede determinar si un punto está dentro del modelo o está fuera con un coste mucho menor que el de determinar un volumen.
En primer lugar el software coloca el modelo dentro de un volumen conocido (por ejemplo, dentro de un cubo de 1 m3 de volumen).
A continuación, genera un punto aleatorio del interior del volumen conocido, y registra si el punto “ha caído” dentro o fuera del modelo. Esto se repite un gran número de veces (miles o millones), consiguiendo un registro muy grande de cuántos puntos han quedado dentro y cuántos fuera.
Como la probabilidad de que caiga dentro es proporcional al volumen del modelo, la proporción de puntos que han caído dentro con respecto al total de puntos generados es la misma proporción de volumen que ocupa el modelo dentro del cubo de 1 \(m^3\).
Si el 50% de los puntos han caído dentro, el modelo ocupa el 50% el volumen total, es decir, 0.5 m3. Evidentemente, cuantos más puntos genere el software, menor será el error de la estimación del volumen.
0.4.2 Prueba de Monte-Carlo de la agrupación espacio-temporal
Ahora todo está listo para que pueda ejecutar la función de prueba de agrupación espacio-temporal. Luego puede graficar los resultados y calcular un valor p para rechazar la hipótesis nula de no agrupación espacio-temporal.
Cualquier agrupación espacio-temporal en un conjunto de datos se eliminará si reordena aleatoriamente las fechas de los puntos de datos. La función stmctest() calcula una estadística de prueba de agrupamiento para sus datos basada en la función K de espacio-tiempo - cuántos puntos están dentro de una ventana espacial y temporal de un punto de los datos. A continuación, realiza una serie de reorganizaciones aleatorias de las fechas entre los puntos y calcula la estadística de agrupación. Después de hacer esto un gran número de veces, puede comparar el estadístico de prueba para sus datos con los valores de los datos aleatorios. Si el estadístico de la prueba para sus datos es suficientemente grande o pequeño, puede rechazar la hipótesis nula de ninguna agrupación espacio-temporal.
La salida de stmctest() es una lista con una sola t0 que es la estadística de prueba para sus datos, y un vector de t de las simulaciones. Convirtiéndolo en un data frame puede alimentarlo a las funciones ggplot.
Debido a que el área de la ventana es un gran número de metros cuadrados, y tenemos alrededor de 400 eventos, el valor numérico de la intensidad es un número muy pequeño. Esto hace que los valores de las diversas funciones K sean muy grandes, ya que son proporcionales a la inversa de la intensidad. No se preocupe si ve 10^10 o más!
El valor p de una prueba de Monte-Carlo como ésta es sólo la proporción de estadísticas de prueba que son mayores que el valor de los datos. Esto se puede calcular a partir de los elementos t y t0 de la salida.
# Run 999 simulations
# install.packages("splancs")
library(splancs)## Loading required package: sp##
## Spatial Point Pattern Analysis Code in S-Plus
##
## Version 2 - Spatial and Space-Time analysislibrary(ggplot2)
sasq_mc <- stmctest(sasq_xy, sasq_t, sasq_poly, tlimits, s, tm, nsim = 999, quiet = TRUE)
names(sasq_mc)## [1] "t0" "t"# Histogram the simulated statistics and add a line at the data value
ggplot(data.frame(sasq_mc), aes(x = t)) +
geom_histogram(binwidth = 1e13) +
geom_vline(aes(xintercept = t0))
# Compute the p-value as the proportion of tests greater than the data
sum(sasq_mc$t > sasq_mc$t0) / 1000## [1] 0.052