Loading dependent variables

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
setwd("C:\\Users\\anami\\OneDrive\\Documents\\Health Disparities\\Project PA")
data1<- read.csv("C:\\Users\\anami\\OneDrive\\Documents\\Health Disparities\\Project PA\\data PLACES.csv")

Independent variable: Walkability Index

dataA <- read.csv("C:\\Users\\anami\\OneDrive\\Documents\\Health Disparities\\Project PA\\WalkabilityIndex_Tract_2019.csv" )
library(dplyr)
dataB<-dataA%>%
filter(StAbbr=="TX")
dataB

Calculating walkability index for Texas census tracts

library(ggplot2)
dataB$rank1 <- as.numeric(cut_number(dataB$D2A_EPHHM,20))
dataB$rank2 <- as.numeric(cut_number(dataB$D2B_E8MIXA,20))
dataB$rank3 <- as.numeric(cut_number(dataB$D3B,20))
library(dplyr)
dataC<-subset(dataB, select=c(TractID2019,Pop2018,rank1,rank2,rank3,D4A,NatWalkInd)) %>%
mutate(d4A =ifelse(D4A ==-99999.00,1, D4A))
dataC
library(ggplot2)
dataB$rank1 <- as.numeric(cut_number(dataB$D2A_EPHHM,20))
dataB$rank2 <- as.numeric(cut_number(dataB$D2B_E8MIXA,20))
dataB$rank3 <- as.numeric(cut_number(dataB$D3B,20))
library(dplyr)
dataC<-subset(dataB, select=c(TractID2019,Pop2018,rank1,rank2,rank3,D4A,NatWalkInd)) %>%
mutate(d4A =ifelse(D4A ==-99999.00,1, D4A))
dataC
dataC$rank4 <- as.numeric( cut(dataC$d4A,20, labels=F))
library(dplyr)
dataD<-dataC %>%
mutate(WI=(rank1/6)+(rank2/6)+(rank3/3)+(rank4/3))
names(dataD)
##  [1] "TractID2019" "Pop2018"     "rank1"       "rank2"       "rank3"      
##  [6] "D4A"         "NatWalkInd"  "d4A"         "rank4"       "WI"
library(dplyr)
dataE<-subset(dataD, select=c(TractID2019,WI)) %>%
group_by(TractID2019)
names(dataE)
## [1] "TractID2019" "WI"

Control Variables

file1<- read.csv("C:\\Users\\anami\\OneDrive\\Documents\\Health Disparities\\Project PA\\ACS2019.csv") 
file2<-subset(file1, select=c(geoid,totpop,fpop,mage,hispan,nh_white,nh_black,nh_asian))
names(file2)
## [1] "geoid"    "totpop"   "fpop"     "mage"     "hispan"   "nh_white" "nh_black"
## [8] "nh_asian"
nrow(file2)
## [1] 5265
file3<- read.csv("C:\\Users\\anami\\OneDrive\\Documents\\Health Disparities\\Project PA\\SVI2.csv") 
file4<-subset(file3, select=c(FIPS,RPL_THEMES))
names(file4)
## [1] "FIPS"       "RPL_THEMES"
nrow(file4)
## [1] 5254
summary(file4)
##       FIPS             RPL_THEMES       
##  Min.   :4.800e+10   Min.   :-999.0000  
##  1st Qu.:4.811e+10   1st Qu.:   0.2435  
##  Median :4.820e+10   Median :   0.4957  
##  Mean   :4.823e+10   Mean   :  -8.0606  
##  3rd Qu.:4.836e+10   3rd Qu.:   0.7478  
##  Max.   :4.851e+10   Max.   :   1.0000
nrow(file4)
## [1] 5254
library(dplyr)
file5<-file4 %>%
mutate(SVI =ifelse(RPL_THEMES ==-999.000,NA, RPL_THEMES))
file5
file6<-subset(file5, select=c(FIPS,SVI))
names(file6)
## [1] "FIPS" "SVI"
nrow(file6)
## [1] 5254
file7<-file6%>%
na.omit(file6)
nrow(file7)
summary(file7)
##       FIPS                SVI      
##  Min.   :4.800e+10   Min.   :0.00  
##  1st Qu.:4.811e+10   1st Qu.:0.25  
##  Median :4.820e+10   Median :0.50  
##  Mean   :4.823e+10   Mean   :0.50  
##  3rd Qu.:4.836e+10   3rd Qu.:0.75  
##  Max.   :4.851e+10   Max.   :1.00
nrow(file7)
## [1] 5209

