SJV CRC PCA

packages <- c('readxl','dplyr','ggplot2')
missing_packages <- setdiff(packages, rownames(installed.packages()))
if (length(missing_packages)) install.packages(missing_packages)
library(readxl)
library(dplyr)
library(ggplot2)

index <-  read_excel('C:/Users/noahs/Desktop/Ch1/SJV/Index.xlsx')

# Keep tract IDs as 11-character FIPS codes
# because Excel may store these as numbers
index <-  read_excel('C:/Users/noahs/Desktop/Ch1/SJV/Index.xlsx')
index$tract <- sprintf('%011.0f', as.numeric(index$tract))

themes <- list(
  LU = c('ImpairedWaterBodies','SmATS','avg_impervious','tree_canopy'),
  HV = c('asthma','Low_Birth_Weight','Cardiovascular'),
  HE = c('ozone','PM2_5','Diesel_PM','Pesticides','traffic','LST_Max','LST_Mean'),
  SV = c('Linguistic_Isol','Poverty','HousBurd','children','elderly',
         'pop_density','Total_outdoor','perc_no_vehicle_access','energy_burden')
)
# Combine all indicators into one list
vars <- unlist(themes, use.names = FALSE)

# Check for missing indicator columns
missing_vars <- setdiff(vars, names(index))

if (length(missing_vars) > 0) {
  stop(paste("Missing columns:",
             paste(missing_vars, collapse = ", ")))
}

# Convert indicators to numeric values
for (v in vars) {
  index[[v]] <- as.numeric(index[[v]])
}

# Check for missing values after conversion
colSums(is.na(index[vars]))

   ImpairedWaterBodies                  SmATS         avg_impervious            tree_canopy                 asthma       Low_Birth_Weight         Cardiovascular                  ozone                  PM2_5 
                     0                      0                      0                      0                      0                      1                      0                      0                      0 
             Diesel_PM             Pesticides                traffic                LST_Max               LST_Mean        Linguistic_Isol                Poverty               HousBurd               children 
                     0                      0                      1                      0                      1                      1                      1                      1                      1 
               elderly            pop_density          Total_outdoor perc_no_vehicle_access          energy_burden 
                     0                      0                      0                      7                      0 

# Check missing values for all indicators
missing_data <- colSums(is.na(index[vars]))

print(missing_data)

   ImpairedWaterBodies                  SmATS         avg_impervious            tree_canopy                 asthma       Low_Birth_Weight         Cardiovascular                  ozone                  PM2_5 
                     0                      0                      0                      0                      0                      1                      0                      0                      0 
             Diesel_PM             Pesticides                traffic                LST_Max               LST_Mean        Linguistic_Isol                Poverty               HousBurd               children 
                     0                      0                      1                      0                      1                      1                      1                      1                      1 
               elderly            pop_density          Total_outdoor perc_no_vehicle_access          energy_burden 
                     0                      0                      0                      7                      0 

# Identify census tracts with missing values
incomplete_tracts <- index[
  !complete.cases(index[vars]),
  c("tract", vars)
]

print(incomplete_tracts)

# A tibble: 9 × 24
  tract  ImpairedWaterBodies SmATS avg_impervious tree_canopy asthma Low_Birth_Weight Cardiovascular  ozone PM2_5 Diesel_PM Pesticides traffic LST_Max LST_Mean Linguistic_Isol Poverty HousBurd children elderly
  <chr>                <dbl> <dbl>          <dbl>       <dbl>  <dbl>            <dbl>          <dbl>  <dbl> <dbl>     <dbl>      <dbl>   <dbl>   <dbl>    <dbl>           <dbl>   <dbl>    <dbl>    <dbl>   <dbl>
