R Markdown

Stats methods: Performed Chi-square tests of independence to determine the potential for an association between two variables. The null hypothesis states that there is no relationship between the two variables. Alpha 0.05.

We then calculated Pearson residuals to determine which categories contribute the most to the total chi-squared calculation. The highest absolute value = most contributing cells to the total Chi-square score.Positive values in cells specify an attraction (positive association) between the corresponding row and column variables.Negative residuals imply a repulsion (negative association) between the corresponding row and column variables.

load(file = "/Users/carissawonkka/Documents/RangeWUI-main/DataSummaries.Rdata")

#function to get contingency tables
make_matrix <- function(df,rownames = NULL){
  my_matrix <-  as.matrix(df)
  if(!is.null(rownames))
    rownames(my_matrix) = rownames
  my_matrix
}

Test to see if total WUI differs by cover classes

propWUI <-
  DataSummaries$OverallRangeWUI %>%
  select(!distance) %>%
  group_by(class) %>%
  summarise(across(everything(), list(sum))) %>%
  mutate(totalWUI = rowSums(.[3:4]))%>%
  mutate(propWUI = totalWUI / sum(totalWUI)) %>%
  select(!NotWUI_1:totalWUI)
# Simulate response data based on multinomial distribution
# Set simulation parameters   
obs.size = 100 # effective N
resp.cats = propWUI$class # category names
# Get observed probabilities but drop categories with < 5% total 
obs.prob <- propWUI %>% 
  filter(propWUI >= 0.05) %>%
  select(propWUI) %>% 
  as_vector() 
exp.prob = 1/length(obs.prob) # null hypothesis (even % across categories)
n.sims = 1000  # total interactions for simulation
# Create empty array into which simulations are added
cat.sims<-array(NA,c(n.sims,length(obs.prob)))
cat.sims <- rmultinom(n.sims, size=obs.size, prob=obs.prob) %>%
  t() %>%
  as.data.frame()
names(cat.sims)<-(resp.cats)
# Create mpty tibble to stash stats results
PearsonsOut <- tibble() 
# Run Chi-squared test on simulated data
for (i in 1:nrow(cat.sims)){  
  # Grab a row of simulated data  
  sim.d <- cat.sims %>% slice(i)
  # run Chi-squared test
  x2 <- chisq.test(sim.d)
  # calculate effect size
  ef <- effectsize::effectsize(x2, type = 'pearsons_c', alternative = 'two')
  # extract residuals, combine with effect size, export results of i-th run
  x2$residuals %>% 
    as_tibble()  %>%
    mutate(levels = resp.cats) %>%
    pivot_wider(names_from = levels, 
                values_from = value) %>%
    add_column(PearsonsC = ef$Pearsons_c, .before = 1) %>%
    bind_rows(PearsonsOut) -> PearsonsOut
} 
# Calculate median and 95% CI from N simulation results 
PearsonsOut %>%
  reframe(across(c(PearsonsC:Rangeland), ~quantile(.x, prob=c(0.025,0.5,0.975))) ) %>%
  add_column(quantile = c('ciL', 'median', 'ciU'), .before = 1) %>%
  pivot_longer(values_to = 'value', 
               names_to = 'factor', 
               -quantile) %>%
  pivot_wider(names_from = quantile, 
              values_from = value) %>%
  pander("Pearson's C = effect size for the overall Chi-squared test. Remaining rows are Pearson's contrasts for each factor level.")
Pearson’s C = effect size for the overall Chi-squared test. Remaining rows are Pearson’s contrasts for each factor level.
factor ciL median ciU
PearsonsC 0.3478 0.4891 0.6139
Developed -4.4 -3.4 -2.2
Farmland -2.2 -0.6 1
Other -1.8 -0.2 1.4
Rangeland 2.4 4.2 6.4

Test to see if the proportion of WUI types differs across cover classes - simulation method

propWUI <-
  DataSummaries$OverallRangeWUI %>%
  select(!distance) %>%
  group_by(class) %>%
  summarise(across(everything(), list(sum))) %>%
  mutate(totalWUI = rowSums(.[3:4]))%>%
  mutate(propWUIinterface = Interface_1 / sum(Interface_1 )) %>%
  mutate(propWUIintermix = Intermix_1 / sum(Intermix_1 )) %>%
  select(class,propWUIinterface, propWUIintermix)
# Simulate response data based on multinomial distribution
# Set simulation parameters   
obs.size = 100 # effective N
propWUI = propWUI %>%
  filter(!class %in% 'Developed')
resp.cats = propWUI$class # category names
# Get observed probabilities but drop categories with < 5% total 
  obs.prob.face <- 
    propWUI %>% 
      select(propWUIinterface) %>% 
        as_vector() 
  obs.prob.mix <- 
    propWUI %>% 
      select(propWUIintermix) %>% 
      as_vector() 
  exp.prob = 1/length(obs.prob.face) # null hypothesis (even % across categories)
  n.sims = 1000  # total interactions for simulation
# Create empty array into which simulations are added for interface
  cat.sims.face <- tibble() 
  for (i in 1:n.sims) { 
  rmultinom(1, size=jitter(obs.size, 10), prob=obs.prob.face) %>%
                      t() %>%
                      as_tibble() %>%
      bind_rows(cat.sims.face) -> cat.sims.face
  }
  names(cat.sims.face)<-(resp.cats)
  cat.sims.face %<>% add_column( type = 'interface', .before = 1)
# Create empty array into which simulations are added for intermix
  cat.sims.mix <- tibble() 
  for (i in 1:n.sims) { 
    rmultinom(1, size=jitter(obs.size, 10), prob=obs.prob.mix) %>%
      t() %>%
      as_tibble() %>%
      bind_rows(cat.sims.mix) -> cat.sims.mix
  }
  names(cat.sims.mix)<-(resp.cats)
  cat.sims.mix %<>% add_column( type = 'intermix', .before = 1)
# Create empty tibble to stash stats results
  PearsonsOut <- tibble() 
# Run Chi-squared test on simulated data
  for (i in 1:n.sims) {  
    # Grab a row of simulated data for interface 
      sim.d.face <- cat.sims.face %>% slice(i)
    # Grab a row of simulated data for intermix
      sim.d.mix <- cat.sims.mix %>% slice(i)
    # Combine intermix and face into one table
      sim.d <- full_join(by = c('type', "Farmland", "Other", "Rangeland"), 
                         sim.d.face, 
                         sim.d.mix) %>% 
        pivot_longer(!type, names_to = "class", values_to = "count") %>%  
        pivot_wider(names_from=type, values_from=count) 
    # run Chi-squared test
      x2 <- chisq.test(sim.d %>% select(-class))
    # calculate effect size
      ef <- effectsize::effectsize(x2, type = 'pearsons_c', alternative = 'two')
    # extract residuals, combine with effect size, export results of i-th run
     x2= x2$residuals %>% 
        as_tibble()  %>%
        mutate(type = sim.d$class) %>%
        pivot_longer(names_to = 'class', 
                     values_to = 'residual', 
                    -type) %>%
        pivot_wider(names_from = c(type, class), 
                    values_from = residual) %>%
        add_column(PearsonsC = ef$Pearsons_c, .before = 1) %>%
        bind_rows(PearsonsOut) -> PearsonsOut
  } 