Combining datasets:

combineA <-data2 %>%
  left_join(dataE, by = c("TractFIPS" = "TractID2019")) %>%
  dplyr::select(TractFIPS,WI,OBESITY_CrudePrev,LPA_CrudePrev)
names(combineA)
## [1] "TractFIPS"         "WI"                "OBESITY_CrudePrev"
## [4] "LPA_CrudePrev"
combineB <-combineA %>%
  left_join(file2, by = c("TractFIPS" = "geoid")) %>%
  dplyr::select(TractFIPS,WI,OBESITY_CrudePrev,LPA_CrudePrev,totpop,fpop,mage,hispan,nh_white,nh_black,nh_asian)
names(combineB)
##  [1] "TractFIPS"         "WI"                "OBESITY_CrudePrev"
##  [4] "LPA_CrudePrev"     "totpop"            "fpop"             
##  [7] "mage"              "hispan"            "nh_white"         
## [10] "nh_black"          "nh_asian"
combine1 <-combineB %>%
  left_join(file7, by = c("TractFIPS" = "FIPS")) %>%
  dplyr::select(TractFIPS,WI,OBESITY_CrudePrev, LPA_CrudePrev,totpop,fpop,mage,hispan,nh_white,nh_black,nh_asian,SVI)
names(combine1)
##  [1] "TractFIPS"         "WI"                "OBESITY_CrudePrev"
##  [4] "LPA_CrudePrev"     "totpop"            "fpop"             
##  [7] "mage"              "hispan"            "nh_white"         
## [10] "nh_black"          "nh_asian"          "SVI"
nrow(combine1)
## [1] 5222
summary(combine1)
##    TractFIPS               WI         OBESITY_CrudePrev LPA_CrudePrev  
##  Min.   :4.800e+10   Min.   : 1.000   Min.   :18.30     Min.   :10.10  
##  1st Qu.:4.811e+10   1st Qu.: 6.000   1st Qu.:32.73     1st Qu.:20.80  
##  Median :4.820e+10   Median : 8.333   Median :36.65     Median :26.50  
##  Mean   :4.823e+10   Mean   : 8.358   Mean   :36.83     Mean   :27.37  
##  3rd Qu.:4.836e+10   3rd Qu.:10.667   3rd Qu.:40.90     3rd Qu.:33.30  
##  Max.   :4.851e+10   Max.   :19.000   Max.   :55.30     Max.   :56.30  
##                                                                        
##      totpop           fpop           mage              hispan         
##  Min.   :   16   Min.   :    0   Length:5222        Length:5222       
##  1st Qu.: 3384   1st Qu.: 1682   Class :character   Class :character  
##  Median : 4812   Median : 2414   Mode  :character   Mode  :character  
##  Mean   : 5412   Mean   : 2724                                        
##  3rd Qu.: 6640   3rd Qu.: 3350                                        
##  Max.   :72041   Max.   :37626                                        
##                                                                       
##    nh_white           nh_black           nh_asian              SVI        
##  Length:5222        Length:5222        Length:5222        Min.   :0.0002  
##  Class :character   Class :character   Class :character   1st Qu.:0.2501  
##  Mode  :character   Mode  :character   Mode  :character   Median :0.5002  
##                                                           Mean   :0.5001  
##                                                           3rd Qu.:0.7501  
##                                                           Max.   :1.0000  
##                                                           NA's   :15

