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
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
data <- read.csv('https://raw.githubusercontent.com/jmhulbert/redhot/main/data/urban-data-modified.csv')
# 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"
data <- data %>% filter(Area=="King County") %>% droplevels()
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
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)
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
binomial.af <- glmmTMB(binary.tree.canopy.symptoms ~ DN_AF1,data=data,family=binomial)
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"
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(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"
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()`).
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
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()`).
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(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"
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 | ** |