# Calculate median and 95% CI from N simulation results 
PearsonsOut %>% 
  reframe(across(c(PearsonsC:Rangeland_intermix), ~quantile(.x, prob=c(0.025,0.5,0.975))) ) %>%
  add_column(quantile = c('ciL', 'median', 'ciU'), .before = 1) %>%
  pivot_longer(values_to = 'value', 
               names_to = 'factor', 
               -quantile) %>%
  pivot_wider(names_from = quantile, 
              values_from = value) %>%
  separate(factor, into=c('type', 'class'), sep = '_') %>%
  pander("Pearson's C = effect size for the overall Chi-squared test. Remaining rows are Pearson's contrasts for each factor level.")
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 1 rows [1].
Pearson’s C = effect size for the overall Chi-squared test. Remaining rows are Pearson’s contrasts for each factor level.
type class ciL median ciU
PearsonsC NA 0.4476 0.5205 0.5785
Farmland interface 3.635 4.709 5.778
Farmland intermix -5.565 -4.729 -3.832
Other interface -3.543 -2.443 -1.233
Other intermix 1.215 2.456 3.57
Rangeland interface -3.817 -2.875 -1.859
Rangeland intermix 1.87 2.86 3.952

Test to see if the proportion of WUI types differs across cover classes broken down by ecoregion. Do all ecoregions follow the broader pattern? simulation

  #Great Plains

propWUI <-
  DataSummaries$EcoregionRangeWUI %>%
  select(!distance) %>%
  filter(L1=="Great Plains")%>%
  select(!L1) %>%
  group_by(class) %>%
  summarise(across(everything(), list(sum))) %>%
  mutate(totalWUI = rowSums(.[3:4]))%>%
  mutate(propWUIinterface = Interface_1 / sum(Interface_1 )) %>%
  mutate(propWUIintermix = Intermix_1 / sum(Intermix_1 )) %>%
  select(class,propWUIinterface, propWUIintermix)
# Simulate response data based on multinomial distribution
# Set simulation parameters   
obs.size = 100 # effective N
propWUI = propWUI %>%
  filter(!class %in% 'Other')
resp.cats = propWUI$class # category names
# Get observed probabilities but drop categories with < 5% total 
  obs.prob.face <- 
    propWUI %>% 
      select(propWUIinterface) %>% 
        as_vector() 
  obs.prob.mix <- 
    propWUI %>% 
      select(propWUIintermix) %>% 
      as_vector() 
  exp.prob = 1/length(obs.prob.face) # null hypothesis (even % across categories)
  n.sims = 1000  # total interactions for simulation
# Create empty array into which simulations are added for interface
  cat.sims.face <- tibble() 
  for (i in 1:n.sims) { 
  rmultinom(1, size=jitter(obs.size, 10), prob=obs.prob.face) %>%
                      t() %>%
                      as_tibble() %>%
      bind_rows(cat.sims.face) -> cat.sims.face
  }
  names(cat.sims.face)<-(resp.cats)
  cat.sims.face %<>% add_column( type = 'face', .before = 1)
# Create empty array into which simulations are added for intermix
  cat.sims.mix <- tibble() 
  for (i in 1:n.sims) { 
    rmultinom(1, size=jitter(obs.size, 10), prob=obs.prob.mix) %>%
      t() %>%
      as_tibble() %>%
      bind_rows(cat.sims.mix) -> cat.sims.mix
  }
  names(cat.sims.mix)<-(resp.cats)
  cat.sims.mix %<>% add_column( type = 'mix', .before = 1)
# Create empty tibble to stash stats results
  PearsonsOut <- tibble() 
# Run Chi-squared test on simulated data
  for (i in 1:n.sims) {  
    # Grab a row of simulated data for interface 
      sim.d.face <- cat.sims.face %>% slice(i)
    # Grab a row of simulated data for intermix
      sim.d.mix <- cat.sims.mix %>% slice(i)
    # Combine intermix and face into one table
      sim.d <- full_join(by = c('type', 'Developed',"Farmland", "Rangeland"), 
                         sim.d.face, 
                         sim.d.mix) %>% 
        pivot_longer(!type, names_to = "class", values_to = "count") %>%  
        pivot_wider(names_from=type, values_from=count) 
    # run Chi-squared test
      x2 <- chisq.test(sim.d %>% select(-class))
    # calculate effect size
      ef <- effectsize::effectsize(x2, type = 'pearsons_c', alternative = 'two')
    # extract residuals, combine with effect size, export results of i-th run
     x2= x2$residuals %>% 
        as_tibble()  %>%
        mutate(type = sim.d$class) %>%
        pivot_longer(names_to = 'class', 
                     values_to = 'residual', 
                    -type) %>%
        pivot_wider(names_from = c(type, class), 
                    values_from = residual) %>%
        add_column(PearsonsC = ef$Pearsons_c, .before = 1) %>%
        bind_rows(PearsonsOut) -> PearsonsOut
  } 
## Warning in chisq.test(sim.d %>% select(-class)): Chi-squared approximation may
## be incorrect

## Warning in chisq.test(sim.d %>% select(-class)): Chi-squared approximation may
## be incorrect
# Calculate median and 95% CI from N simulation results 
PearsonsOut %>% 
  reframe(across(c(PearsonsC:Rangeland_mix), ~quantile(.x, prob=c(0.025,0.5,0.975))) ) %>%
  add_column(quantile = c('ciL', 'median', 'ciU'), .before = 1) %>%
  pivot_longer(values_to = 'value', 
               names_to = 'factor', 
               -quantile) %>%
  pivot_wider(names_from = quantile, 
              values_from = value) %>%
  separate(factor, into=c('type', 'class'), sep = '_') %>%
  pander("Pearson's C = effect size for the overall Chi-squared test. Remaining rows are Pearson's contrasts for each factor level.")
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 1 rows [1].
Pearson’s C = effect size for the overall Chi-squared test. Remaining rows are Pearson’s contrasts for each factor level.
type class ciL median ciU
PearsonsC NA 0.436 0.5198 0.5852
Developed face 0.5234 1.629 2.812
Developed mix -2.709 -1.644 -0.5418
Farmland face 2.985 4.009 5.12
Farmland mix -4.892 -3.983 -2.974
Rangeland face -5.036 -4.183 -3.241
Rangeland mix 3.058 4.158 5.124
 #North American Deserts

propWUI <-
  DataSummaries$EcoregionRangeWUI %>%
  select(!distance) %>%
  filter(L1=="North American Deserts")%>%
  select(!L1) %>%
  group_by(class) %>%
  summarise(across(everything(), list(sum))) %>%
  mutate(totalWUI = rowSums(.[3:4]))%>%
  mutate(propWUIinterface = Interface_1 / sum(Interface_1 )) %>%
  mutate(propWUIintermix = Intermix_1 / sum(Intermix_1 )) %>%
  select(class,propWUIinterface, propWUIintermix)
# Simulate response data based on multinomial distribution
# Set simulation parameters   
obs.size = 100 # effective N
propWUI = propWUI %>%
  filter(!class %in% 'Developed')
resp.cats = propWUI$class # category names
# Get observed probabilities but drop categories with < 5% total 
  obs.prob.face <- 
    propWUI %>% 
      select(propWUIinterface) %>% 
        as_vector() 
  obs.prob.mix <- 
    propWUI %>% 
      select(propWUIintermix) %>% 
      as_vector() 
  exp.prob = 1/length(obs.prob.face) # null hypothesis (even % across categories)
  n.sims = 1000  # total interactions for simulation