Variables description

library(ggplot2)
combine1$group <- as.numeric(cut_number(combine1$WI, 4))
combine1$group<-factor(combine1$group, levels = c(1, 2, 3,4), 
                      labels = c("Q1", "Q2", "Q3","Q4"))
table(combine1$group)
## 
##   Q1   Q2   Q3   Q4 
## 1332 1346 1246 1298
combine1a <- combine1 %>%
  rename(
    LPA = LPA_CrudePrev,
    Obesity = OBESITY_CrudePrev)
library(dplyr)
final1<-subset(combine1, select=c(WI,totpop,fpop,mage,hispan,nh_white,nh_black,nh_asian,group))%>%
  mutate(percentfemale=(fpop*100)/totpop)
names(final1)
##  [1] "WI"            "totpop"        "fpop"          "mage"         
##  [5] "hispan"        "nh_white"      "nh_black"      "nh_asian"     
##  [9] "group"         "percentfemale"
total_sum <- sum(final1$totpop)
print(total_sum)
## [1] 28260633
library(dplyr)
final1<-subset(combine1, select=c(WI,totpop,fpop,mage,hispan,nh_white,nh_black,nh_asian,SVI,group))%>%
  mutate(percentfemale=(fpop*100)/totpop)
names(final1)
##  [1] "WI"            "totpop"        "fpop"          "mage"         
##  [5] "hispan"        "nh_white"      "nh_black"      "nh_asian"     
##  [9] "SVI"           "group"         "percentfemale"
final1$mage <- as.numeric(final1$mage)
final1$hispan <- as.numeric(final1$hispan)
final1$nh_white <- as.numeric(final1$nh_white)
final1$nh_black <- as.numeric(final1$nh_black)
final1$nh_asian <- as.numeric(final1$nh_asian)
library(dplyr)
final2<-subset(final1, select=c(WI,percentfemale,mage,hispan,nh_white,nh_black,nh_asian,SVI,group))%>%
group_by(group)
names(final2)
## [1] "WI"            "percentfemale" "mage"          "hispan"       
## [5] "nh_white"      "nh_black"      "nh_asian"      "SVI"          
## [9] "group"
library(data.table)
## 
## Attaching package: 'data.table'
## The following objects are masked from 'package:dplyr':
## 
##     between, first, last
setnames(final2, old = c('percentfemale','mage','hispan','nh_white','nh_black','nh_asian'), new = c('%Female','Median Age','%Hispanics','%Non-hispanic White','%Non-hispanic Black','%Non-hispanic Asian'))
final2
## # A tibble: 5,222 × 9
## # Groups:   group [4]
##       WI `%Female` `Median Age` `%Hispanics` `%Non-hispanic White`
##    <dbl>     <dbl>        <dbl>        <dbl>                 <dbl>
##  1  6.83      51.7         21.3         20.3                  66.9
##  2  3.33      52.6         39.2         55.4                  43.7
##  3  7.67      49.1         37.4         29.7                  66.2
##  4  4         50.1         52           15.4                  80.6
##  5  8.17      47.2         40.5         42.6                  46  
##  6  9         52.4         31.7         58.9                  30.5
##  7 14.2       50.9         42           29.6                  64.4
##  8  7.33      48.7         34.3         94.9                   2.9
##  9  4.67      52.0         31.7         90.1                   7.1
## 10  7.83      51.3         31.5         23                    44  
## # ℹ 5,212 more rows
## # ℹ 4 more variables: `%Non-hispanic Black` <dbl>, `%Non-hispanic Asian` <dbl>,
## #   SVI <dbl>, group <fct>
final2 <- na.omit(final2)
library(arsenal) 
## Warning: package 'arsenal' was built under R version 4.4.2
table_1 <- tableby(group ~ ., data =final2) 
summary(table_1, title = "Descriptive Analysis I")
## 
## 
## Table: Descriptive Analysis I
## 
## |                            |   Q1 (N=1323)   |   Q2 (N=1343)   |   Q3 (N=1244)   |   Q4 (N=1297)   | Total (N=5207)  | p value|
## |:---------------------------|:---------------:|:---------------:|:---------------:|:---------------:|:---------------:|-------:|
## |**WI**                      |                 |                 |                 |                 |                 | < 0.001|
## |&nbsp;&nbsp;&nbsp;Mean (SD) |  4.276 (1.172)  |  7.222 (0.693)  |  9.504 (0.657)  | 12.629 (1.419)  |  8.365 (3.241)  |        |
## |&nbsp;&nbsp;&nbsp;Range     |  1.000 - 6.000  |  6.000 - 8.333  | 8.500 - 10.667  | 10.833 - 19.000 | 1.000 - 19.000  |        |
## |**%Female**                 |                 |                 |                 |                 |                 | < 0.001|
## |&nbsp;&nbsp;&nbsp;Mean (SD) | 49.655 (5.656)  | 50.485 (4.187)  | 50.807 (3.932)  | 50.390 (4.001)  | 50.327 (4.527)  |        |
## |&nbsp;&nbsp;&nbsp;Range     | 1.145 - 66.900  | 4.658 - 63.081  | 17.187 - 71.375 | 26.809 - 68.496 | 1.145 - 71.375  |        |
## |**Median Age**              |                 |                 |                 |                 |                 | < 0.001|
## |&nbsp;&nbsp;&nbsp;Mean (SD) | 38.300 (7.832)  | 35.766 (6.341)  | 34.846 (5.931)  | 35.514 (6.175)  | 36.127 (6.749)  |        |
## |&nbsp;&nbsp;&nbsp;Range     | 18.800 - 73.700 | 19.800 - 60.300 | 19.600 - 64.700 | 19.800 - 61.700 | 18.800 - 73.700 |        |
## |**%Hispanics**              |                 |                 |                 |                 |                 | < 0.001|
## |&nbsp;&nbsp;&nbsp;Mean (SD) | 33.477 (28.184) | 38.632 (26.977) | 44.939 (28.221) | 41.626 (27.782) | 39.575 (28.093) |        |
## |&nbsp;&nbsp;&nbsp;Range     | 0.000 - 100.000 | 0.600 - 100.000 | 0.000 - 100.000 | 0.000 - 99.900  | 0.000 - 100.000 |        |
## |**%Non-hispanic White**     |                 |                 |                 |                 |                 | < 0.001|
## |&nbsp;&nbsp;&nbsp;Mean (SD) | 54.310 (28.583) | 42.348 (26.484) | 34.523 (25.301) | 38.847 (25.518) | 42.646 (27.517) |        |
## |&nbsp;&nbsp;&nbsp;Range     | 0.000 - 100.000 | 0.000 - 94.400  | 0.000 - 89.800  | 0.000 - 98.400  | 0.000 - 100.000 |        |
## |**%Non-hispanic Black**     |                 |                 |                 |                 |                 | < 0.001|
## |&nbsp;&nbsp;&nbsp;Mean (SD) | 8.646 (13.445)  | 12.925 (16.526) | 14.036 (17.738) | 11.210 (13.402) | 11.676 (15.500) |        |
## |&nbsp;&nbsp;&nbsp;Range     | 0.000 - 95.600  | 0.000 - 96.800  | 0.000 - 92.100  | 0.000 - 81.300  | 0.000 - 96.800  |        |
## |**%Non-hispanic Asian**     |                 |                 |                 |                 |                 | < 0.001|
## |&nbsp;&nbsp;&nbsp;Mean (SD) |  1.502 (3.414)  |  3.955 (6.992)  |  4.520 (7.157)  |  6.031 (8.701)  |  3.984 (7.023)  |        |
## |&nbsp;&nbsp;&nbsp;Range     | 0.000 - 42.800  | 0.000 - 61.200  | 0.000 - 54.700  | 0.000 - 75.800  | 0.000 - 75.800  |        |
## |**SVI**                     |                 |                 |                 |                 |                 | < 0.001|
## |&nbsp;&nbsp;&nbsp;Mean (SD) |  0.470 (0.254)  |  0.504 (0.288)  |  0.554 (0.296)  |  0.475 (0.308)  |  0.500 (0.289)  |        |
## |&nbsp;&nbsp;&nbsp;Range     |  0.000 - 0.996  |  0.003 - 0.999  |  0.002 - 1.000  |  0.000 - 1.000  |  0.000 - 1.000  |        |
final3<-subset(combine1, select=c(OBESITY_CrudePrev,LPA_CrudePrev,group))
names(final3)
## [1] "OBESITY_CrudePrev" "LPA_CrudePrev"     "group"
final3 <- na.omit(final3)
options(digits=2)
library(arsenal) 
table_2 <- tableby(group ~ ., data =final3) 
summary(table_2)
## 
## 
## |                            |   Q1 (N=1332)   |   Q2 (N=1346)   |   Q3 (N=1246)   |   Q4 (N=1298)   | Total (N=5222)  | p value|
## |:---------------------------|:---------------:|:---------------:|:---------------:|:---------------:|:---------------:|-------:|
## |**OBESITY_CrudePrev**       |                 |                 |                 |                 |                 | < 0.001|
## |&nbsp;&nbsp;&nbsp;Mean (SD) | 37.376 (4.641)  | 36.907 (5.820)  | 37.690 (6.372)  | 35.365 (6.778)  | 36.830 (6.009)  |        |
## |&nbsp;&nbsp;&nbsp;Range     | 19.100 - 52.700 | 22.100 - 52.900 | 21.600 - 55.300 | 18.300 - 53.100 | 18.300 - 55.300 |        |
## |**LPA_CrudePrev**           |                 |                 |                 |                 |                 | < 0.001|
## |&nbsp;&nbsp;&nbsp;Mean (SD) | 27.345 (6.641)  | 27.133 (8.374)  | 28.904 (8.882)  | 26.167 (9.695)  | 27.369 (8.509)  |        |
## |&nbsp;&nbsp;&nbsp;Range     | 11.500 - 51.900 | 11.500 - 52.200 | 10.900 - 51.600 | 10.100 - 56.300 | 10.100 - 56.300 |        |

