Welcome Data Analyses Discussion Maps
Analyses
Urban Heat All Areas Portland Tacoma King County iNaturalist City Inventories
library(tidyverse)
library(corrplot)
library(knitr)
library(kableExtra)
library(gghalves)
library(patchwork)
library(scales)
library(glmmTMB)
library(bbmle)
library(DHARMa)

#KING COUNTY

Purpose

The purpose of this page is to summarize an investigation of tree health categories as a binary response to urban heat (temperatures) for trees within King County.

Response variables Binary dead top: no top dieback, top dieback

Read Data

data <- read.csv('https://raw.githubusercontent.com/jmhulbert/redhot/main/data/urban-data-modified.csv')

Wrangle Data

Remove dead trees

# Remove Dead Trees
data.w.dead <- data 
data <- data %>% filter(field.tree.canopy.symptoms!="Tree is dead") %>% filter(field.dieback.percent<100) %>% droplevels()
data$binary.tree.canopy.symptoms <- as.factor(data$binary.tree.canopy.symptoms)
levels(data$binary.tree.canopy.symptoms)
## [1] "Healthy"   "Unhealthy"

Filter to King County Data

data <- data %>% filter(Area=="King County") %>% droplevels()

Visualize King County Data

ggplot(data,aes(binary.tree.canopy.symptoms,DN_AF1))+geom_violin()+coord_flip()+theme_bw()

p1 <- ggplot(data,aes(binary.tree.canopy.symptoms))+geom_histogram(stat="count")+theme_bw()+coord_flip()
## Warning in geom_histogram(stat = "count"): Ignoring unknown parameters:
## `binwidth`, `bins`, and `pad`
p2 <- ggplot(data,aes(reclassified.tree.canopy.symptoms))+geom_histogram(stat="count")+theme_bw()+coord_flip()
## Warning in geom_histogram(stat = "count"): Ignoring unknown parameters:
## `binwidth`, `bins`, and `pad`
p1 + p2

p3 <-ggplot(data,aes(binary.tree.canopy.symptoms,DN_AF1, fill=binary.tree.canopy.symptoms ))+geom_boxplot(alpha=0.5)+theme_bw()+coord_flip() + scale_fill_manual(name="Tree Condition",values=c("#7fcdbb","#fe9929"))+guides(fill=FALSE) +labs(x="Tree Condition",y=NULL)
## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
p4 <-ggplot(data,aes(DN_AF1,fill=binary.tree.canopy.symptoms))+geom_density(alpha=0.5) + scale_fill_manual(name="Tree Condition",values=c("#7fcdbb","#fe9929")) +theme_bw() +labs(x="Afternoon Temperature") +guides(fill=FALSE)
p3 / p4

Analyses

Time Series Model Fit

data$top.dieback <- as.factor(data$top.dieback)
data$thinning <- as.factor(data$thinning)
binomial.daily <- glmmTMB(binary.tree.canopy.symptoms ~ mean.temp.daily ,family=binomial,data=data)
binomial.am <- glmmTMB(binary.tree.canopy.symptoms ~ DN_AM1 ,data=data,family=binomial)
binomial.af <- glmmTMB(binary.tree.canopy.symptoms ~ DN_AF1 ,data=data,family=binomial)
binomial.pm <- glmmTMB(binary.tree.canopy.symptoms ~ DN_PM1,data=data,family=binomial)

Best Fit Comparison

AICtab(binomial.daily,binomial.am,binomial.af,binomial.pm)
##                dAIC df
## binomial.af     0.0 2 
## binomial.am     8.9 2 
## binomial.pm    11.4 2 
## binomial.daily 11.7 2

Afternoon temperature was the best fit model

Afternoon Temperatures

Response 1.0: healthy / unhealthy

binomial.af <- glmmTMB(binary.tree.canopy.symptoms ~ DN_AF1,data=data,family=binomial)

Summary