# Create empty array into which simulations are added for interface
  cat.sims.face <- tibble() 
  for (i in 1:n.sims) { 
  rmultinom(1, size=jitter(obs.size, 10), prob=obs.prob.face) %>%
                      t() %>%
                      as_tibble() %>%
      bind_rows(cat.sims.face) -> cat.sims.face
  }
  names(cat.sims.face)<-(resp.cats)
  cat.sims.face %<>% add_column( type = 'face', .before = 1)
# Create empty array into which simulations are added for intermix
  cat.sims.mix <- tibble() 
  for (i in 1:n.sims) { 
    rmultinom(1, size=jitter(obs.size, 10), prob=obs.prob.mix) %>%
      t() %>%
      as_tibble() %>%
      bind_rows(cat.sims.mix) -> cat.sims.mix
  }
  names(cat.sims.mix)<-(resp.cats)
  cat.sims.mix %<>% add_column( type = 'mix', .before = 1)
# Create empty tibble to stash stats results
  PearsonsOut <- tibble() 
# Run Chi-squared test on simulated data
  for (i in 1:n.sims) {  
    # Grab a row of simulated data for interface 
      sim.d.face <- cat.sims.face %>% slice(i)
    # Grab a row of simulated data for intermix
      sim.d.mix <- cat.sims.mix %>% slice(i)
    # Combine intermix and face into one table
      sim.d <- full_join(by = c('type',"Farmland",'Other', "Rangeland"), 
                         sim.d.face, 
                         sim.d.mix) %>% 
        pivot_longer(!type, names_to = "class", values_to = "count") %>%  
        pivot_wider(names_from=type, values_from=count) 
    # run Chi-squared test
      x2 <- chisq.test(sim.d %>% select(-class))
    # calculate effect size
      ef <- effectsize::effectsize(x2, type = 'pearsons_c', alternative = 'two')
    # extract residuals, combine with effect size, export results of i-th run
     x2= x2$residuals %>% 
        as_tibble()  %>%
        mutate(type = sim.d$class) %>%
        pivot_longer(names_to = 'class', 
                     values_to = 'residual', 
                    -type) %>%
        pivot_wider(names_from = c(type, class), 
                    values_from = residual) %>%
        add_column(PearsonsC = ef$Pearsons_c, .before = 1) %>%
        bind_rows(PearsonsOut) -> PearsonsOut
  } 
# Calculate median and 95% CI from N simulation results 
PearsonsOut %>% 
  reframe(across(c(PearsonsC:Rangeland_mix), ~quantile(.x, prob=c(0.025,0.5,0.975))) ) %>%
  add_column(quantile = c('ciL', 'median', 'ciU'), .before = 1) %>%
  pivot_longer(values_to = 'value', 
               names_to = 'factor', 
               -quantile) %>%
  pivot_wider(names_from = quantile, 
              values_from = value) %>%
  separate(factor, into=c('type', 'class'), sep = '_') %>%
  pander("Pearson's C = effect size for the overall Chi-squared test. Remaining rows are Pearson's contrasts for each factor level.")
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 1 rows [1].
Pearson’s C = effect size for the overall Chi-squared test. Remaining rows are Pearson’s contrasts for each factor level.
type class ciL median ciU
PearsonsC NA 0.4382 0.5141 0.5754
Farmland face 3.646 4.696 5.809
Farmland mix -5.524 -4.72 -3.789
Other face -2.565 -1.393 -0.158
Other mix 0.1599 1.392 2.609
Rangeland face -4.286 -3.321 -2.362
Rangeland mix 2.314 3.349 4.33
#North American Deserts

propWUI <-
  DataSummaries$EcoregionRangeWUI %>%
  select(!distance) %>%
  filter(L1=="Northwestern Forested Mountains")%>%
  select(!L1) %>%
  group_by(class) %>%
  summarise(across(everything(), list(sum))) %>%
  mutate(totalWUI = rowSums(.[3:4]))%>%
  mutate(propWUIinterface = Interface_1 / sum(Interface_1 )) %>%
  mutate(propWUIintermix = Intermix_1 / sum(Intermix_1 )) %>%
  select(class,propWUIinterface, propWUIintermix)
# Simulate response data based on multinomial distribution
# Set simulation parameters   
obs.size = 100 # effective N
propWUI = propWUI %>%
  filter(!class %in% 'Developed')
resp.cats = propWUI$class # category names
# Get observed probabilities but drop categories with < 5% total 
  obs.prob.face <- 
    propWUI %>% 
      select(propWUIinterface) %>% 
        as_vector() 
  obs.prob.mix <- 
    propWUI %>% 
      select(propWUIintermix) %>% 
      as_vector() 
  exp.prob = 1/length(obs.prob.face) # null hypothesis (even % across categories)
  n.sims = 1000  # total interactions for simulation
# Create empty array into which simulations are added for interface
  cat.sims.face <- tibble() 
  for (i in 1:n.sims) { 
  rmultinom(1, size=jitter(obs.size, 10), prob=obs.prob.face) %>%
                      t() %>%
                      as_tibble() %>%
      bind_rows(cat.sims.face) -> cat.sims.face
  }
  names(cat.sims.face)<-(resp.cats)
  cat.sims.face %<>% add_column( type = 'face', .before = 1)
# Create empty array into which simulations are added for intermix
  cat.sims.mix <- tibble() 
  for (i in 1:n.sims) { 
    rmultinom(1, size=jitter(obs.size, 10), prob=obs.prob.mix) %>%
      t() %>%
      as_tibble() %>%
      bind_rows(cat.sims.mix) -> cat.sims.mix
  }
  names(cat.sims.mix)<-(resp.cats)
  cat.sims.mix %<>% add_column( type = 'mix', .before = 1)
# Create empty tibble to stash stats results
  PearsonsOut <- tibble() 
# Run Chi-squared test on simulated data
  for (i in 1:n.sims) {  
    # Grab a row of simulated data for interface 
      sim.d.face <- cat.sims.face %>% slice(i)
    # Grab a row of simulated data for intermix
      sim.d.mix <- cat.sims.mix %>% slice(i)
    # Combine intermix and face into one table
      sim.d <- full_join(by = c('type',"Farmland",'Other', "Rangeland"), 
                         sim.d.face, 
                         sim.d.mix) %>% 
        pivot_longer(!type, names_to = "class", values_to = "count") %>%  
        pivot_wider(names_from=type, values_from=count) 
    # run Chi-squared test
      x2 <- chisq.test(sim.d %>% select(-class))
    # calculate effect size
      ef <- effectsize::effectsize(x2, type = 'pearsons_c', alternative = 'two')
    # extract residuals, combine with effect size, export results of i-th run
     x2= x2$residuals %>% 
        as_tibble()  %>%
        mutate(type = sim.d$class) %>%
        pivot_longer(names_to = 'class', 
                     values_to = 'residual', 
                    -type) %>%
        pivot_wider(names_from = c(type, class), 
                    values_from = residual) %>%
        add_column(PearsonsC = ef$Pearsons_c, .before = 1) %>%
        bind_rows(PearsonsOut) -> PearsonsOut
  } 