1 06077…                   6   0.9          47.5        1.43    69.7             8.3            21.5 0.0441  8.65    0.0181      7.55     NA      45.6     42.6           19.5    19.9     14.8     12.4     7.32
2 06077…                   9  63.8           1.30       2.10    52.8             6.44           15.1 0.0441  7.67    0.0389    147.      171.     56.1     NA              1.64   16.3      8.75    10.3    19.1 
3 06077…                   4 153.           30.6        1.98    47.2             6.24           20.5 0.0454  8.47    0.0978    112.      791.     58.0     43.1           10.0    24.9     21.6     17.3     5.69
4 06077…                   0  70.0          40.7        3.99    44.8             6              20.1 0.0459  8.66    0.0880      0.140   266.     50.3     43.1            5.31   22.1     16        8.14   11.9 
5 06077…                   0 120.           54.3        5.41    40.4             5.11           20.1 0.0463  8.64    0.226       2.33    553.     49.7     45.6            5.95   36.0     20.1     11.8    13.8 
6 06077…                   0  15            60.2        2.45    40.5             7.86           17.7 0.0440  5.54    0.0175      0       125.     49.9     45.9            6.30   12.8     17.4     12.9    13.5 
7 06077…                   0  15.3          35.5        2.36    49.8             7.02           18.4 0.0443  6.25    0.0492      0.463   235.     52.3     47.2            1.88   24.9     10        7.25   11.2 
8 06077…                   3   0.2          48.1        4.99    28.1             4.95           14.4 0.0436  8.07    0.0238     32.6     394.     46.7     41.8            2.00    8.12     3.40     9.96    5.96
9 06077…                   3   0.2          29.0        0.502   41.5            NA              15.6 0.0453  7.50    0.0597     62.6      18.3    46.9     42.7           NA      NA       NA       NA      13.1 
# ℹ 4 more variables: pop_density <dbl>, Total_outdoor <dbl>, perc_no_vehicle_access <dbl>, energy_burden <dbl>

# Min-max normalization
normalize <- function(x) {
  min_x <- min(x, na.rm = TRUE)
  max_x <- max(x, na.rm = TRUE)

  if (max_x == min_x) {
    return(rep(NA_real_, length(x)))
  }

  (x - min_x) / (max_x - min_x)
}

# Create normalized dataset
normalized <- index

# Normalize all 23 indicators
for (v in vars) {
  normalized[[v]] <- normalize(index[[v]])
}

# Reverse tree canopy
# Lower canopy = higher vulnerability
normalized$tree_canopy <- 1 - normalized$tree_canopy

