Leuke128_exp

Getting Start

Understand “county_leukemia_128”

head(county_leukemia_128)

Simple feature collection with 6 features and 19 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 63.42027 ymin: -107.3468 xmax: 168.2946 ymax: 64.20241
Geodetic CRS:  WGS 84
  id leuk_cases     AverageAge GreenNeighborhoodIndex  CarsPerHH    PopDens PctPopGrowth Dist2FracLoc_1 Dist2FracLoc_2
1  1         33   50 and older             0.55333137 0.62988964 0.15221254  -0.80870577       6.274089       7.075374
2  2         41   50 and older             0.09487596 0.20972911 0.16145311   0.07389212       6.285877       4.642242
3  3         41          15-50             0.51076218 0.02426262 0.02274151  -0.55979462       6.930900       3.074290
4  4         38   50 and older             0.22416091 0.34061953 0.16278211   0.54392735       7.172685       5.044313
5  5         40 15 and younger             0.30390560 0.51989769 0.33524316  -0.59737725       6.498872       4.031625
6  6         38          15-50             0.21691189 0.78509375 0.39659359  -2.62724807       4.481083       7.930914
  Dist2FracLoc_3 Dist2FracLoc_4 Dist2FracLoc_5 Dist2FracLoc_6 Dist2FracLoc_7 Dist2FracLoc_8 Dist2FracLoc_9 Dist2FracLoc_10
1       6.872283       7.252645       2.956615      7.0120166      0.9967184       2.764238      0.2184568        3.270855
2       5.240281       5.767233       1.207414      7.1390090      3.8212552       5.500097      2.7641495        6.296465
3       1.652301       6.552304       8.038954      0.6169963      5.4261200       3.570710      6.6205182        3.655810
4       5.352336       7.079574       1.945761      6.7390486      2.3888972       4.099856      1.6237244        4.887874
5       4.669190       5.737782       1.711124      6.6417114      4.1801288       5.701428      3.2698725        6.542110
6       8.318726       5.189049       1.904093      8.9759740      3.1181019       4.646180      1.7453820        4.946390
        lon       lat                       geometry