# Calculate median and 95% CI from N simulation results 
PearsonsOut %>% 
  reframe(across(c(PearsonsC:Rangeland_mix), ~quantile(.x, prob=c(0.025,0.5,0.975))) ) %>%
  add_column(quantile = c('ciL', 'median', 'ciU'), .before = 1) %>%
  pivot_longer(values_to = 'value', 
               names_to = 'factor', 
               -quantile) %>%
  pivot_wider(names_from = quantile, 
              values_from = value) %>%
  separate(factor, into=c('type', 'class'), sep = '_') %>%
  pander("Pearson's C = effect size for the overall Chi-squared test. Remaining rows are Pearson's contrasts for each factor level.")
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 1 rows [1].
Pearson’s C = effect size for the overall Chi-squared test. Remaining rows are Pearson’s contrasts for each factor level.
type class ciL median ciU
PearsonsC NA 0.474 0.5383 0.5906
Farmland face 3.873 4.933 5.983
Farmland mix -5.678 -4.926 -4.027
Other face -4.344 -3.381 -2.347
Other mix 2.346 3.35 4.49
Rangeland face -3.145 -2.064 -0.9645
Rangeland mix 0.9525 2.053 3.215
#Temperate Sierras

propWUI <-
  DataSummaries$EcoregionRangeWUI %>%
  select(!distance) %>%
  filter(L1=="Temperate Sierras") %>%
  select(!L1) %>%
  replace_na(list(Interface = 0, Intermix=0, NotWUI=0)) %>%
  group_by(class) %>%
  summarise(across(everything(), list(sum))) %>%
  mutate(totalWUI = rowSums(.[3:4]))%>%
  mutate(propWUIinterface = Interface_1 / sum(Interface_1 )) %>%
  mutate(propWUIintermix = Intermix_1 / sum(Intermix_1 )) %>%
  select(class,propWUIinterface, propWUIintermix)
# Simulate response data based on multinomial distribution
# Set simulation parameters   
obs.size = 100 # effective N
propWUI = propWUI %>%
  filter(!class %in% c('Developed', 'Farmland'))
resp.cats = propWUI$class # category names
# Get observed probabilities but drop categories with < 5% total 
  obs.prob.face <- 
    propWUI %>% 
      select(propWUIinterface) %>% 
        as_vector() 
  obs.prob.mix <- 
    propWUI %>% 
      select(propWUIintermix) %>% 
      as_vector() 
  exp.prob = 1/length(obs.prob.face) # null hypothesis (even % across categories)
  n.sims = 1000  # total interactions for simulation
# Create empty array into which simulations are added for interface
  cat.sims.face <- tibble() 
  for (i in 1:n.sims) { 
  rmultinom(1, size=jitter(obs.size, 10), prob=obs.prob.face) %>%
                      t() %>%
                      as_tibble() %>%
      bind_rows(cat.sims.face) -> cat.sims.face
  }
  names(cat.sims.face)<-(resp.cats)
  cat.sims.face %<>% add_column( type = 'face', .before = 1)
# Create empty array into which simulations are added for intermix
  cat.sims.mix <- tibble() 
  for (i in 1:n.sims) { 
    rmultinom(1, size=jitter(obs.size, 10), prob=obs.prob.mix) %>%
      t() %>%
      as_tibble() %>%
      bind_rows(cat.sims.mix) -> cat.sims.mix
  }
  names(cat.sims.mix)<-(resp.cats)
  cat.sims.mix %<>% add_column( type = 'mix', .before = 1)
# Create empty tibble to stash stats results
  PearsonsOut <- tibble() 
# Run Chi-squared test on simulated data
  for (i in 1:n.sims) {  
    # Grab a row of simulated data for interface 
      sim.d.face <- cat.sims.face %>% slice(i)
    # Grab a row of simulated data for intermix
      sim.d.mix <- cat.sims.mix %>% slice(i)
    # Combine intermix and face into one table
      sim.d <- full_join(by = c('type','Other', "Rangeland"), 
                         sim.d.face, 
                         sim.d.mix) %>% 
        pivot_longer(!type, names_to = "class", values_to = "count") %>%  
        pivot_wider(names_from=type, values_from=count) 
    # run Chi-squared test
      x2 <- chisq.test(sim.d %>% select(-class))
    # calculate effect size
      ef <- effectsize::effectsize(x2, type = 'pearsons_c', alternative = 'two')
    # extract residuals, combine with effect size, export results of i-th run
     x2= x2$residuals %>% 
        as_tibble()  %>%
        mutate(type = sim.d$class) %>%
        pivot_longer(names_to = 'class', 
                     values_to = 'residual', 
                    -type) %>%
        pivot_wider(names_from = c(type, class), 
                    values_from = residual) %>%
        add_column(PearsonsC = ef$Pearsons_c, .before = 1) %>%
        bind_rows(PearsonsOut) -> PearsonsOut
  } 
# Calculate median and 95% CI from N simulation results 
PearsonsOut %>% 
  reframe(across(c(PearsonsC:Rangeland_mix), ~quantile(.x, prob=c(0.025,0.5,0.975))) ) %>%
  add_column(quantile = c('ciL', 'median', 'ciU'), .before = 1) %>%
  pivot_longer(values_to = 'value', 
               names_to = 'factor', 
               -quantile) %>%
  pivot_wider(names_from = quantile, 
              values_from = value) %>%
  separate(factor, into=c('type', 'class'), sep = '_') %>%
  pander("Pearson's C = effect size for the overall Chi-squared test. Remaining rows are Pearson's contrasts for each factor level.")
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 1 rows [1].
Pearson’s C = effect size for the overall Chi-squared test. Remaining rows are Pearson’s contrasts for each factor level.
type class ciL median ciU
PearsonsC NA 0.1941 0.3196 0.4235
Other face 1.318 2.371 3.503
Other mix -3.327 -2.397 -1.339
Rangeland face -3.333 -2.344 -1.369
Rangeland mix 1.309 2.349 3.31

Test if the total amount of WUI differs across urban/rural categories - simulation method

###chi squared testing if the proportion of WUI differs across rural/urban

propWUI <-
  DataSummaries$OverallRuralWUI %>%
  mutate(totalWUI = Interface + Intermix, 
         propWUI = totalWUI / sum(totalWUI) )  %>%
  select(RuralUrban, propWUI )
# Simulate response data based on multinomial distribution
# Set simulation parameters   
obs.size = 100 # effective N
resp.cats = propWUI$RuralUrban # category names
# Get observed probabilities but drop categories with < 5% total 
obs.prob <- propWUI %>% 
  filter(propWUI >= 0.05) %>%
  select(propWUI) %>% 
  as_vector() 
exp.prob = 1/length(obs.prob) # null hypothesis (even % across categories)
n.sims = 1000  # total interactions for simulation
# Create empty array into which simulations are added
cat.sims<-array(NA,c(n.sims,length(obs.prob)) )
cat.sims <- rmultinom(n.sims, size=obs.size, prob=obs.prob) %>%
  t() %>%
  as.data.frame() 