# Check normalized values
summary(normalized[vars])

 ImpairedWaterBodies     SmATS         avg_impervious    tree_canopy         asthma       Low_Birth_Weight Cardiovascular       ozone            PM2_5          Diesel_PM         Pesticides       
 Min.   :0.0000      Min.   :0.00000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.00000   Min.   :0.0000000  
 1st Qu.:0.0000      1st Qu.:0.04781   1st Qu.:0.4158   1st Qu.:0.4808   1st Qu.:0.1479   1st Qu.:0.4322   1st Qu.:0.3413   1st Qu.:0.3553   1st Qu.:0.5368   1st Qu.:0.02400   1st Qu.:0.0000000  
 Median :0.1579      Median :0.12426   Median :0.6373   Median :0.6952   Median :0.2619   Median :0.5491   Median :0.5998   Median :0.4824   Median :0.6564   Median :0.04408   Median :0.0005198  
 Mean   :0.2093      Mean   :0.19339   Mean   :0.5402   Mean   :0.6454   Mean   :0.2924   Mean   :0.5577   Mean   :0.5549   Mean   :0.4631   Mean   :0.6571   Mean   :0.09086   Mean   :0.0518297  
 3rd Qu.:0.3553      3rd Qu.:0.24496   3rd Qu.:0.6973   3rd Qu.:0.8310   3rd Qu.:0.4135   3rd Qu.:0.6847   3rd Qu.:0.7853   3rd Qu.:0.5299   3rd Qu.:0.8104   3rd Qu.:0.08704   3rd Qu.:0.0265091  
 Max.   :1.0000      Max.   :1.00000   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1.00000   Max.   :1.0000000  
                                                                                          NA's   :1                                                                                                
    traffic          LST_Max          LST_Mean      Linguistic_Isol      Poverty          HousBurd         children         elderly        pop_density     Total_outdoor    perc_no_vehicle_access
 Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.00000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.00000       
 1st Qu.:0.1365   1st Qu.:0.1667   1st Qu.:0.5493   1st Qu.:0.07489   1st Qu.:0.2059   1st Qu.:0.1891   1st Qu.:0.3056   1st Qu.:0.1362   1st Qu.:0.1705   1st Qu.:0.2337   1st Qu.:0.04503       
 Median :0.2492   Median :0.2236   Median :0.6896   Median :0.17736   Median :0.3425   Median :0.2892   Median :0.4308   Median :0.2674   Median :0.3327   Median :0.3240   Median :0.08628       
 Mean   :0.2957   Mean   :0.2632   Mean   :0.6324   Mean   :0.21637   Mean   :0.3907   Mean   :0.3134   Mean   :0.4345   Mean   :0.2844   Mean   :0.3248   Mean   :0.3413   Mean   :0.13101       
 3rd Qu.:0.3995   3rd Qu.:0.3308   3rd Qu.:0.7899   3rd Qu.:0.27681   3rd Qu.:0.5332   3rd Qu.:0.4066   3rd Qu.:0.5646   3rd Qu.:0.3818   3rd Qu.:0.4570   3rd Qu.:0.4169   3rd Qu.:0.17308       
 Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1.00000   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1.00000       
 NA's   :1                         NA's   :1        NA's   :1         NA's   :1        NA's   :1        NA's   :1                                                           NA's   :7             
 energy_burden   
 Min.   :0.0000  
 1st Qu.:0.3333  
 Median :0.5000  
 Mean   :0.4751  
 3rd Qu.:0.5000  
 Max.   :1.0000  
                 

# Land Use & Environmental Features
normalized$LU <- rowMeans(
  normalized[themes$LU],
  na.rm = FALSE
)

# Health Vulnerability
normalized$HV <- rowMeans(
  normalized[themes$HV],
  na.rm = FALSE
)

# Heat & Environmental Exposure
normalized$HE <- rowMeans(
  normalized[themes$HE],
  na.rm = FALSE
)

# Socioeconomic Vulnerability
normalized$SV <- rowMeans(
  normalized[themes$SV],
  na.rm = FALSE
)

# Calculate Opportunity Index
normalized$Opportunity_Index <- rowMeans(
  normalized[c("LU", "HV", "HE", "SV")],
  na.rm = FALSE
)

# View results
summary(normalized$Opportunity_Index)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
 0.2100  0.3193  0.3752  0.3866  0.4451  0.6396       9 

# Create lists for PCA results
pca_results <- list()
pca_loadings <- list()
pca_variance <- list()

# Run PCA for each theme
for (theme in names(themes)) {

  # Select indicators
  data <- index[themes[[theme]]]

  # Remove rows with missing values
  data <- na.omit(data)

  # Standardize and run PCA
  pca <- prcomp(
    data,
    center = TRUE,
    scale. = TRUE
  )

  # Save PCA results
  pca_results[[theme]] <- pca

  # Extract loadings
  pca_loadings[[theme]] <- pca$rotation

  # Calculate explained variance
  variance <- pca$sdev^2 / sum(pca$sdev^2)

  pca_variance[[theme]] <- variance

  # Print results
  cat("\n========================\n")
  cat("PCA Theme:", theme, "\n")
  cat("========================\n")

  print(summary(pca))
  print(pca$rotation)
}

========================
PCA Theme: LU 
========================
Importance of components:
                          PC1    PC2    PC3    PC4
Standard deviation     1.2292 1.0492 0.9534 0.6922
Proportion of Variance 0.3777 0.2752 0.2272 0.1198
Cumulative Proportion  0.3777 0.6530 0.8802 1.0000
                           PC1        PC2        PC3        PC4