Histogram

library(ggplot2)

ggplot(combine1, aes(x = WI, fill = group)) + 
  geom_histogram(color = "black", bins = 30) + 
  scale_fill_manual(
    values = c("#d9f0d3", "#a6dba0", "#5aae61", "#1b7837"),
    name = "Group",
    labels = c("Q1", "Q2", "Q3", "Q4")
  ) +
  theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.title = element_text(size = 12),
    legend.text = element_text(size = 10),
    legend.spacing.x = unit(1, 'cm'),    # <-- Increase space between legend keys
    axis.title.x = element_text(size = 14),
    axis.title.y = element_text(size = 14),
    axis.text = element_text(size = 10),
    plot.title = element_text(size = 16, face = "bold")
  ) +
  labs(
    x = "Walkability Index",
    y = "Number of Census Tracts",
    title = "Distribution of Walkability Index by Quartiles"
  )

Scatterplots

library(ggplot2)

P1 <- ggplot(combine1a, aes(WI, Obesity)) + 
  geom_point() +
  geom_smooth(method = 'lm') +
  labs(x = "Walkability Index", y = "Obesity") +
  theme(
    axis.title.x = element_text(size = 12.5),  # Axis title size
    axis.title.y = element_text(size = 12.5),
    axis.text.x = element_text(size = 11),     # Axis number (tick) size
    axis.text.y = element_text(size = 11)
  )