names(cat.sims)<-levels(resp.cats)
# Create mpty tibble to stash stats results
PearsonsOut <- tibble() 
# Run Chi-squared test on simulated data
for (i in 1:nrow(cat.sims)){  
  # Grab a row of simulated data  
  sim.d <- cat.sims %>% slice(i)
  # run Chi-squared test
  x2 <- chisq.test(sim.d)
  # calculate effect size
  ef <- effectsize::effectsize(x2, type = 'pearsons_c', alternative = 'two')
  # extract residuals, combine with effect size, export results of i-th run
  x2$residuals %>% 
    as_tibble()  %>%
    mutate(levels = resp.cats) %>%
    pivot_wider(names_from = levels, 
                values_from = value) %>%
    add_column(PearsonsC = ef$Pearsons_c, .before = 1) %>%
    bind_rows(PearsonsOut) -> PearsonsOut
} 
# Calculate median and 95% CI from N simulation results 
PearsonsOut %>%
  reframe(across(c(PearsonsC:Rural), ~quantile(.x, prob=c(0.025,0.5,0.975))) ) %>%
  add_column(quantile = c('ciL', 'median', 'ciU'), .before = 1) %>%
  pivot_longer(values_to = 'value', 
               names_to = 'factor', 
               -quantile) %>%
  pivot_wider(names_from = quantile, 
              values_from = value) %>%
  pander("Pearson's C = effect size for the overall Chi-squared test. Remaining rows are Pearson's contrasts for each factor level.")
Pearson’s C = effect size for the overall Chi-squared test. Remaining rows are Pearson’s contrasts for each factor level. Testing if the proportion of WUI types differs across urban/rural categories.
factor ciL median ciU
PearsonsC 0.4262 0.5371 0.6326
Urban core -2.46 -0.8944 0.6708
Peri-urban 0.8944 2.683 4.919
Small town -3.578 -2.46 -1.118
Rural commuting -4.025 -3.13 -2.012
Rural 1.789 3.578 5.814
 #make contingency table
  wui.dat=DataSummaries$OverallRuralWUI %>%
    mutate_at(vars(Interface:NotWUI), ~ as.integer(round(.x)))

Within ecoregions test for differences in total WUI in different urban/rural classes.

#Great Plains
propWUI <-
  DataSummaries$EcoRuralWUI %>%
   filter(L1=="Great Plains") %>%
  mutate(totalWUI = Interface + Intermix, 
         propWUI = totalWUI / sum(totalWUI) )  %>%
  select(RuralUrban, propWUI )
# Simulate response data based on multinomial distribution
# Set simulation parameters   
obs.size = 100 # effective N
resp.cats = propWUI$RuralUrban # category names
# Get observed probabilities but drop categories with < 5% total 
obs.prob <- propWUI %>% 
  filter(propWUI >= 0.05) %>%
  select(propWUI) %>% 
  as_vector() 
exp.prob = 1/length(obs.prob) # null hypothesis (even % across categories)
n.sims = 1000  # total interactions for simulation
# Create empty array into which simulations are added
cat.sims<-array(NA,c(n.sims,length(obs.prob)) )
cat.sims <- rmultinom(n.sims, size=obs.size, prob=obs.prob) %>%
  t() %>%
  as.data.frame() 
names(cat.sims)<-levels(resp.cats)
# Create mpty tibble to stash stats results
PearsonsOut <- tibble() 
# Run Chi-squared test on simulated data
for (i in 1:nrow(cat.sims)){  
  # Grab a row of simulated data  
  sim.d <- cat.sims %>% slice(i)
  # run Chi-squared test
  x2 <- chisq.test(sim.d)
  # calculate effect size
  ef <- effectsize::effectsize(x2, type = 'pearsons_c', alternative = 'two')
  # extract residuals, combine with effect size, export results of i-th run
  x2$residuals %>% 
    as_tibble()  %>%
    mutate(levels = resp.cats) %>%
    pivot_wider(names_from = levels, 
                values_from = value) %>%
    add_column(PearsonsC = ef$Pearsons_c, .before = 1) %>%
    bind_rows(PearsonsOut) -> PearsonsOut
} 
# Calculate median and 95% CI from N simulation results 
PearsonsOut %>%
  reframe(across(c(PearsonsC:Rural), ~quantile(.x, prob=c(0.025,0.5,0.975))) ) %>%
  add_column(quantile = c('ciL', 'median', 'ciU'), .before = 1) %>%
  pivot_longer(values_to = 'value', 
               names_to = 'factor', 
               -quantile) %>%
  pivot_wider(names_from = quantile, 
              values_from = value) %>%
  pander("Pearson's C = effect size for the overall Chi-squared test. Remaining rows are Pearson's contrasts for each factor level.")
Pearson’s C = effect size for the overall Chi-squared test. Remaining rows are Pearson’s contrasts for each factor level.
factor ciL median ciU
PearsonsC 0.4852 0.5936 0.6763
Urban core -3.13 -1.789 -0.4472
Peri-urban 0.6708 2.683 4.701
Small town -3.578 -2.46 -1.118
Rural commuting -4.249 -3.354 -2.236
Rural 2.683 4.919 7.155
#North American Deserts
propWUI <-
  DataSummaries$EcoRuralWUI %>%
   filter(L1=="North American Deserts") %>%
  mutate(totalWUI = Interface + Intermix, 
         propWUI = totalWUI / sum(totalWUI) )  %>%
  select(RuralUrban, propWUI )
# Simulate response data based on multinomial distribution
# Set simulation parameters   
obs.size = 100 # effective N
resp.cats = propWUI$RuralUrban # category names
# Get observed probabilities but drop categories with < 5% total 
obs.prob <- propWUI %>% 
  filter(propWUI >= 0.05) %>%
  select(propWUI) %>% 
  as_vector() 
exp.prob = 1/length(obs.prob) # null hypothesis (even % across categories)
n.sims = 1000  # total interactions for simulation
# Create empty array into which simulations are added
cat.sims<-array(NA,c(n.sims,length(obs.prob)) )
cat.sims <- rmultinom(n.sims, size=obs.size, prob=obs.prob) %>%
  t() %>%
  as.data.frame() 
names(cat.sims)<-levels(resp.cats)
# Create mpty tibble to stash stats results
PearsonsOut <- tibble() 
# Run Chi-squared test on simulated data
for (i in 1:nrow(cat.sims)){  
  # Grab a row of simulated data  
  sim.d <- cat.sims %>% slice(i)
  # run Chi-squared test
  x2 <- chisq.test(sim.d)
  # calculate effect size
  ef <- effectsize::effectsize(x2, type = 'pearsons_c', alternative = 'two')
  # extract residuals, combine with effect size, export results of i-th run
  x2$residuals %>% 
    as_tibble()  %>%
    mutate(levels = resp.cats) %>%
    pivot_wider(names_from = levels, 
                values_from = value) %>%
    add_column(PearsonsC = ef$Pearsons_c, .before = 1) %>%
    bind_rows(PearsonsOut) -> PearsonsOut
} 
# Calculate median and 95% CI from N simulation results 
PearsonsOut %>%
  reframe(across(c(PearsonsC:Rural), ~quantile(.x, prob=c(0.025,0.5,0.975))) ) %>%
  add_column(quantile = c('ciL', 'median', 'ciU'), .before = 1) %>%
  pivot_longer(values_to = 'value', 
               names_to = 'factor', 
               -quantile) %>%
  pivot_wider(names_from = quantile, 
              values_from = value) %>%
  pander("Pearson's C = effect size for the overall Chi-squared test. Remaining rows are Pearson's contrasts for each factor level.")