1  80.88125 -53.05634 POLYGON ((80.77247 -44.0823...
2 161.11282 -33.45208 POLYGON ((168.2946 -30.8436...
3  66.71317  58.53052 POLYGON ((73.07079 54.84957...
4 131.93339 -36.47916 POLYGON ((142.0258 -35.2617...
5 162.71433 -23.36775 POLYGON ((164.7341 -16.6514...
6 125.53991 -87.01094 POLYGON ((141.8736 -90.6235...

summary(county_leukemia_128)

       id          leuk_cases             AverageAge  GreenNeighborhoodIndex   CarsPerHH           PopDens         
 Min.   :  1.0   Min.   :29.00   15 and younger:186   Min.   :0.0002333      Min.   :0.001174   Min.   :0.0000152  
 1st Qu.:125.8   1st Qu.:37.00   15-50         :167   1st Qu.:0.2346412      1st Qu.:0.256101   1st Qu.:0.0600216  
 Median :250.5   Median :38.00   50 and older  :147   Median :0.4859838      Median :0.502197   Median :0.1753587  
 Mean   :250.5   Mean   :38.31                        Mean   :0.4862802      Mean   :0.497339   Mean   :0.2427995  
 3rd Qu.:375.2   3rd Qu.:40.00                        3rd Qu.:0.7251946      3rd Qu.:0.749654   3rd Qu.:0.3783928  
 Max.   :500.0   Max.   :48.00                        Max.   :0.9967076      Max.   :0.999106   Max.   :0.9666323  
  PctPopGrowth     Dist2FracLoc_1   Dist2FracLoc_2   Dist2FracLoc_3   Dist2FracLoc_4   Dist2FracLoc_5   Dist2FracLoc_6  
 Min.   :-6.9749   Min.   : 0.000   Min.   : 0.000   Min.   : 0.000   Min.   : 0.000   Min.   : 0.000   Min.   : 0.000  
 1st Qu.:-1.4347   1st Qu.: 3.995   1st Qu.: 3.049   1st Qu.: 2.291   1st Qu.: 3.930   1st Qu.: 2.758   1st Qu.: 2.505  
 Median :-0.2529   Median : 5.060   Median : 5.270   Median : 5.084   Median : 5.023   Median : 4.823   Median : 5.226  
 Mean   :-0.4161   Mean   : 5.032   Mean   : 5.282   Mean   : 4.888   Mean   : 5.003   Mean   : 4.920   Mean   : 5.111  
 3rd Qu.: 0.6439   3rd Qu.: 6.272   3rd Qu.: 7.728   3rd Qu.: 7.321   3rd Qu.: 6.192   3rd Qu.: 7.114   3rd Qu.: 7.721  
 Max.   : 4.9987   Max.   :10.000   Max.   :10.000   Max.   :10.000   Max.   :10.000   Max.   :10.000   Max.   :10.000  
 Dist2FracLoc_7   Dist2FracLoc_8   Dist2FracLoc_9   Dist2FracLoc_10       lon                lat                   geometry  
 Min.   : 0.000   Min.   : 0.000   Min.   : 0.000   Min.   : 0.000   Min.   :-179.340   Min.   :-89.535   POLYGON      :500  
 1st Qu.: 3.331   1st Qu.: 3.615   1st Qu.: 2.506   1st Qu.: 3.835   1st Qu.: -83.564   1st Qu.:-48.626   epsg:4326    :  0  
 Median : 4.990   Median : 5.091   Median : 5.026   Median : 5.204   Median :  11.977   Median : -2.612   +proj=long...:  0  
 Mean   : 5.160   Mean   : 4.995   Mean   : 5.006   Mean   : 5.138   Mean   :   5.288   Mean   : -1.318                      
 3rd Qu.: 7.169   3rd Qu.: 6.306   3rd Qu.: 7.465   3rd Qu.: 6.392   3rd Qu.:  94.009   3rd Qu.: 45.897                      
 Max.   :10.000   Max.   :10.000   Max.   :10.000   Max.   :10.000   Max.   : 179.861   Max.   : 88.581                      

min_value <- min(county_leukemia_128$leuk_cases)
max_value <- max(county_leukemia_128$leuk_cases)
mean_value <- mean(county_leukemia_128$leuk_cases)

boxplot(county_leukemia_128$leuk_cases, 
        main = "Leukemia cases", 
        ylab = "Cases", 
        col = "lightblue",
        ylim = c(20,60))

text(x = 1.2, y = min_value - 5, labels = paste("Min:", round(min_value, 2)), pos = 4, col = "red")

text(x = 1.2, y = max_value + 5, labels = paste("Max:", round(max_value, 2)), pos = 4, col = "red")
text(x = 1.2, y = mean_value, labels = paste("Mean:", round(mean_value, 2)), pos = 4, col = "blue")

points(x = 1, y = mean_value, col = "blue", pch = 18, cex = 1.5)  # pch = 19 for a filled circle

# Create multiple boxplots for the selected variables
boxplot(list(county_leukemia_128$GreenNeighborhoodIndex, 
             county_leukemia_128$CarsPerHH, 
             county_leukemia_128$PopDens),
        names = c("Green Index", 
                  "Cars/Household", 
                  "Population Density"),
        main = "Comparing Green Neighborhood Index, 
        Cars per Household, and Population Density",
        col = c("lightgreen", "lightpink", "lightyellow"),
        ylab = "Values")

# Add the mean values on top of each boxplot
means <- c(mean(county_leukemia_128$GreenNeighborhoodIndex, na.rm = TRUE),
           mean(county_leukemia_128$CarsPerHH, na.rm = TRUE),
           mean(county_leukemia_128$PopDens, na.rm = TRUE))

# Loop through each variable and add the mean as a label
for (i in 1:3) {
  text(x = i, y = means[i], labels = paste("Mean:", round(means[i], 2)), 
       pos = 3, col = "black")
  
  # Add mean dot on the boxplot
  points(x = i, y = means[i], col = "black", pch = 18, cex = 1.5)  # pch = 19 for a solid circle
}

# Create a list of the selected columns
data_for_boxplots <- county_leukemia_128[, c("Dist2FracLoc_1", "Dist2FracLoc_2", "Dist2FracLoc_3", 
                                             "Dist2FracLoc_4", "Dist2FracLoc_5", "Dist2FracLoc_6", 
                                             "Dist2FracLoc_7", "Dist2FracLoc_8", "Dist2FracLoc_9", 
                                             "Dist2FracLoc_10")]

# Create the boxplots for each of the columns
boxplot(data_for_boxplots, 
        main = "Boxplots of Distances to Fracking Sites", 
        xlab = "Distance Category", 
        ylab = "Distance", 
        col = rainbow(10),  # Color the boxes differently
        names = paste("Dst", 1:10, sep = ""),  # Label each boxplot with its corresponding column name
        outline = TRUE)  # Show outliers

Error in sort.int(x, na.last = na.last, decreasing = decreasing, ...) : 
  'x' must be atomic

min_value <- min(county_leukemia_128$PctPopGrowth)
max_value <- max(county_leukemia_128$PctPopGrowth)
mean_value <- mean(county_leukemia_128$PctPopGrowth)

boxplot(county_leukemia_128$PctPopGrowth, 
        main = "PctPopGrowth", 
        ylab = "Values", 
        col = "lightblue",
        ylim = c(-10,10))

text(x = 1.2, y = min_value, labels = paste("Min:", round(min_value, 2)), pos = 4, col = "red")

text(x = 1.2, y = max_value, labels = paste("Max:", round(max_value, 2)), pos = 4, col = "red")
text(x = 1.2, y = mean_value, labels = paste("Mean:", round(mean_value, 2)), pos = 4, col = "blue")

points(x = 1, y = mean_value, col = "blue", pch = 18, cex = 1.5)  # pch = 19 for a filled circle

library(ggplot2)

# Create a scatter plot of longitude vs. latitude
ggplot(county_leukemia_128, aes(x = lon, y = lat)) +
  geom_point(color = "blue", size = 2) +
  labs(title = "Locations of Fracking Sites", x = "Longitude", y = "Latitude") +
  theme_minimal()

Understand “Fracking_location_128”

head(fracking_locations_128)

summary(fracking_locations_128)

 geometry.Length  geometry.Class  geometry.Mode
 2        XY       numeric                     
 2        XY       numeric                     
 2        XY       numeric                     
 2        XY       numeric                     
 2        XY       numeric                     
 2        XY       numeric                     
 2        XY       numeric                     
 2        XY       numeric                     
 2        XY       numeric                     
 2        XY       numeric                     
    FracLoc     
 Min.   : 1.00  
 1st Qu.: 3.25  
 Median : 5.50  
 Mean   : 5.50  
 3rd Qu.: 7.75  
 Max.   :10.00  
                
                
                
                

library(sf)

Linking to GEOS 3.12.2, GDAL 3.9.3, PROJ 9.4.1; sf_use_s2()
is TRUE

library(ggplot2)

# Assuming you have an 'sf' object
fracking_locations <- st_as_sf(fracking_locations_128, coords = c("longitude", "latitude"), crs = 4326)

ggplot(data = fracking_locations) +
  geom_sf() +  # Creates a map plot
  ggtitle("Fracking Locations") +
  theme_minimal()

Experiment Spatial 1

# Assuming you have an 'sf' object for spatial data
library(ggplot2)
library(sf)

# Convert fracking location data into sf object if it's not already
fracking_locations <- st_as_sf(fracking_locations_128, coords = c("longitude", "latitude"), crs = 4326)

# Create a map using ggplot with sf and add scatter plot points
ggplot(data = fracking_locations) +
  geom_sf() +  # Plot the geometry of locations (using 'sf' object)
  ggtitle("Fracking Locations") +
  theme_minimal() +
  labs(x = "Longitude", y = "Latitude")

ggplot(county_leukemia_128, aes(x = lon, y = lat)) +
  geom_point(color = "blue", size = 2) +
  labs(title = "Locations of Fracking Sites", x = "Longitude", y = "Latitude") +
  theme_minimal()

plot(st_geometry(county_leukemia_128), main = "Fracking Locations")

Warning: bounding box has potentially an invalid value range for longlat data

ggplot() +
  # Plot the geometry of fracking locations using 'geom_sf'
  geom_sf(data = fracking_locations, fill = "lightgreen", color = "red", alpha = 1, size = 2) + 
  # Plot the scatter points for the fracking sites
  geom_point(data = county_leukemia_128, aes(x = lon, y = lat), color = "blue", size = 1, alpha = 0.5) + 
  ggtitle("Countries Location and Fracking Sites") +
  labs(x = "Longitude", y = "Latitude") +
  theme_minimal()

# Convert fracking location data into sf object if it's not already
county_leukemia_sf <- st_as_sf(county_leukemia_128, coords = c("lon", "lat"), crs = 4326)

# Create a map using ggplot with sf and add scatter plot points
ggplot(data = county_leukemia_sf) +
  geom_sf() +  # Plot the geometry of locations (using 'sf' object)
  geom_point(aes(x = lon, y = lat), color = "blue", size = 2) +  # Scatter plot for points
  ggtitle("Fracking Locations") +
  theme_minimal() +
  labs(x = "Longitude", y = "Latitude")

Finally get map show Boarders, countries, and Fracking Sites.

county_leukemia_sf <- st_as_sf(county_leukemia_128, coords = c("lon", "lat"), crs = 4326)

ggplot(county_leukemia_sf) +
  geom_sf() +
  # Plot the geometry of fracking locations using 'geom_sf'
  geom_sf(data = fracking_locations, fill = "lightgreen", color = "red", alpha = 1, size = 2) + 
  # Plot the scatter points for the fracking sites
  geom_point(data = county_leukemia_128, aes(x = lon, y = lat), color = "blue", size = 1, alpha = 0.5) +
  ggtitle("Countries Location and Fracking Sites") +
  labs(x = "Longitude", y = "Latitude") +
  theme_minimal()

NA
NA

# Extract coordinates for fracking locations
fracking_coords <- st_coordinates(fracking_locations)

# Add the coordinates as columns to the fracking_locations object
fracking_locations$lon <- fracking_coords[, 1]
fracking_locations$lat <- fracking_coords[, 2]

ggplot(county_leukemia_sf) +
  geom_sf() +
  # Plot the geometry of fracking locations using 'geom_sf'
  geom_sf(data = fracking_locations, fill = "lightgreen", color = "red", alpha = 1, size = 2) + 
  # Plot the scatter points for the fracking sites
  geom_point(data = county_leukemia_128, aes(x = lon, y = lat), color = "blue", size = 1, alpha = 0.5) +
  geom_text(data = fracking_locations, aes(x = lon, y = lat, label = as.character(FracLoc)), 
            color = "black", size = 3, vjust = 0, hjust = 0) +  # Adjust position for clarity
  ggtitle("Countries Location and Fracking Sites") +
  labs(x = "Longitude", y = "Latitude") +
  theme_minimal()

clean Data

ls()

 [1] "all_objects"            "coords"                
 [3] "county_leukemia_128"    "county_leukemia_sf"    
 [5] "county_leukemia_sp"     "county_locations"      
 [7] "county_ppp"             "data_for_boxplots"     
 [9] "dist_test"              "envcase"               
[11] "file_name"              "fracking_coords"       
[13] "fracking_locations"     "fracking_locations_128"
[15] "gam_model"              "glm_sqrt_model"        
[17] "i"                      "invalid_geometries"    
[19] "K"                      "max_value"             
[21] "mean_value"             "means"                 
[23] "min_value"              "missing_data"          
[25] "nndistances"            "objects_to_delete"     
[27] "poisson_model"          "r"                     
[29] "seed_number"            "selected_cars"         

There is no missing Data

# Identify columns with missing values
missing_data <- sapply(county_leukemia_128, function(x) sum(is.na(x)))
print(missing_data)

                    id             leuk_cases 
                     0                      0 
            AverageAge GreenNeighborhoodIndex 
                     0                      0 
             CarsPerHH                PopDens 
                     0                      0 
          PctPopGrowth         Dist2FracLoc_1 
                     0                      0 
        Dist2FracLoc_2         Dist2FracLoc_3 
                     0                      0 
        Dist2FracLoc_4         Dist2FracLoc_5 
                     0                      0 
        Dist2FracLoc_6         Dist2FracLoc_7 
                     0                      0 
        Dist2FracLoc_8         Dist2FracLoc_9 
                     0                      0 
       Dist2FracLoc_10                    lon 
                     0                      0 
                   lat               geometry 
                     0                      0 

missing_data <- sapply(fracking_locations_128, function(x) sum(is.na(x)))
print(missing_data)

geometry  FracLoc 
       0        0 

Data Type Check : Format & Factor

# Convert 'AverageAge' to a factor with levels
county_leukemia_128$AverageAge <- factor(county_leukemia_128$AverageAge, 
                                          levels = c("15 and younger", "15-50", "50 and older"))

# Verify the conversion
summary(county_leukemia_128$AverageAge)

15 and younger          15-50   50 and older 
           186            167            147 

county_leukemia_128$leuk_cases <- as.numeric(county_leukemia_128$leuk_cases)

dist_test <- county_leukemia_128$Dist2FracLoc_1

I Couldn’t come up with any reason why I find any Outliers unreasonable, so no elimination or cap.

Feature Selection

SP-X-Y

SP-X-Y

Spatial : lon, lat, geometry, Dist2FracLoc

Outcome (Y) : leuk_cases

Covariates (X) : PopDens, GreenNeighborhoodIndex, AverageAge, PctPopGrowth, CarsPerHH

ggplot(county_leukemia_128) +
  geom_sf(aes(fill = leuk_cases)) +
  theme_minimal()

Scale filled Gradient.

# Extract coordinates for fracking locations
fracking_coords <- st_coordinates(fracking_locations)

# Add the coordinates as columns to the fracking_locations object
fracking_locations$lon <- fracking_coords[, 1]
fracking_locations$lat <- fracking_coords[, 2]

ggplot() +
  geom_sf(data = county_leukemia_sf, aes(fill = county_leukemia_sf$leuk_cases), color = "black") +
  # Plot the geometry of fracking locations using 'geom_sf'
  geom_sf(data = fracking_locations, fill = "lightgreen", color = "red", alpha = 1, size = 2) + 
  # Plot the scatter points for the fracking sites
  geom_point(data = county_leukemia_128, aes(x = lon, y = lat), color = "blue", size = 1, alpha = 0.5) +
  geom_text(data = fracking_locations, aes(x = lon, y = lat, label = as.character(FracLoc)), 
            color = "black", size = 3, vjust = 0, hjust = 0) +  # Adjust position for clarity
  ggtitle("Countries Location and Fracking Sites") +
  labs(x = "Longitude", y = "Latitude") +
  # Use a gradient scale to fill with darker colors for higher 'leuk_cases'
  scale_fill_gradient(low = "lightblue", high = "darkblue") +  # Light color for lower values, dark for higher
  theme_minimal() +
  theme(legend.position = "bottom")

Modeling

Poisson Regression Test

Try with general co variate. — PctPopGrowth looks Correlate.

# Load necessary library
library(MASS)  # For Poisson regression
library(ggplot2)

# Fit a Poisson regression model
poisson_model <- glm(leuk_cases ~ PopDens + AverageAge + GreenNeighborhoodIndex + CarsPerHH + PctPopGrowth,
                     family = poisson(link = "log"),  # log link is default in Poisson regression
                     data = county_leukemia_128)

# Check model summary
summary(poisson_model)

Call:
glm(formula = leuk_cases ~ PopDens + AverageAge + GreenNeighborhoodIndex + 
    CarsPerHH + PctPopGrowth, family = poisson(link = "log"), 
    data = county_leukemia_128)

Coefficients:
                       Estimate Std. Error z value Pr(>|z|)
(Intercept)             3.64119    0.02099 173.436  < 2e-16
PopDens                -0.01640    0.04485  -0.366  0.71464
AverageAge15-50         0.02384    0.01731   1.378  0.16825
AverageAge50 and older  0.02493    0.01791   1.391  0.16410
GreenNeighborhoodIndex -0.02388    0.02533  -0.943  0.34586
CarsPerHH               0.02000    0.03399   0.588  0.55622
PctPopGrowth            0.01296    0.00407   3.183  0.00146
                          
(Intercept)            ***
PopDens                   
AverageAge15-50           
AverageAge50 and older    
GreenNeighborhoodIndex    
CarsPerHH                 
PctPopGrowth           ** 
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 109.788  on 499  degrees of freedom
Residual deviance:  95.743  on 493  degrees of freedom
AIC: 2852.3

Number of Fisher Scoring iterations: 3

Try with Fraction Location — Non of these look correlate.

# Fit a Poisson regression model
poisson_model <- glm(leuk_cases ~ Dist2FracLoc_1 + Dist2FracLoc_2 + Dist2FracLoc_3 + Dist2FracLoc_4 + Dist2FracLoc_5 + Dist2FracLoc_6 + Dist2FracLoc_7 + Dist2FracLoc_8 + Dist2FracLoc_9 + Dist2FracLoc_10,
                     family = poisson(link = "log"),  # log link is default in Poisson regression
                     data = county_leukemia_128)

# Check model summary
summary(poisson_model)

Call:
glm(formula = leuk_cases ~ Dist2FracLoc_1 + Dist2FracLoc_2 + 
    Dist2FracLoc_3 + Dist2FracLoc_4 + Dist2FracLoc_5 + Dist2FracLoc_6 + 
    Dist2FracLoc_7 + Dist2FracLoc_8 + Dist2FracLoc_9 + Dist2FracLoc_10, 
    family = poisson(link = "log"), data = county_leukemia_128)

Coefficients:
                Estimate Std. Error z value Pr(>|z|)    
(Intercept)      2.60448    0.44367   5.870 4.35e-09 ***
Dist2FracLoc_1   0.02288    0.04129   0.554    0.579    
Dist2FracLoc_2   0.02291    0.03875   0.591    0.554    
Dist2FracLoc_3   0.01574    0.04653   0.338    0.735    
Dist2FracLoc_4   0.05249    0.06724   0.781    0.435    
Dist2FracLoc_5  -0.01305    0.02534  -0.515    0.607    
Dist2FracLoc_6   0.00448    0.03226   0.139    0.890    
Dist2FracLoc_7  -0.01735    0.05560  -0.312    0.755    
Dist2FracLoc_8   0.14237    0.09916   1.436    0.151    
Dist2FracLoc_9   0.05339    0.04030   1.325    0.185    
Dist2FracLoc_10 -0.07460    0.07589  -0.983    0.326    
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 109.788  on 499  degrees of freedom
Residual deviance:  58.585  on 489  degrees of freedom
AIC: 2823.2

Number of Fisher Scoring iterations: 3

GLM_Sqrt Model – Found Age, Neighborhood Index, and Pop Growth Correlated.

# Fit a GLM with square root transformation
glm_sqrt_model <- glm(sqrt(leuk_cases) ~ PopDens + AverageAge + GreenNeighborhoodIndex + CarsPerHH + PctPopGrowth,
                       family = gaussian(link = "identity"),  # Identity link for the square root transformation
                       data = county_leukemia_128)

# Check model summary
summary(glm_sqrt_model)

Call:
glm(formula = sqrt(leuk_cases) ~ PopDens + AverageAge + GreenNeighborhoodIndex + 
    CarsPerHH + PctPopGrowth, family = gaussian(link = "identity"), 
    data = county_leukemia_128)

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)
(Intercept)             6.172825   0.028703 215.061  < 2e-16
PopDens                -0.049021   0.061499  -0.797  0.42578
AverageAge15-50         0.073641   0.023673   3.111  0.00197
AverageAge50 and older  0.076896   0.024512   3.137  0.00181
GreenNeighborhoodIndex -0.074678   0.034684  -2.153  0.03180
CarsPerHH               0.060235   0.046622   1.292  0.19697
PctPopGrowth            0.039404   0.005555   7.094 4.57e-12
                          
(Intercept)            ***
PopDens                   
AverageAge15-50        ** 
AverageAge50 and older ** 
GreenNeighborhoodIndex *  
CarsPerHH                 
PctPopGrowth           ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 0.04896861)

    Null deviance: 27.572  on 499  degrees of freedom
Residual deviance: 24.142  on 493  degrees of freedom
AIC: -80.399

Number of Fisher Scoring iterations: 2

summary(county_leukemia_128$AverageAge)

15 and younger          15-50   50 and older 
           186            167            147 

GLM With Spatial – somehow related to 4, 8, 9, 10

# Fit a GLM with square root transformation
glm_sqrt_model <- glm(sqrt(leuk_cases) ~ Dist2FracLoc_1 + Dist2FracLoc_2 + 
    Dist2FracLoc_3 + Dist2FracLoc_4 + Dist2FracLoc_5 + Dist2FracLoc_6 + 
    Dist2FracLoc_7 + Dist2FracLoc_8 + Dist2FracLoc_9 + Dist2FracLoc_10,
    family = gaussian(link = "identity"),  # Identity link for the square root transformation
                       data = county_leukemia_128)

# Check model summary
summary(glm_sqrt_model)

Call:
glm(formula = sqrt(leuk_cases) ~ Dist2FracLoc_1 + Dist2FracLoc_2 + 
    Dist2FracLoc_3 + Dist2FracLoc_4 + Dist2FracLoc_5 + Dist2FracLoc_6 + 
    Dist2FracLoc_7 + Dist2FracLoc_8 + Dist2FracLoc_9 + Dist2FracLoc_10, 
    family = gaussian(link = "identity"), data = county_leukemia_128)

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)      2.89931    0.47536   6.099 2.17e-09 ***
Dist2FracLoc_1   0.06941    0.04420   1.570 0.116986    
Dist2FracLoc_2   0.07103    0.04148   1.712 0.087478 .  
Dist2FracLoc_3   0.04999    0.04986   1.003 0.316585    
Dist2FracLoc_4   0.17061    0.07200   2.370 0.018195 *  
Dist2FracLoc_5  -0.04118    0.02715  -1.517 0.130035    
Dist2FracLoc_6   0.01328    0.03453   0.385 0.700633    
Dist2FracLoc_7  -0.05163    0.05930  -0.871 0.384442    
Dist2FracLoc_8   0.45581    0.10778   4.229 2.80e-05 ***
Dist2FracLoc_9   0.16531    0.04296   3.848 0.000135 ***
Dist2FracLoc_10 -0.24125    0.08236  -2.929 0.003558 ** 
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 0.02989407)

    Null deviance: 27.572  on 499  degrees of freedom
Residual deviance: 14.618  on 489  degrees of freedom
AIC: -323.23

Number of Fisher Scoring iterations: 2

Try GAM – only related to PctPopGrowth

# Load the mgcv package for GAM
library(mgcv)

Loading required package: nlme
This is mgcv 1.9-1. For overview type 'help("mgcv-package")'.

# Fit a GAM model
gam_model <- gam(leuk_cases ~ s(PopDens) + (AverageAge) + s(GreenNeighborhoodIndex) + s(CarsPerHH) + s(PctPopGrowth),
                 family = poisson(link = "log"),  # Poisson regression for count data
                 data = county_leukemia_128)

# Check model summary
summary(gam_model)

Family: poisson 
Link function: log 

Formula:
leuk_cases ~ s(PopDens) + (AverageAge) + s(GreenNeighborhoodIndex) + 
    s(CarsPerHH) + s(PctPopGrowth)

Parametric coefficients:
                       Estimate Std. Error z value Pr(>|z|)
(Intercept)             3.63016    0.01196 303.469   <2e-16
AverageAge15-50         0.02384    0.01730   1.378    0.168
AverageAge50 and older  0.02493    0.01791   1.391    0.164
                          
(Intercept)            ***
AverageAge15-50           
AverageAge50 and older    
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                          edf Ref.df Chi.sq p-value   
s(PopDens)                  1      1  0.134 0.71465   
s(GreenNeighborhoodIndex)   1      1  0.889 0.34587   
s(CarsPerHH)                1      1  0.346 0.55624   
s(PctPopGrowth)             1      1 10.133 0.00146 **
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.119   Deviance explained = 12.8%
UBRE = -0.78051  Scale est. = 1         n = 500

Retry with Spatial information – show no sign of relationship

# Fit a GAM model
gam_model <- gam(leuk_cases ~ s(Dist2FracLoc_1) + s(Dist2FracLoc_2) + 
    s(Dist2FracLoc_3) + s(Dist2FracLoc_4) + s(Dist2FracLoc_5) + s(Dist2FracLoc_6) + 
    s(Dist2FracLoc_7) + s(Dist2FracLoc_8) + s(Dist2FracLoc_9) + s(Dist2FracLoc_10),,
                 family = poisson(link = "log"),  # Poisson regression for count data
                 data = county_leukemia_128)

# Check model summary
summary(gam_model)

Family: poisson 
Link function: log 

Formula:
leuk_cases ~ s(Dist2FracLoc_1) + s(Dist2FracLoc_2) + s(Dist2FracLoc_3) + 
    s(Dist2FracLoc_4) + s(Dist2FracLoc_5) + s(Dist2FracLoc_6) + 
    s(Dist2FracLoc_7) + s(Dist2FracLoc_8) + s(Dist2FracLoc_9) + 
    s(Dist2FracLoc_10)

Parametric coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) 3.644458   0.007235   503.7   <2e-16 ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                     edf Ref.df Chi.sq p-value
s(Dist2FracLoc_1)  1.000  1.000  0.362   0.547
s(Dist2FracLoc_2)  1.000  1.000  0.420   0.517
s(Dist2FracLoc_3)  1.000  1.000  0.094   0.759
s(Dist2FracLoc_4)  1.000  1.000  0.547   0.460
s(Dist2FracLoc_5)  1.000  1.000  0.244   0.621
s(Dist2FracLoc_6)  1.000  1.000  0.023   0.879
s(Dist2FracLoc_7)  1.003  1.005  0.112   0.739
s(Dist2FracLoc_8)  1.000  1.000  1.977   0.160
s(Dist2FracLoc_9)  1.663  2.007  3.255   0.195
s(Dist2FracLoc_10) 1.000  1.000  0.889   0.346

R-sq.(adj) =  0.478   Deviance explained = 48.8%
UBRE = -0.84087  Scale est. = 1         n = 500

Modeling SPPA

SPPA #1

# install.packages("spatstat")
library(spatstat)

Loading required package: spatstat.data
Loading required package: spatstat.univar
spatstat.univar 3.1-1
Loading required package: spatstat.geom
spatstat.geom 3.3-4

Attaching package: ‘spatstat.geom’

The following object is masked from ‘package:MASS’:

    area

Loading required package: spatstat.random
spatstat.random 3.3-2
Loading required package: spatstat.explore
spatstat.explore 3.3-3
Loading required package: spatstat.model
Loading required package: rpart
spatstat.model 3.3-3
Loading required package: spatstat.linnet
spatstat.linnet 3.2-3

spatstat 3.3-0 
For an introduction to spatstat, type ‘beginner’ 

coords <- st_coordinates(county_leukemia_sf)


# Create a point pattern object (ppp) from the coordinates
county_ppp <- ppp(x = coords[, 1], y = coords[, 2], 
                  window = owin(c(min(coords[, 1]), max(coords[, 1])), 
                                c(min(coords[, 2]), max(coords[, 2]))))

Warning: data contain duplicated points

set.seed(120109)
r <- seq(0, .5, by = 0.001)# number and values of distances r
envcase <- envelope(county_ppp, fun=Gest, r=r, nrank=2, nsim=99)

Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45,
46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75,
76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90,
91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.

Ripley’s K — Above expected — Cluster

# Calculate Ripley’s K function
K <- Kest(county_ppp)

number of data points exceeds 3000 - computing border correction estimate only

# Plot the results
plot(K)

Nearest Neighbor Analysis – ??

# Calculate nearest neighbor distances
nndistances <- nndist(county_ppp)

# Plot the nearest neighbor distances
hist(nndistances, main = "Nearest Neighbor Distances", xlab = "Distance")

Try Sven Code

library(sp)
library(spatstat)

county_leukemia_sp <- as(county_leukemia_sf, "Spatial")
county_leukemia_sp_young <- county_leukemia_sp[county_leukemia_sp$AverageAge == "15 and younger"]

Error in `[.data.frame`(x@data, i, j, ..., drop = FALSE) : 
  undefined columns selected

Try GG Pair

library(GGally)

# Convert AverageAge to a factor if it is not already
# county_leukemia_sf$AverageAge <- factor(county_leukemia_sf$AverageAge,
# levels = c("15 and younger", "15-50", "50 and older"))


selected_vars <- county_leukemia_sf[, c("leuk_cases", "GreenNeighborhoodIndex", "CarsPerHH") ]

selected_vars <- st_drop_geometry(selected_vars)

ggpairs(selected_vars, title = "Scatterplot Matrix of Selected Variables")

Can we GG-pair with Location?

library(GGally)

# Convert AverageAge to a factor if it is not already
# county_leukemia_sf$AverageAge <- factor(county_leukemia_sf$AverageAge,
# levels = c("15 and younger", "15-50", "50 and older"))


selected_location_vars <- county_leukemia_sf[, c("leuk_cases", "Dist2FracLoc_1", "Dist2FracLoc_2") ]

selected_location_vars <- st_drop_geometry(selected_location_vars)

ggpairs(selected_location_vars, title = "Scatterplot Matrix of Selected Variables")