P1
## `geom_smooth()` using formula = 'y ~ x'

library(ggplot2)

P2 <- ggplot(combine1a, aes(WI, LPA)) + 
  geom_point() +
  geom_smooth(method = 'lm') +
  labs(x = "Walkability Index", y = "Lack of Physical Activity") +
  theme(
    axis.title.x = element_text(size = 12.5),  # Axis title size
    axis.title.y = element_text(size = 12.5),
    axis.text.x = element_text(size = 11),     # Axis number (tick) size
    axis.text.y = element_text(size = 11)
  )

P2
## `geom_smooth()` using formula = 'y ~ x'

library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
plot_list <- list(P1,P2) 
do.call("grid.arrange", c(plot_list, ncol = 1))   
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

Spearman’s Correlation Test

cor.test(combine1$OBESITY_CrudePrev, combine1$WI, method = 'spearman',exact=FALSE)
## 
##  Spearman's rank correlation rho
## 
## data:  combine1$OBESITY_CrudePrev and combine1$WI
## S = 3e+10, p-value <2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##   rho 
## -0.11
cor.test(combine1$LPA_CrudePrev, combine1$WI, method = 'spearman',exact=FALSE)
## 
##  Spearman's rank correlation rho
## 
## data:  combine1$LPA_CrudePrev and combine1$WI
## S = 3e+10, p-value = 6e-06
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##    rho 
## -0.063