Pearson’s C = effect size for the overall Chi-squared test. Remaining rows are Pearson’s contrasts for each factor level.
factor ciL median ciU
PearsonsC 0.3887 0.4987 0.58
Urban core -0.6708 1.118 3.13
Peri-urban 0.8944 2.683 4.919
Small town -3.578 -2.46 -1.118
Rural commuting -4.249 -3.354 -2.236
Rural -0.2236 1.789 3.801
#Northwestern Forested Mountains
propWUI <-
  DataSummaries$EcoRuralWUI %>%
   filter(L1=="Northwestern Forested Mountains") %>%
  mutate(totalWUI = Interface + Intermix, 
         propWUI = totalWUI / sum(totalWUI) )  %>%
  select(RuralUrban, propWUI )
# Simulate response data based on multinomial distribution
# Set simulation parameters   
obs.size = 100 # effective N
resp.cats = propWUI$RuralUrban # category names
# Get observed probabilities but drop categories with < 5% total 
obs.prob <- propWUI %>% 
  filter(propWUI >= 0.05) %>%
  select(propWUI) %>% 
  as_vector() 
exp.prob = 1/length(obs.prob) # null hypothesis (even % across categories)
n.sims = 1000  # total interactions for simulation
# Create empty array into which simulations are added
cat.sims<-array(NA,c(n.sims,length(obs.prob)) )
cat.sims <- rmultinom(n.sims, size=obs.size, prob=obs.prob) %>%
  t() %>%
  as.data.frame() 
names(cat.sims)<-levels(resp.cats)
# Create mpty tibble to stash stats results
PearsonsOut <- tibble() 
# Run Chi-squared test on simulated data
for (i in 1:nrow(cat.sims)){  
  # Grab a row of simulated data  
  sim.d <- cat.sims %>% slice(i)
  # run Chi-squared test
  x2 <- chisq.test(sim.d)
  # calculate effect size
  ef <- effectsize::effectsize(x2, type = 'pearsons_c', alternative = 'two')
  # extract residuals, combine with effect size, export results of i-th run
  x2$residuals %>% 
    as_tibble()  %>%
    mutate(levels = resp.cats) %>%
    pivot_wider(names_from = levels, 
                values_from = value) %>%
    add_column(PearsonsC = ef$Pearsons_c, .before = 1) %>%
    bind_rows(PearsonsOut) -> PearsonsOut
} 
# Calculate median and 95% CI from N simulation results 
PearsonsOut %>%
  reframe(across(c(PearsonsC:Rural), ~quantile(.x, prob=c(0.025,0.5,0.975))) ) %>%
  add_column(quantile = c('ciL', 'median', 'ciU'), .before = 1) %>%
  pivot_longer(values_to = 'value', 
               names_to = 'factor', 
               -quantile) %>%
  pivot_wider(names_from = quantile, 
              values_from = value) %>%
  pander("Pearson's C = effect size for the overall Chi-squared test. Remaining rows are Pearson's contrasts for each factor level.")
Pearson’s C = effect size for the overall Chi-squared test. Remaining rows are Pearson’s contrasts for each factor level.
factor ciL median ciU
PearsonsC 0.4604 0.5642 0.6483
Urban core -3.13 -1.789 -0.2236
Peri-urban 1.118 2.907 5.143
Small town -3.801 -2.907 -1.565
Rural commuting -3.801 -2.683 -1.342
Rural 2.012 4.025 6.261
#Southern Semi-Arid Highlands
propWUI <-
  DataSummaries$EcoRuralWUI %>%
   filter(L1=="Southern Semi-Arid Highlands") %>%
  mutate(totalWUI = Interface + Intermix, 
         propWUI = totalWUI / sum(totalWUI) )  %>%
  select(RuralUrban, propWUI )
# Simulate response data based on multinomial distribution
# Set simulation parameters   
obs.size = 100 # effective N
propWUI = propWUI %>%
  filter(!RuralUrban %in% 'Small town') 
resp.cats = propWUI$RuralUrban 

# Get observed probabilities but drop categories with < 5% total 
obs.prob <- propWUI %>% 
  filter(propWUI >= 0.05) %>%
  select(propWUI) %>% 
  as_vector() 
exp.prob = 1/length(obs.prob) # null hypothesis (even % across categories)
n.sims = 1000  # total interactions for simulation
# Create empty array into which simulations are added
cat.sims<-array(NA,c(n.sims,length(obs.prob)) )
cat.sims <- rmultinom(n.sims, size=obs.size, prob=obs.prob) %>%
  t() %>%
  as.data.frame() 
names(cat.sims)<-resp.cats
# Create mpty tibble to stash stats results
PearsonsOut <- tibble() 
# Run Chi-squared test on simulated data
for (i in 1:nrow(cat.sims)){  
  # Grab a row of simulated data  
  sim.d <- cat.sims %>% slice(i)
  # run Chi-squared test
  x2 <- chisq.test(sim.d)
  # calculate effect size
  ef <- effectsize::effectsize(x2, type = 'pearsons_c', alternative = 'two')
  # extract residuals, combine with effect size, export results of i-th run
  x2$residuals %>% 
    as_tibble()  %>%
    mutate(levels = resp.cats) %>%
    pivot_wider(names_from = levels, 
                values_from = value) %>%
    add_column(PearsonsC = ef$Pearsons_c, .before = 1) %>%
    bind_rows(PearsonsOut) -> PearsonsOut
} 
# Calculate median and 95% CI from N simulation results 
PearsonsOut %>%
  reframe(across(c(PearsonsC:Rural), ~quantile(.x, prob=c(0.025,0.5,0.975))) ) %>%
  add_column(quantile = c('ciL', 'median', 'ciU'), .before = 1) %>%
  pivot_longer(values_to = 'value', 
               names_to = 'factor', 
               -quantile) %>%
  pivot_wider(names_from = quantile, 
              values_from = value) %>%
  pander("Pearson's C = effect size for the overall Chi-squared test. Remaining rows are Pearson's contrasts for each factor level.")
Pearson’s C = effect size for the overall Chi-squared test. Remaining rows are Pearson’s contrasts for each factor level.
factor ciL median ciU
PearsonsC 0.3273 0.4417 0.5496
Urban core -1.8 -0.3 1.4
Peri-urban 1.2 3.2 5.2
Rural commuting -4.4 -3.4 -2.4
Rural -1.2 0.6 2.4
#Temperate Sierras
propWUI <-
  DataSummaries$EcoRuralWUI %>%
   filter(L1=="Temperate Sierras") %>%
  mutate(totalWUI = Interface + Intermix, 
         propWUI = totalWUI / sum(totalWUI) )  %>%
  select(RuralUrban, propWUI )
# Simulate response data based on multinomial distribution
# Set simulation parameters   
obs.size = 100 # effective N
propWUI = propWUI %>%
  filter(!RuralUrban %in% 'Rural commuting') 
resp.cats = propWUI$RuralUrban  # category names
# Get observed probabilities but drop categories with < 5% total 
obs.prob <- propWUI %>% 
  filter(propWUI >= 0.05) %>%
  select(propWUI) %>% 
  as_vector() 