summary(binomial.af)
##  Family: binomial  ( logit )
## Formula:          binary.tree.canopy.symptoms ~ DN_AF1
## Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##    568.1    576.1   -282.0    564.1      418 
## 
## 
## Conditional model:
##              Estimate Std. Error z value Pr(>|z|)   
## (Intercept) -25.95480    8.10408  -3.203  0.00136 **
## DN_AF1        0.28792    0.09074   3.173  0.00151 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
levels(data$binary.tree.canopy.symptoms)
## [1] "Healthy"   "Unhealthy"

Response 2.0: healthy / top-dieback

data.top <- data %>% filter(reclassified.tree.canopy.symptoms=="Healthy"|reclassified.tree.canopy.symptoms=="Dead Top") %>% droplevels()
top.dieback.binomial.af <- glmmTMB(top.dieback ~ DN_AF1,data=data.top,family=binomial)

Summary

summary(top.dieback.binomial.af)
##  Family: binomial  ( logit )
## Formula:          top.dieback ~ DN_AF1
## Data: data.top
## 
##      AIC      BIC   logLik deviance df.resid 
##    296.7    304.1   -146.4    292.7      297 
## 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -51.4452    14.5502  -3.536 0.000407 ***
## DN_AF1        0.5607     0.1625   3.451 0.000558 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
levels(data$top.dieback)
## [1] "No"  "Yes"

Visualize effects of predictors

First, generate predictions from model

# Create a data frame for predictions
new.data.2 <- expand.grid(
  DN_AF1 = seq(min(80), max(150), length.out = 1000)  # keeping predictor2 constant at its mean
)
new.data.2$predicted_response <- predict(top.dieback.binomial.af, newdata = new.data.2, type = "response")
# Create a data frame for predictions
new.data.2.inf <- expand.grid(
  DN_AF1 = seq(min(80), max(150), length.out = 1000)  # keeping predictor2 constant at its mean
)
new.data.2.inf$predicted_response <- predict(top.dieback.binomial.af, newdata = new.data.2.inf, type = "response")
head(new.data.2)
##     DN_AF1 predicted_response
## 1 80.00000        0.001371565
## 2 80.07007        0.001426445
## 3 80.14014        0.001483516
## 4 80.21021        0.001542868
## 5 80.28028        0.001604590
## 6 80.35035        0.001668777
p <- ggplot(new.data.2.inf, aes(x = DN_AF1, y = predicted_response)) + geom_line() + labs(x = "Afternoon temperature (F)", y = "Dieback Probability") + scale_x_continuous(limits = c(80, 115)) +theme_minimal()

p
## Warning: Removed 500 rows containing missing values or values outside the scale range
## (`geom_line()`).

Calculate Dieback Probability Categories

head(new.data.2[new.data.2$DN_AF1<85.5 & new.data.2$DN_AF1>84.75,])
##      DN_AF1 predicted_response
## 69 84.76476         0.01947612
## 70 84.83483         0.02024072
## 71 84.90490         0.02103468
## 72 84.97497         0.02185910
## 73 85.04505         0.02271508
## 74 85.11512         0.02360378

King County Trees at 85 degrees F had about a 0.023 probability

head(new.data.2[new.data.2$predicted_response>0.024 & new.data.2$predicted_response<0.026,])
##      DN_AF1 predicted_response
## 75 85.18519         0.02452636
## 76 85.25526         0.02548407

King County trees at 85.22 had about a 0.025 probablilty

head(new.data.2[new.data.2$predicted_response>0.048 & new.data.2$predicted_response<0.052,])
##      DN_AF1 predicted_response
## 93 86.44645         0.04852177
## 94 86.51652         0.05036807

King County trees at 86.5 had a 0.05 probability

head(new.data.2[new.data.2$predicted_response>0.095 & new.data.2$predicted_response<0.105,,])
##       DN_AF1 predicted_response
## 112 87.77778         0.09712838
## 113 87.84785         0.10062860
## 114 87.91792         0.10424040

trees at about 87.8 degrees had a 0.1 predicted response.

head(new.data.2[new.data.2$predicted_response>0.24,])
##       DN_AF1 predicted_response
## 140 89.73974          0.2442526
## 141 89.80981          0.2515774
## 142 89.87988          0.2590466
## 143 89.94995          0.2666585
## 144 90.02002          0.2744113
## 145 90.09009          0.2823027