ImpairedWaterBodies  0.5886010 -0.4650469 -0.2592655 -0.6083270
SmATS               -0.3402914 -0.6963567 -0.4820635  0.4085387
avg_impervious      -0.6422848  0.2187652 -0.3901179 -0.6224307
tree_canopy          0.3538656  0.5009596 -0.7404081  0.2749808

========================
PCA Theme: HV 
========================
Importance of components:
                          PC1    PC2    PC3
Standard deviation     1.2562 0.9345 0.7406
Proportion of Variance 0.5261 0.2911 0.1828
Cumulative Proportion  0.5261 0.8172 1.0000
                        PC1        PC2        PC3
asthma           -0.6578798  0.1342572 -0.7410595
Low_Birth_Weight -0.5987070  0.5037045  0.6227614
Cardiovascular   -0.4568852 -0.8533796  0.2509963

========================
PCA Theme: HE 
========================
Importance of components:
                          PC1    PC2    PC3    PC4     PC5     PC6     PC7
Standard deviation     1.4550 1.1505 1.0360 0.9797 0.82166 0.71578 0.58194
Proportion of Variance 0.3024 0.1891 0.1533 0.1371 0.09645 0.07319 0.04838
Cumulative Proportion  0.3024 0.4915 0.6449 0.7820 0.87843 0.95162 1.00000
                  PC1         PC2        PC3         PC4        PC5         PC6         PC7
ozone      -0.2787390  0.52456281 -0.2993920 -0.16178113 -0.7056835 -0.11950813  0.13805178
PM2_5       0.4281421  0.39919658  0.2774168 -0.18939281  0.1852795 -0.71189920  0.05811841
Diesel_PM   0.2792124  0.62588852 -0.1996177 -0.05958408  0.3832975  0.56760427  0.13347248
Pesticides -0.5318277  0.31420723  0.2086590 -0.11608873  0.2476442 -0.02470387 -0.70673901
traffic     0.2782188  0.16462861  0.6688910  0.37736956 -0.4395314  0.29204597 -0.16489431
LST_Max    -0.2774607  0.21540812 -0.1175404  0.84488840  0.2090034 -0.24488033  0.21287085
LST_Mean    0.4730841 -0.03885984 -0.5373897  0.25433480 -0.1502244 -0.10428429 -0.62270771

========================
PCA Theme: SV 
========================
Importance of components:
                          PC1    PC2    PC3     PC4     PC5     PC6     PC7     PC8     PC9
Standard deviation     2.0083 1.1023 1.0244 0.88733 0.75612 0.66756 0.62318 0.55096 0.45377
Proportion of Variance 0.4481 0.1350 0.1166 0.08748 0.06352 0.04951 0.04315 0.03373 0.02288
Cumulative Proportion  0.4481 0.5831 0.6997 0.78720 0.85073 0.90024 0.94339 0.97712 1.00000
                              PC1          PC2         PC3         PC4         PC5         PC6        PC7         PC8         PC9
Linguistic_Isol        -0.4149005  0.107479556 -0.20230411 -0.03343097  0.15527883 -0.02821378 -0.3380019 -0.77460524 -0.18733657
Poverty                -0.4417129  0.135626913  0.04462317  0.11720996 -0.20729694  0.07704066 -0.1015468  0.03481641  0.84281339
HousBurd               -0.3670935  0.186538771  0.28056627 -0.13192789 -0.17447356 -0.74035172  0.3636362  0.02485198 -0.15136913
children               -0.2218364 -0.527728099  0.13493684  0.69505130 -0.08484864  0.10828159  0.3359050 -0.15841301 -0.11889655
elderly                 0.2676716  0.592532037  0.20587422  0.13815726 -0.39947177  0.33375232  0.3451782 -0.34973825 -0.05790471
pop_density            -0.2001110 -0.292895328  0.64238418 -0.53287431  0.07853711  0.38867041  0.0880527 -0.11054547 -0.01868365
Total_outdoor          -0.3039878 -0.005944897 -0.59412342 -0.30374637  0.07354248  0.26962389  0.6153776  0.07152041 -0.02003120
perc_no_vehicle_access -0.2862734  0.462671200  0.22007226  0.29841112  0.65062428  0.18668151  0.0242996  0.29501148 -0.14393663
energy_burden          -0.4033915  0.076012865 -0.06083070  0.01701178 -0.54846142  0.25183488 -0.3524175  0.38259397 -0.43897672