exp.prob = 1/length(obs.prob) # null hypothesis (even % across categories)
n.sims = 1000  # total interactions for simulation
# Create empty array into which simulations are added
cat.sims<-array(NA,c(n.sims,length(obs.prob)) )
cat.sims <- rmultinom(n.sims, size=obs.size, prob=obs.prob) %>%
  t() %>%
  as.data.frame() 
names(cat.sims)<-resp.cats
# Create mpty tibble to stash stats results
PearsonsOut <- tibble() 
# Run Chi-squared test on simulated data
for (i in 1:nrow(cat.sims)){  
  # Grab a row of simulated data  
  sim.d <- cat.sims %>% slice(i)
  # run Chi-squared test
  x2 <- chisq.test(sim.d)
  # calculate effect size
  ef <- effectsize::effectsize(x2, type = 'pearsons_c', alternative = 'two')
  # extract residuals, combine with effect size, export results of i-th run
  x2$residuals %>% 
    as_tibble()  %>%
    mutate(levels = resp.cats) %>%
    pivot_wider(names_from = levels, 
                values_from = value) %>%
    add_column(PearsonsC = ef$Pearsons_c, .before = 1) %>%
    bind_rows(PearsonsOut) -> PearsonsOut
} 
# Calculate median and 95% CI from N simulation results 
PearsonsOut %>%
  reframe(across(c(PearsonsC:Rural), ~quantile(.x, prob=c(0.025,0.5,0.975))) ) %>%
  add_column(quantile = c('ciL', 'median', 'ciU'), .before = 1) %>%
  pivot_longer(values_to = 'value', 
               names_to = 'factor', 
               -quantile) %>%
  pivot_wider(names_from = quantile, 
              values_from = value) %>%
  pander("Pearson's C = effect size for the overall Chi-squared test. Remaining rows are Pearson's contrasts for each factor level.")
Pearson’s C = effect size for the overall Chi-squared test. Remaining rows are Pearson’s contrasts for each factor level. Test to see if the amount of WUI differs across woody encroachment classes for all interior west
factor ciL median ciU
PearsonsC 0.2106 0.3669 0.4938
Urban core -1.2 0.4 2.4
Peri-urban 0.8 2.6 4.6
Small town -3.6 -2.4 -1
Rural -2.4 -0.8 1
propWUI <-
  DataSummaries$OverallWoodyWUI %>%
  filter(scope=="All Interior West")%>%
  select(class:NotWUI) %>%
  mutate(propWUI = WUI / sum(WUI)) 
 
# Simulate response data based on multinomial distribution
# Set simulation parameters   
obs.size = 100 # effective N
propWUI = propWUI %>%
  filter(!class %in% 'Tree-encroached') 
resp.cats = propWUI$class # category names
# Get observed probabilities but drop categories with < 5% total 
obs.prob <- propWUI %>% 
  filter(propWUI >= 0.05) %>%
  select(propWUI) %>% 
  as_vector() 
exp.prob = 1/length(obs.prob.face) # null hypothesis (even % across categories)
n.sims = 1000  # total interactions for simulation
# Create empty array into which simulations are added
cat.sims<-array(NA,c(n.sims,length(obs.prob)) )
cat.sims <- rmultinom(n.sims, size=obs.size, prob=obs.prob) %>%
  t() %>%
  as.data.frame() 

names(cat.sims)<-resp.cats
# Create mpty tibble to stash stats results
PearsonsOut <- tibble() 
# Run Chi-squared test on simulated data
for (i in 1:nrow(cat.sims)){  
  # Grab a row of simulated data  
  sim.d <- cat.sims %>% slice(i)
  # run Chi-squared test
  x2 <- chisq.test(sim.d)
  # calculate effect size
  ef <- effectsize::effectsize(x2, type = 'pearsons_c', alternative = 'two')
  # extract residuals, combine with effect size, export results of i-th run
  x2$residuals %>% 
    as_tibble()  %>%
    mutate(levels = resp.cats) %>%
    pivot_wider(names_from = levels, 
                values_from = value) %>%
    add_column(PearsonsC = ef$Pearsons_c, .before = 1) %>%
    bind_rows(PearsonsOut) -> PearsonsOut
} 
# Calculate median and 95% CI from N simulation results 
PearsonsOut %>%
  reframe(across(c(PearsonsC:Rangeland), ~quantile(.x, prob=c(0.025,0.5,0.975))) ) %>%
  add_column(quantile = c('ciL', 'median', 'ciU'), .before = 1) %>%
  pivot_longer(values_to = 'value', 
               names_to = 'factor', 
               -quantile) %>%
  pivot_wider(names_from = quantile, 
              values_from = value) %>%
  pander("Pearson's C = effect size for the overall Chi-squared test. Remaining rows are Pearson's contrasts for each factor level.") 
Pearson’s C = effect size for the overall Chi-squared test. Remaining rows are Pearson’s contrasts for each factor level. Test to see if the amount of WUI differs across woody encroachment classes for tribal areas
factor ciL median ciU
PearsonsC 0.6051 0.6606 0.6925
Rangeland 5.374 6.223 6.788
#####
propWUI <-
  DataSummaries$OverallWoodyWUI %>%
  filter(scope=="Tribal areas only")%>%
  select(class:NotWUI) %>%
  mutate(propWUI = WUI / sum(WUI)) 

Within ecoregions test for differences in total WUI in different woody encroachment classes for all interior west

#Great Plains
propWUI <-
  DataSummaries$EcoRegionWoodyWUI %>%
  filter(L1=="Great Plains")%>%
  select(class:NotWUI) %>%
  mutate(propWUI = WUI / sum(WUI)) 
 
# Simulate response data based on multinomial distribution
# Set simulation parameters   
obs.size = 100 # effective N
propWUI = propWUI %>%
  filter(!class %in% 'Tree-encroached') 
resp.cats = propWUI$class # category names
# Get observed probabilities but drop categories with < 5% total 
obs.prob <- propWUI %>% 
  filter(propWUI >= 0.05) %>%
  select(propWUI) %>% 
  as_vector() 
exp.prob = 1/length(obs.prob.face) # null hypothesis (even % across categories)
n.sims = 1000  # total interactions for simulation
# Create empty array into which simulations are added
cat.sims<-array(NA,c(n.sims,length(obs.prob)) )
cat.sims <- rmultinom(n.sims, size=obs.size, prob=obs.prob) %>%
  t() %>%
  as.data.frame() 

names(cat.sims)<-resp.cats
# Create mpty tibble to stash stats results
PearsonsOut <- tibble() 
# Run Chi-squared test on simulated data
for (i in 1:nrow(cat.sims)){  
  # Grab a row of simulated data  
  sim.d <- cat.sims %>% slice(i)
  # run Chi-squared test
  x2 <- chisq.test(sim.d)
  # calculate effect size
  ef <- effectsize::effectsize(x2, type = 'pearsons_c', alternative = 'two')
  # extract residuals, combine with effect size, export results of i-th run
  x2$residuals %>% 
    as_tibble()  %>%
    mutate(levels = resp.cats) %>%
    pivot_wider(names_from = levels, 
                values_from = value) %>%
    add_column(PearsonsC = ef$Pearsons_c, .before = 1) %>%
    bind_rows(PearsonsOut) -> PearsonsOut
} 
# Calculate median and 95% CI from N simulation results 
PearsonsOut %>%
  reframe(across(c(PearsonsC:Rangeland), ~quantile(.x, prob=c(0.025,0.5,0.975))) ) %>%
  add_column(quantile = c('ciL', 'median', 'ciU'), .before = 1) %>%
  pivot_longer(values_to = 'value', 
               names_to = 'factor', 
               -quantile) %>%
  pivot_wider(names_from = quantile, 
              values_from = value) %>%
  pander("Pearson's C = effect size for the overall Chi-squared test. Remaining rows are Pearson's contrasts for each factor level.") 