Regression models

library(broom)
fit<-lm(OBESITY_CrudePrev~WI, data=combine1)
summary(fit)
## 
## Call:
## lm(formula = OBESITY_CrudePrev ~ WI, data = combine1)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -18.86  -4.18  -0.44   4.04  18.63 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  38.9129     0.2278   170.8   <2e-16 ***
## WI           -0.2492     0.0254    -9.8   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6 on 5220 degrees of freedom
## Multiple R-squared:  0.0181, Adjusted R-squared:  0.0179 
## F-statistic: 96.1 on 1 and 5220 DF,  p-value: <2e-16
library(broom)
fit<-lm(LPA_CrudePrev~WI, data=combine1)
summary(fit)
## 
## Call:
## lm(formula = LPA_CrudePrev ~ WI, data = combine1)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -16.46  -6.54  -1.01   6.00  29.57 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  28.7148     0.3250   88.36  < 2e-16 ***
## WI           -0.1610     0.0362   -4.44  9.2e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.5 on 5220 degrees of freedom
## Multiple R-squared:  0.00376,    Adjusted R-squared:  0.00357 
## F-statistic: 19.7 on 1 and 5220 DF,  p-value: 9.16e-06

Multivariate Analysis

library(dplyr)
combine2<-combine1%>%
  mutate(percentfemale=(fpop*100)/totpop)