# Create explained variance table
variance_table <- data.frame()

for (theme in names(pca_variance)) {

  variance <- pca_variance[[theme]]

  temp <- data.frame(
    Theme = theme,
    Component = paste0("PC", seq_along(variance)),
    Explained_Variance = variance * 100,
    Cumulative_Variance = cumsum(variance) * 100
  )

  variance_table <- rbind(
    variance_table,
    temp
  )
}

print(variance_table)

   Theme Component Explained_Variance Cumulative_Variance
1     LU       PC1          37.775198            37.77520
2     LU       PC2          27.521826            65.29702
3     LU       PC3          22.723343            88.02037
4     LU       PC4          11.979633           100.00000
5     HV       PC1          52.604857            52.60486
6     HV       PC2          29.111320            81.71618
7     HV       PC3          18.283823           100.00000
8     HE       PC1          30.243669            30.24367
9     HE       PC2          18.908299            49.15197
10    HE       PC3          15.333320            64.48529
11    HE       PC4          13.712935            78.19822
12    HE       PC5           9.644661            87.84288
13    HE       PC6           7.319116            95.16200
14    HE       PC7           4.838000           100.00000
15    SV       PC1          44.812051            44.81205
16    SV       PC2          13.499525            58.31158
17    SV       PC3          11.660407            69.97198
18    SV       PC4           8.748383            78.72037
19    SV       PC5           6.352444            85.07281
20    SV       PC6           4.951458            90.02427
21    SV       PC7           4.315062            94.33933
22    SV       PC8           3.372812            97.71214
23    SV       PC9           2.287858           100.00000

# Extract PC1 loadings
loading_table <- data.frame()

for (theme in names(pca_results)) {

  pca <- pca_results[[theme]]

  temp <- data.frame(
    Theme = theme,
    Indicator = rownames(pca$rotation),
    PC1_Loading = pca$rotation[, 1]
  )

  loading_table <- rbind(
    loading_table,
    temp
  )
}

print(loading_table)

                       Theme              Indicator PC1_Loading
ImpairedWaterBodies       LU    ImpairedWaterBodies   0.5886010
SmATS                     LU                  SmATS  -0.3402914
avg_impervious            LU         avg_impervious  -0.6422848
tree_canopy               LU            tree_canopy   0.3538656
asthma                    HV                 asthma  -0.6578798
Low_Birth_Weight          HV       Low_Birth_Weight  -0.5987070
Cardiovascular            HV         Cardiovascular  -0.4568852
ozone                     HE                  ozone  -0.2787390
PM2_5                     HE                  PM2_5   0.4281421
Diesel_PM                 HE              Diesel_PM   0.2792124
Pesticides                HE             Pesticides  -0.5318277
traffic                   HE                traffic   0.2782188
LST_Max                   HE                LST_Max  -0.2774607
LST_Mean                  HE               LST_Mean   0.4730841
Linguistic_Isol           SV        Linguistic_Isol  -0.4149005
Poverty                   SV                Poverty  -0.4417129
HousBurd                  SV               HousBurd  -0.3670935
children                  SV               children  -0.2218364
elderly                   SV                elderly   0.2676716
pop_density               SV            pop_density  -0.2001110
Total_outdoor             SV          Total_outdoor  -0.3039878
perc_no_vehicle_access    SV perc_no_vehicle_access  -0.2862734
energy_burden             SV          energy_burden  -0.4033915