Pearson’s C = effect size for the overall Chi-squared test. Remaining rows are Pearson’s contrasts for each factor level.
factor ciL median ciU
PearsonsC 0.615 0.669 0.6999
Rangeland 5.515 6.364 6.93
  #North American Deserts
 propWUI <-
  DataSummaries$EcoRegionWoodyWUI %>%
  filter(L1=="North American Deserts")%>%
  select(class:NotWUI) %>%
  mutate(propWUI = WUI / sum(WUI)) 
 

  #Northwestern Forested Mountains
 propWUI <-
  DataSummaries$EcoRegionWoodyWUI %>%
  filter(L1=="Northwestern Forested Mountains")%>%
  select(class:NotWUI) %>%
  mutate(propWUI = WUI / sum(WUI)) 
 
# Simulate response data based on multinomial distribution
# Set simulation parameters   
obs.size = 100 # effective N
propWUI = propWUI %>%
  filter(!class %in% 'Tree-encroached') 
resp.cats = propWUI$class # category names
# Get observed probabilities but drop categories with < 5% total 
obs.prob <- propWUI %>% 
  filter(propWUI >= 0.05) %>%
  select(propWUI) %>% 
  as_vector() 
exp.prob = 1/length(obs.prob.face) # null hypothesis (even % across categories)
n.sims = 1000  # total interactions for simulation
# Create empty array into which simulations are added
cat.sims<-array(NA,c(n.sims,length(obs.prob)) )
cat.sims <- rmultinom(n.sims, size=obs.size, prob=obs.prob) %>%
  t() %>%
  as.data.frame() 

names(cat.sims)<-resp.cats
# Create mpty tibble to stash stats results
PearsonsOut <- tibble() 
# Run Chi-squared test on simulated data
for (i in 1:nrow(cat.sims)){  
  # Grab a row of simulated data  
  sim.d <- cat.sims %>% slice(i)
  # run Chi-squared test
  x2 <- chisq.test(sim.d)
  # calculate effect size
  ef <- effectsize::effectsize(x2, type = 'pearsons_c', alternative = 'two')
  # extract residuals, combine with effect size, export results of i-th run
  x2$residuals %>% 
    as_tibble()  %>%
    mutate(levels = resp.cats) %>%
    pivot_wider(names_from = levels, 
                values_from = value) %>%
    add_column(PearsonsC = ef$Pearsons_c, .before = 1) %>%
    bind_rows(PearsonsOut) -> PearsonsOut
} 
# Calculate median and 95% CI from N simulation results 
PearsonsOut %>%
  reframe(across(c(PearsonsC:Rangeland), ~quantile(.x, prob=c(0.025,0.5,0.975))) ) %>%
  add_column(quantile = c('ciL', 'median', 'ciU'), .before = 1) %>%
  pivot_longer(values_to = 'value', 
               names_to = 'factor', 
               -quantile) %>%
  pivot_wider(names_from = quantile, 
              values_from = value) %>%
  pander("Pearson's C = effect size for the overall Chi-squared test. Remaining rows are Pearson's contrasts for each factor level.") 
Pearson’s C = effect size for the overall Chi-squared test. Remaining rows are Pearson’s contrasts for each factor level.
factor ciL median ciU
PearsonsC 0.4751 0.5623 0.6341
Rangeland 3.818 4.808 5.798
 #Temperate Sierras 
   propWUI <-
  DataSummaries$EcoRegionWoodyWUI %>%
  filter(L1=="Temperate Sierras")%>%
  select(class:NotWUI) %>%
  mutate(propWUI = WUI / sum(WUI)) 
 


  #Southern Semi-Arid Highlands 
   propWUI <-
  DataSummaries$EcoRegionWoodyWUI %>%
  filter(L1=="Southern Semi-Arid Highlands")%>%
  select(class:NotWUI) %>%
  mutate(propWUI = WUI / sum(WUI)) 
propWUI <-
  DataSummaries$OverallRangeWUI %>%
  select(!class) %>%
  group_by(distance) %>%
  summarise(across(everything(), list(sum))) %>%
  mutate(totalWUI = rowSums(.[3:4]))%>%
  mutate(propWUI = totalWUI / sum(totalWUI)) %>%
  select(!NotWUI_1:totalWUI)
# Simulate response data based on multinomial distribution
# Set simulation parameters   
obs.size = 100 # effective N
propWUI = propWUI %>%
  filter(!distance %in% '75+') 
resp.cats = propWUI$distance # category names
# Get observed probabilities but drop categories with < 5% total 
obs.prob <- propWUI %>% 
  filter(propWUI >= 0.05) %>%
  select(propWUI) %>% 
  as_vector() 
exp.prob = 1/length(obs.prob) # null hypothesis (even % across categories)
n.sims = 1000  # total interactions for simulation
# Create empty array into which simulations are added
cat.sims<-array(NA,c(n.sims,length(obs.prob)))
cat.sims <- rmultinom(n.sims, size=obs.size, prob=obs.prob) %>%
  t() %>%
  as.data.frame()
names(cat.sims)<-(resp.cats)
# Create mpty tibble to stash stats results
PearsonsOut <- tibble() 
# Run Chi-squared test on simulated data
for (i in 1:nrow(cat.sims)){  
  # Grab a row of simulated data  
  sim.d <- cat.sims %>% slice(i)
  # run Chi-squared test
  x2 <- chisq.test(sim.d)
  # calculate effect size
  ef <- effectsize::effectsize(x2, type = 'pearsons_c', alternative = 'two')
  # extract residuals, combine with effect size, export results of i-th run
  x2$residuals %>% 
    as_tibble()  %>%
    mutate(levels = resp.cats) %>%
    pivot_wider(names_from = levels, 
                values_from = value) %>%
    add_column(PearsonsC = ef$Pearsons_c, .before = 1) %>%
    bind_rows(PearsonsOut) -> PearsonsOut
} 
# Calculate median and 95% CI from N simulation results 
PearsonsOut %>%
  reframe(across(c(PearsonsC:'50-75'), ~quantile(.x, prob=c(0.025,0.5,0.975))) ) %>%
  add_column(quantile = c('ciL', 'median', 'ciU'), .before = 1) %>%
  pivot_longer(values_to = 'value', 
               names_to = 'factor', 
               -quantile) %>%
  pivot_wider(names_from = quantile, 
              values_from = value) %>%
  pander("Pearson's C = effect size for the overall Chi-squared test. Remaining rows are Pearson's contrasts for each factor level.")
Pearson’s C = effect size for the overall Chi-squared test. Remaining rows are Pearson’s contrasts for each factor level.
factor ciL median ciU
PearsonsC 0.3037 0.4382 0.5586
0-10 1.6 3.4 5.405
10-25 -1.405 0.2 2
25-50 -2.2 -0.6 1
50-75 -4.2 -3.2 -1.8