King County Trees at 89.78 had about a 0.25 probability

  • King County Top-Dieback Probability Categories (10%)
    • Most Concern (>10%): >87.8
    • Moderate Concern (2.5-10%): 85.22-87.8
    • Lowest Concern (<2.5%): <85.22
  • King County Top-Dieback Probability Categories (25%)
    • Most Concern (>25%): 89.78
    • Moderate Concern (5-25%): 86.5 - 89.78
    • Lowest Concern (<5%): <86.5
p <- ggplot(new.data.2, aes(x = DN_AF1, y = predicted_response)) + geom_line() + labs(x = "Afternoon temperature (F)", y = "Dieback Probability") + scale_x_continuous(limits = c(80, 92.5),breaks=c(80,82.5,85,87.5,90,92.5)) + scale_y_continuous(limits = c(0, .3)) +theme_minimal()

p2 <- p + annotate(geom="rect",ymin = -Inf, ymax = Inf, 
              xmin = -Inf, xmax = 85.22, fill = alpha('#DDA0DD',0.5)) +
    annotate(geom="rect",ymin = -Inf, ymax = Inf, 
              xmin = 85.22, xmax = 87.8, fill = alpha('#7fcdbb',0.5)) + 
    annotate(geom="rect",ymin = -Inf, ymax = Inf, 
          xmin = 87.8, xmax = Inf, fill = alpha('#fe9929',0.5))
p3 <- p2 + annotate("text", x=82, y=.15, label= "Lowest\nConcern") + annotate("text", x=86.44, y=.15, label= "Moderate\nConcern") + annotate("text", x=91, y=.15, label= "Most\nConcern")
p3
## Warning: Removed 853 rows containing missing values or values outside the scale range
## (`geom_line()`).

The probability of a tree having top dieback increases with afternoon temperature

p <- ggplot(data.top, aes(x=DN_AF1))+geom_histogram() +theme_minimal() +labs(y="Number of Trees",x="Afternoon Temperature (F)") + scale_x_continuous(limits = c(80, 92.5),breaks=c(80,82.5,85,87.5,90,92.5))

p2 <- p + annotate(geom="rect",ymin = -Inf, ymax = Inf, 
              xmin = -Inf, xmax = 85.22, fill = alpha('#DDA0DD',0.5)) +
    annotate(geom="rect",ymin = -Inf, ymax = Inf, 
              xmin = 85.22, xmax = 87.8, fill = alpha('#7fcdbb',0.5)) + 
    annotate(geom="rect",ymin = -Inf, ymax = Inf, 
          xmin = 87.8, xmax = Inf, fill = alpha('#fe9929',0.5))
p3 <- p2 + annotate("text", x=82, y=30, label= "Lowest\nConcern") + annotate("text", x=86.44, y=30, label= "Moderate\nConcern") + annotate("text", x=91, y=30, label= "Most\nConcern")
p3
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_bar()`).

Response 3.0: healthy / thinning

data.thin <- data %>% filter(reclassified.tree.canopy.symptoms=="Healthy"|reclassified.tree.canopy.symptoms=="Thinning Canopy") %>% droplevels()
thinning.binomial.af <- glmmTMB(thinning ~ DN_AF1,data=data.thin,family=binomial)

Summary

summary(thinning.binomial.af)
##  Family: binomial  ( logit )
## Formula:          thinning ~ DN_AF1
## Data: data.thin
## 
##      AIC      BIC   logLik deviance df.resid 
##    302.2    309.6   -149.1    298.2      297 
## 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept) -38.9374    13.3769  -2.911  0.00361 **
## DN_AF1        0.4211     0.1495   2.817  0.00485 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
levels(data$thinning)
## [1] "No"  "Yes"

Overall Summary

The probability of a tree having top dieback increases with afternoon temperature

Area Response Explanatory df.resid Estimate Std.Error z-val Pr(>z) Signif.
King County healthy/unhealthy DN_AF 418 0.28792 0.09074 3.173 0.00151 **
King County no/dead-top DN_AF 297 0.5607 0.1625 3.451 0.000558 ***
King County no/thinning DN_AF 297 0.4211 0.1495 2.817 0.00485 **