# Correlation matrices for each theme
for (theme in names(themes)) {

  cat("\nCorrelation Matrix:", theme, "\n")

  correlation <- cor(
    index[themes[[theme]]],
    use = "pairwise.complete.obs"
  )

  print(round(correlation, 3))
}

Correlation Matrix: LU 
                    ImpairedWaterBodies  SmATS avg_impervious tree_canopy
ImpairedWaterBodies               1.000  0.048         -0.410       0.153
SmATS                             0.048  1.000          0.212      -0.188
avg_impervious                   -0.410  0.212          1.000      -0.042
tree_canopy                       0.153 -0.188         -0.042       1.000

Correlation Matrix: HV 
                 asthma Low_Birth_Weight Cardiovascular
asthma            1.000            0.428          0.275
Low_Birth_Weight  0.428            1.000          0.142
Cardiovascular    0.275            0.142          1.000

Correlation Matrix: HE 
            ozone  PM2_5 Diesel_PM Pesticides traffic LST_Max LST_Mean
ozone       1.000 -0.078     0.131      0.333  -0.140   0.143   -0.124
PM2_5      -0.078  1.000     0.380     -0.205   0.308  -0.215    0.209
Diesel_PM   0.131  0.380     1.000     -0.065   0.102  -0.018    0.252
Pesticides  0.333 -0.205    -0.065      1.000  -0.173   0.263   -0.571
traffic    -0.140  0.308     0.102     -0.173   1.000  -0.018    0.040
LST_Max     0.143 -0.215    -0.018      0.263  -0.018   1.000   -0.064
LST_Mean   -0.124  0.209     0.252     -0.571   0.040  -0.064    1.000

Correlation Matrix: SV 
                       Linguistic_Isol Poverty HousBurd children elderly pop_density Total_outdoor perc_no_vehicle_access energy_burden
Linguistic_Isol                  1.000   0.697    0.531    0.256  -0.400       0.182         0.545                  0.473         0.614
Poverty                          0.697   1.000    0.638    0.359  -0.299       0.268         0.467                  0.531         0.748
HousBurd                         0.531   0.638    1.000    0.210  -0.224       0.334         0.302                  0.446         0.519
children                         0.256   0.359    0.210    1.000  -0.398       0.186         0.109                  0.124         0.301
elderly                         -0.400  -0.299   -0.224   -0.398   1.000      -0.268        -0.359                 -0.043        -0.261
pop_density                      0.182   0.268    0.334    0.186  -0.268       1.000         0.042                  0.143         0.228
Total_outdoor                    0.545   0.467    0.302    0.109  -0.359       0.042         1.000                  0.202         0.486
perc_no_vehicle_access           0.473   0.531    0.446    0.124  -0.043       0.143         0.202                  1.000         0.359
energy_burden                    0.614   0.748    0.519    0.301  -0.261       0.228         0.486                  0.359         1.000

# Create output folder
outdir <- dirname(
  "C:/Users/noahs/Desktop/Ch1/SJV/Index.xlsx"
)

# Export Opportunity Index
write.csv(
  normalized,
  file.path(outdir, "Opportunity_Index_Results.csv"),
  row.names = FALSE
)

# Export explained variance
write.csv(
  variance_table,
  file.path(outdir, "PCA_Explained_Variance.csv"),
  row.names = FALSE
)

# Export PC1 loadings
write.csv(
  loading_table,
  file.path(outdir, "PCA_Loadings.csv"),
  row.names = FALSE
)

# Export missing data report
write.csv(
  data.frame(
    Indicator = names(missing_data),
    Missing_Count = as.integer(missing_data)
  ),
  file.path(outdir, "Missing_Data_Report.csv"),
  row.names = FALSE
)

hist(
  normalized$Opportunity_Index,
  breaks = 20,
  main = "Central Valley Opportunity Index",
  xlab = "Opportunity Index",
  ylab = "Number of Census Tracts"
)