names(combine2)
##  [1] "TractFIPS"         "WI"                "OBESITY_CrudePrev"
##  [4] "LPA_CrudePrev"     "totpop"            "fpop"             
##  [7] "mage"              "hispan"            "nh_white"         
## [10] "nh_black"          "nh_asian"          "SVI"              
## [13] "group"             "percentfemale"
combine2$mage <- as.numeric(combine2$mage)
combine2$hispan <- as.numeric(combine2$hispan)
combine2$nh_white <- as.numeric(combine2$nh_white)
combine2$nh_black <- as.numeric(combine2$nh_black)
combine2$nh_asian <- as.numeric(combine2$nh_asian)
library(broom)
fitA<-lm(OBESITY_CrudePrev~WI+mage+percentfemale+SVI+hispan+nh_white+nh_black, data=combine2)
summary(fitA)
## 
## Call:
## lm(formula = OBESITY_CrudePrev ~ WI + mage + percentfemale + 
##     SVI + hispan + nh_white + nh_black, data = combine2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -11.094  -1.719   0.072   1.804   9.583 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   13.18849    0.68368   19.29   <2e-16 ***
## WI            -0.23721    0.01235  -19.21   <2e-16 ***
## mage           0.07765    0.00686   11.32   <2e-16 ***
## percentfemale -0.01449    0.00834   -1.74    0.082 .  
## SVI            9.83863    0.20824   47.25   <2e-16 ***
## hispan         0.22225    0.00568   39.10   <2e-16 ***
## nh_white       0.15197    0.00575   26.43   <2e-16 ***
## nh_black       0.28808    0.00638   45.17   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.7 on 5199 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.797,  Adjusted R-squared:  0.797 
## F-statistic: 2.92e+03 on 7 and 5199 DF,  p-value: <2e-16
library(broom)
fitB<-lm(LPA_CrudePrev~WI+mage+percentfemale+SVI+hispan+nh_white+nh_black, data=combine2)
summary(fitB)
## 
## Call:
## lm(formula = LPA_CrudePrev ~ WI + mage + percentfemale + SVI + 
##     hispan + nh_white + nh_black, data = combine2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -13.771  -2.239  -0.037   2.181  15.633 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    4.22342    0.89626    4.71  2.5e-06 ***
## WI            -0.16318    0.01618  -10.08  < 2e-16 ***
## mage           0.13163    0.00899   14.63  < 2e-16 ***
## percentfemale -0.01347    0.01093   -1.23     0.22    
## SVI           20.83245    0.27298   76.31  < 2e-16 ***
## hispan         0.14229    0.00745   19.10  < 2e-16 ***
## nh_white       0.06440    0.00754    8.54  < 2e-16 ***
## nh_black       0.14044    0.00836   16.80  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.5 on 5199 degrees of freedom
##   (15 observations deleted due to missingness)
## Multiple R-squared:  0.827,  Adjusted R-squared:  0.827 
## F-statistic: 3.55e+03 on 7 and 5199 DF,  p-value: <2e-16

Mediation Analysis:

library(broom)
model1<-lm(LPA_CrudePrev~WI, data=combine1)
summary(model1)
## 
## Call:
## lm(formula = LPA_CrudePrev ~ WI, data = combine1)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -16.46  -6.54  -1.01   6.00  29.57 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  28.7148     0.3250   88.36  < 2e-16 ***
## WI           -0.1610     0.0362   -4.44  9.2e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.5 on 5220 degrees of freedom
## Multiple R-squared:  0.00376,    Adjusted R-squared:  0.00357 
## F-statistic: 19.7 on 1 and 5220 DF,  p-value: 9.16e-06
library(broom)
model2<-lm(OBESITY_CrudePrev~WI, data=combine1)
summary(model2)
## 
## Call:
## lm(formula = OBESITY_CrudePrev ~ WI, data = combine1)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -18.86  -4.18  -0.44   4.04  18.63 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  38.9129     0.2278   170.8   <2e-16 ***
## WI           -0.2492     0.0254    -9.8   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6 on 5220 degrees of freedom
## Multiple R-squared:  0.0181, Adjusted R-squared:  0.0179 
## F-statistic: 96.1 on 1 and 5220 DF,  p-value: <2e-16
library(broom)
model3<-lm(OBESITY_CrudePrev~LPA_CrudePrev, data=combine1)
summary(model3)
## 
## Call:
## lm(formula = OBESITY_CrudePrev ~ LPA_CrudePrev, data = combine1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -11.971  -1.759  -0.195   1.571   9.135 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   19.45791    0.12274     159   <2e-16 ***
## LPA_CrudePrev  0.63474    0.00428     148   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.6 on 5220 degrees of freedom
## Multiple R-squared:  0.808,  Adjusted R-squared:  0.808 
## F-statistic: 2.2e+04 on 1 and 5220 DF,  p-value: <2e-16
library(psych)
## Warning: package 'psych' was built under R version 4.4.2
## 
## Attaching package: 'psych'
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
modelA <- mediate(Obesity ~WI+ (LPA), data =combine1a)

plot(modelA)

MAPPING

library(dplyr)
library(tigris)
## To enable caching of data, set `options(tigris_use_cache = TRUE)`
## in your R script or .Rprofile.
library(leaflet)
tx_censustracts <- tracts(state = "TX")
## Retrieving data for the year 2022
tx_censustracts_cb <- tracts("TX", cb = TRUE)
## Retrieving data for the year 2022
plot(tx_censustracts$geometry)

plot(tx_censustracts_cb$geometry)

dmap <- tx_censustracts_cb %>%  
    mutate(GEOID = as.numeric(GEOID))
library(dplyr)  # Make sure dplyr is loaded

dmap1 <- dmap %>%
  left_join(combine1a, by = c("GEOID" = "TractFIPS")) %>%
  dplyr::select(GEOID, WI, Obesity, LPA)

names(dmap1)
## [1] "GEOID"    "WI"       "Obesity"  "LPA"      "geometry"
library(tmap)
## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')
m1 <- tm_shape(dmap1) + 
  tm_fill(
    col = "Obesity", 
    style = "quantile", 
    n = 4, 
    palette = "Reds", 
    colorNA = "lightgray",    
    textNA = "Missing",
    title = "% Obesity"
  ) +
  tm_borders(col = "gray50", lwd = 0.5) +
  
  tm_scale_bar(
    position = c(0.78, 0.04),  
    breaks = c(0, 100, 200)
  ) +
  
  tm_compass(
    type = "4star", 
    size = 2.5,              
    position = c(0.8, 0.12) 
  ) +
  
  tm_layout(
    title = "Obesity",
    title.position = c(0.65, 0.95),
    title.size = 1.2,
    title.fontface = "bold",
    frame = TRUE,
    legend.position = c("left", "bottom"),
    legend.title.size = 0.9,
    legend.text.size = 0.8
  )

m1

library(tmap)

m2 <- tm_shape(dmap1) + 
  tm_fill(
    col = "LPA", 
    style = "quantile", 
    n = 4, 
    palette = "Blues", 
    colorNA = "lightgray",    
    textNA = "Missing",
    title = "% Lack of Physical Activity"
  ) +
  tm_borders(col = "gray50", lwd = 0.5) +
  
  tm_scale_bar(
    position = c(0.78, 0.04),  # Moved scale bar more to the right
    breaks = c(0, 100, 200)
  ) +
  
  tm_compass(
    type = "4star", 
    size = 2.5,              # Bigger compass
    position = c(0.8, 0.12) # Slightly higher than scale bar, more spacing
  ) +
  
  tm_layout(
    title = "Lack of Physical Activity",
    title.position = c(0.52, 0.95),
    title.size = 1.1,
    title.fontface = "bold",
    frame = TRUE,
    legend.position = c("left", "bottom"),
    legend.title.size = 0.9,
    legend.text.size = 0.8
  )

m2

library(tmap)

m3 <- tm_shape(dmap1) + 
  tm_fill(
    col = "WI", 
    style = "quantile", 
    n = 4, 
    palette = "Greens", 
    colorNA = "lightgray",    
    textNA = "Missing",
    title = "Walkability Index"
  ) +
  tm_borders(col = "gray50", lwd = 0.5) +
  
  tm_scale_bar(
    position = c(0.78, 0.04),  # Moved scale bar more to the right
    breaks = c(0, 100, 200)
  ) +
  
  tm_compass(
    type = "4star", 
    size = 2.5,              # Bigger compass
    position = c(0.8, 0.12) # Slightly higher than scale bar, more spacing
  ) +
  
  tm_layout(
    title = "Walkability",
    title.position = c(0.65, 0.95),
    title.size = 1.2,
    title.fontface = "bold",
    frame = TRUE,
    legend.position = c("left", "bottom"),
    legend.title.size = 0.9,
    legend.text.size = 0.8
  )

m3