Compendium sections
- Title page
- Data preparation and
formatting
- The influence of temperatures on
P. xanthii conidia germination
- The influence of VPD P. xanthii conidia germination and germtube formation
Load R libraries required for analysis
library(data.table) # data structure, cleaning and formatting
library(mgcv) # gam statistical models
library(rstanarm) # Bayesian statistical models
library(ggplot2) # graphing
library(MASS) # logistical ordinal regression functions, among other things
library(splines) # non-linear predictors
library(flextable)
library(marginaleffects) # produce predicted probabilities for polr
source("R/p_star.R")
Impact of vapour pressure deficit on P. xanthii germination
Read in data
g_vpd <- fread("cache/germination_vpd.csv")
Peek at the data
str(g_vpd)
Classes 'data.table' and 'data.frame': 288 obs. of 7 variables:
$ RH : int 99 99 99 99 99 99 99 99 99 99 ...
$ Tm : int 28 28 28 28 28 28 28 28 28 28 ...
$ vpd : num 0.038 0.038 0.038 0.038 0.038 0.038 0.038 0.038 0.038 0.038 ...
$ i_period : int 12 12 12 12 24 24 24 24 36 36 ...
$ rep : int 1 2 3 4 1 2 3 4 1 2 ...
$ germ_n : int 56 54 52 50 64 65 60 63 83 80 ...
$ non_germ_n: int 44 46 48 50 36 35 40 37 17 20 ...
- attr(*, ".internal.selfref")=<externalptr> # inspect environmental treatments
unique(g_vpd[,.(RH,Tm,vpd)])
RH Tm vpd
1: 99 28 0.038
2: 95 28 0.189
3: 85 28 0.567
4: 75 28 0.945
5: 55 28 1.701
6: 32 28 2.570
7: 99 25 0.032
8: 95 25 0.158
9: 85 25 0.475
10: 75 25 0.792
11: 55 25 1.425
12: 32 25 2.154
13: 99 22 0.026
14: 95 22 0.132
15: 85 22 0.397
16: 75 22 0.661
17: 55 22 1.190
18: 32 22 1.798This table shows there are 18 VPD treatments consisting of six relative humidity and three temperatures.
Sampling times in hours.
i_period reps
1: 12 1 - 4
2: 24 1 - 4
3: 36 1 - 4
4: 48 1 - 4Visualise the raw data
Let’s plot the germination over time split by relative humidity.
g_vpd |>
ggplot(aes(x = i_period, y = germ_n, colour = as.factor(Tm))) +
geom_point()+
facet_wrap(facets = "RH")+
geom_smooth(method = "gam",formula = y ~ s(x, k = 3))+
scale_color_viridis_d()+
xlab("Time in hours")+
ylab("Germinated conidia (%)")+
labs(colour = "Temperature")+
theme_classic() +
ggtitle("Germinated conidia over time, subset by relative humidity")+
scale_x_continuous(limits = c(12,48), breaks = c(12,24,36,48))+
scale_y_continuous(limits = c(0,100), n.breaks = 5)
Germination increases in a linear fashion at all the relative humidity treatments with exception to 99 % humidity, which shows a somewhat non-linear trend in the 22\(^\circ\)C and 25\(^\circ\)C temperature treatments.
Let’s look at the same data in terms of VPD split by time.
i_p_labs <-c("12 hours","24 hours","36 hours","48 hours")
names(i_p_labs) <- c(12,24,36,48)
g_vpd[, lty := fcase(Tm == 22, "dotted",
Tm == 25, "dashed",
Tm == 28, "twodash")]
Figure_4 <-
g_vpd |>
ggplot(aes(x = vpd, y = germ_n,colour = factor(Tm))) +
geom_point()+
theme_bw()+
scale_color_grey(start = 0.75,end = 0)+
geom_smooth(method = "gam",
formula = y ~ bs(x,k=0.2),
size = 0.7,se = FALSE)+
facet_wrap(facets = c("i_period"),
labeller = labeller(i_period = i_p_labs))+
xlab("Vapour pressure deficit (kPa)")+
ylab("Germinated conidia (%)")+
labs(colour = "Temperature")+
theme_bw()+
theme(strip.background = element_rect(fill = "white"))
Figure_4
Let’s look at the same data in terms of VPD split by temperature.
g_vpd |>
ggplot(aes(x = vpd, y = germ_n, colour = as.factor(i_period))) +
geom_point()+
facet_wrap(facets = c("Tm"),nrow = 1)+
geom_smooth(method = "gam",formula = y ~ bs(x,k=0.2), alpha = 0.3)+
scale_color_viridis_d()+
xlab("VPD")+
ylab("Germinated conidia (%)")+
labs(colour = "Incubation time")+
theme_bw()+
theme(strip.background = element_rect(fill = "white"))
This data is very hard to fit a curve to and seems like it would benefit with two functions. Lets plot again without the saturated 99% humidity treatment.
i_p_labs <-c("12 hours","24 hours","36 hours","48 hours")
names(i_p_labs) <- c(12,24,36,48)
g_vpd[g_vpd$RH != 99] |>
ggplot(aes(x = vpd, y = germ_n, colour = as.factor(Tm))) +
geom_point()+
facet_wrap(facets = c("i_period"),
labeller = labeller(i_period = i_p_labs))+
geom_smooth(method = "gam",formula = y ~ s(x,k=5), alpha = 0.1)+
scale_color_viridis_d()+
xlab("VPD")+
ylab("Germinated conidia (%)")+
labs(colour = "Temperature")+
theme_classic()
From these plots it seems that the influence of VPD on germination is bimodal, with an increase in germination between 1.5 and 2 kPa. Particularly in the warmer treatments.
VPD influence on germination - analysis
Lets start the analysis with a simple model with smoothing terms on
VPD, temperature and infection period. We are using a binomial general
additive model for all the germination analysis, because it is binary
count data. As we saw from plotting the raw data, we likely have
non-linear predictors and a gam model is a great way to
estimate with splines on the predictors.
m1v <- gam(cbind(germ_n,non_germ_n) ~ s(vpd, k = 5) + s(Tm, k = 3) + s(i_period, k = 3),
data = g_vpd,
family = "binomial")
summary(m1v)
Family: binomial
Link function: logit
Formula:
cbind(germ_n, non_germ_n) ~ s(vpd, k = 5) + s(Tm, k = 3) + s(i_period,
k = 3)
Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.97609 0.02141 -92.31 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Approximate significance of smooth terms:
edf Ref.df Chi.sq p-value
s(vpd) 3.995 4.000 4368.0 <2e-16 ***
s(Tm) 1.911 1.992 160.1 <2e-16 ***
s(i_period) 1.966 1.999 679.5 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R-sq.(adj) = 0.834 Deviance explained = 82%
UBRE = 3.5896 Scale est. = 1 n = 288plot(m1v)
Overall the model fits very well with a R\(^2\) of 0.8337. The splines on
Tm and i_period were statistically
significant, and seem to reflect the raw data. Warmer temperatures, from
the three temperatures included (22, 25 and 28\(^\circ\)C) lead to more conidia
germinating. Incubation time also showed a positive relationship with
the number of germinated conidia, i.e.. the longer the incubation the
more conidia were likely to have germinated. VPD also was a significant
predictor however the wiggliness might be able to be improved.
As earlier, lets try and reduce the willingness by using P-splines.
m2v <- gam(cbind(germ_n,non_germ_n) ~ s(vpd, k = 5, bs = "ps") + s(Tm, k = 3) + s(i_period, k = 3),
data = g_vpd,
family = "binomial")
summary(m2v)
Family: binomial
Link function: logit
Formula:
cbind(germ_n, non_germ_n) ~ s(vpd, k = 5, bs = "ps") + s(Tm,
k = 3) + s(i_period, k = 3)
Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.98330 0.02141 -92.65 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Approximate significance of smooth terms:
edf Ref.df Chi.sq p-value
s(vpd) 3.996 4.000 4652.2 <2e-16 ***
s(Tm) 1.903 1.991 170.4 <2e-16 ***
s(i_period) 1.968 1.999 693.6 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R-sq.(adj) = 0.889 Deviance explained = 86.5%
UBRE = 2.4425 Scale est. = 1 n = 288plot(m2v)
This shows only a small improvement and a decreasing trend in the VPD plot.
Next lets try more VPD spline knots.
m3v <- gam(cbind(germ_n,non_germ_n) ~ s(vpd, k = 7, bs = "ps") + s(Tm, k = 3) + s(i_period, k = 3),
data = g_vpd,
family = "binomial")
summary(m3v)
Family: binomial
Link function: logit
Formula:
cbind(germ_n, non_germ_n) ~ s(vpd, k = 7, bs = "ps") + s(Tm,
k = 3) + s(i_period, k = 3)
Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.96325 0.02119 -92.63 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Approximate significance of smooth terms:
edf Ref.df Chi.sq p-value
s(vpd) 5.971 5.999 5167.6 <2e-16 ***
s(Tm) 1.726 1.924 154.1 <2e-16 ***
s(i_period) 1.971 1.999 717.8 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R-sq.(adj) = 0.95 Deviance explained = 93%
UBRE = 0.82596 Scale est. = 1 n = 288plot(m3v)
While there exists a lot of willingness in the gam, the VPD plot shows the increase in germination around 1.5 kPa and a rapid decrease in germination between the saturated treatments (0.026 - 0.038 kPa), and the next highest humidity (0.158 - 0.189 kPa). This model also provides a high R-squared (0.9502) and low UBRE 0.826, indicating good predictive ability while being less affected by drop one, tests.
Temperature does not seem like it needs a non-linear term. The next model will test this assumption.
m4v <- gam(cbind(germ_n,non_germ_n) ~ s(vpd, k = 7) + Tm + s(i_period,k=3),
data = g_vpd,
family = "binomial")
summary(m4v)
Family: binomial
Link function: logit
Formula:
cbind(germ_n, non_germ_n) ~ s(vpd, k = 7) + Tm + s(i_period,
k = 3)
Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.009685 0.190921 -21.00 <2e-16 ***
Tm 0.080869 0.007478 10.81 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Approximate significance of smooth terms:
edf Ref.df Chi.sq p-value
s(vpd) 5.990 6.000 4868.3 <2e-16 ***
s(i_period) 1.969 1.999 705.5 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R-sq.(adj) = 0.928 Deviance explained = 90.1%
UBRE = 1.5636 Scale est. = 1 n = 288plot(m4v)
There is a slight reduction in the model fit, as shown by the difference in the R\(^2\) and UBRE values between m3v and m4v.
In addition the wiggliness of the VPD curve seems more pronounced. We will next evaluate if there is an interaction between VPD and temperature, which should also help us decide if temperature should be a linear or non-linear term.
m5v <- gam(cbind(germ_n,non_germ_n) ~ te(vpd,Tm, k = 3)+ s(i_period,k=3),
data = g_vpd,
family = "binomial")
summary(m5v)
Family: binomial
Link function: logit
Formula:
cbind(germ_n, non_germ_n) ~ te(vpd, Tm, k = 3) + s(i_period,
k = 3)
Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.95016 0.02132 -91.47 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Approximate significance of smooth terms:
edf Ref.df Chi.sq p-value
te(vpd,Tm) 7.867 7.987 3652.8 <2e-16 ***
s(i_period) 1.965 1.999 650.1 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R-sq.(adj) = 0.697 Deviance explained = 71%
UBRE = 6.3679 Scale est. = 1 n = 288plot(m5v)
anova(m4v,m5v)
Analysis of Deviance Table
Model 1: cbind(germ_n, non_germ_n) ~ s(vpd, k = 7) + Tm + s(i_period,
k = 3)
Model 2: cbind(germ_n, non_germ_n) ~ te(vpd, Tm, k = 3) + s(i_period,
k = 3)
Resid. Df Resid. Dev Df Deviance
1 278.00 718.41
2 277.01 2100.29 0.98715 -1381.9This comparison shows while the interaction between temperature is significant, the m5v model does not fit the data as well as m4v .
Because of the rapid drop in germination in the lower VPD range, and little change between the higher VPD treatments, it is difficult to fit a non-linear model with evenly spaced splines. Also we don’t have enough evenly spaced VPD treatments to increase the number of knots without causing an over fit model. We can override this default behaviour of splines placed at set points by specifying our own knot points.
Here we will specify two knots at 0.1 kPa and 0.4 kPa, in addition to the boundary knots.
b_k <- range(g_vpd$vpd) # set spline boundary knots
m6v <- gam(cbind(germ_n,non_germ_n) ~ s(vpd, k = 4,bs = "cs") + Tm + s(i_period,k=3),
knots = list(vpd = c(b_k[1],0.1,0.4,b_k[2])),
data = g_vpd,
family = "binomial")
summary(m6v)
Family: binomial
Link function: logit
Formula:
cbind(germ_n, non_germ_n) ~ s(vpd, k = 4, bs = "cs") + Tm + s(i_period,
k = 3)
Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.48276 0.19318 -23.20 <2e-16 ***
Tm 0.10138 0.00757 13.39 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Approximate significance of smooth terms:
edf Ref.df Chi.sq p-value
s(vpd) 2.999 3.000 5302.4 <2e-16 ***
s(i_period) 1.971 1.999 723.7 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R-sq.(adj) = 0.958 Deviance explained = 94.1%
UBRE = 0.52704 Scale est. = 1 n = 288plot(m6v)
This model presents the best fit model so far, as shown by the R\(^2\) value of 0.9582, 94.12% deviance explained and 0.5270411UBRE.
Let’s compare the models side-by-side using the ‘Akaike Information Criterion’ (AIC). A lower number is generally viewed as indicating a better model fit.
AIC(m1v,m2v,m3v, m4v, m5v, m6v)
df AIC
m1v 8.872701 2477.644
m2v 8.867062 2147.279
m3v 10.667524 1681.719
m4v 9.958923 1894.166
m5v 10.831486 3277.794
m6v 6.970012 1595.630Model six (m6v) shows the best fit to the raw data
plotted in Figure 4.
# subset by incubation time
par(mfrow = c(2,2),
oma = c(0,0,0,0),
mar = c(0.1,0.1,1,0.1))
vis.gam(m6v,
phi = 20,
theta = 45,
view = c("vpd","Tm"),
cond = list(i_period = 12),
main = "12 hours")
vis.gam(m6v,
phi = 20,
theta = 45,
view = c("vpd","Tm"),
cond = list(i_period = 24),
main = "24 hours")
vis.gam(m6v,
phi = 20,
theta = 45,
view = c("vpd","Tm"),
cond = list(i_period = 36),
main = "36 hours")
vis.gam(m6v,
phi = 20,
theta = 45,
view = c("vpd","Tm"),
cond = list(i_period = 48),
main = "48 hours")
Visualise model against the raw data
To simulate data from the model using predict we need to
provide some values which to estimate over.
# Make a data.frame of test values to feed in the model
test_vpd <- rep(seq(from = 0, to = max(g_vpd$vpd), by = 0.01), 3*4)
test_dat <- data.table(vpd = test_vpd, # between 5 and 8 degrees for each incubation period
Tm = rep(unique(g_vpd$Tm),
each = length(test_vpd),
times = 4),
i_period = rep(c(12,24,36,48), each = length(test_vpd)*3))
# create a data.frame with predicted values and reference them to the respective temperature
m6_pred_list <-
predict(m6v,
test_dat,
type = "response",
se.fit = TRUE)
m6_pred <- data.table(vpd = test_dat$vpd,
Tm = as.factor(test_dat$Tm),
i_period = as.factor(test_dat$i_period),
germinated = m6_pred_list$fit,
germ_se = m6_pred_list$se.fit,
CI_low = m6_pred_list$fit + (m6_pred_list$se.fit * -1.94),
CI_high = m6_pred_list$fit + (m6_pred_list$se.fit * 1.94))
Figure_4 <-
m6_pred |>
ggplot(aes(x = vpd, y = germinated*100, colour = i_period))+
facet_wrap(facets = c("Tm"))+
#labeller = labeller(i_period = i_p_labs))+
geom_line(size = 0.8)+
scale_color_grey(start = 0.75,end = 0)+
theme_bw()+
geom_line(aes(x = vpd, y = CI_low*100), linetype = "dashed", size = 0.6)+
geom_line(aes(x = vpd, y = CI_high*100), linetype = "dashed", size = 0.6)+
xlab("Vapour pressure deficit (kPa)")+
ylab("Germinated conidia (%)")+
labs(colour = "Incubation\nperiod (hours)")+
geom_point(data = g_vpd,
aes(x = vpd, y = (germ_n),
colour = as.factor(i_period)))+
theme(strip.background = element_rect(fill = "white"))
Figure_4
Impact of VPD on P. xanthii branching
Import the clean data
There were three data-sets which recorded germ-tube formation at a range of VPD.
Lets look at the summaries of the data.
rmarkdown::paged_table(b_vpd[, list(
Mean = mean(conidia),
Std_dev = round(sd(conidia), 3),
Median = median(conidia),
Range = paste0(range(conidia), collapse = " - "),
Replicates = .N),
by = .(vpd, Tm, germtubes)])
rmarkdown::paged_table(b_vpd2[, list(
Mean = mean(conidia),
Std_dev = round(sd(conidia), 3),
Median = median(conidia),
Range = paste0(range(conidia), collapse = " - "),
Replicates = .N),
by = .(vpd, Tm, germtubes)])
rmarkdown::paged_table(b_vpd3[, list(
Mean = mean(conidia),
Std_dev = round(sd(conidia), 3),
Median = median(conidia),
Range = paste0(range(conidia), collapse = " - "),
Replicates = .N),
by = .(vpd, Tm, germtubes)])
We will combine the data into one object to interrogate.
b_v <-rbind(b_vpd,b_vpd2,b_vpd3)
Visualise the raw data
First we want observe the effect of VPD on the number of germtubes at 36 hours for each experimental rep.
b_v[,percent := (conidia/sum(conidia))*100,
by = .(vpd,exp)]
b_v |>
ggplot(aes( x = as.factor(vpd),
y = percent ,
fill = as.factor(germtubes))) +
geom_col()+
facet_grid(rows = vars(exp))+
scale_fill_brewer(palette = "Blues", name = "Germ tubes")+
xlab("Vapour pressure deficit")+
ylab("Percent germinated")+
labs(colour = "Germ tubes")+
theme_classic()+
theme(axis.text.x = element_text(angle = 70, vjust = 1.7, hjust = 2))+
geom_text(aes(x = as.factor(vpd), y = 115, label = paste0(Tm,"°C")),size = 2.5)+
scale_y_continuous(limits = c(0,120), breaks = c(0,25,50,75,100))
Each experiment has similar trends and different variation for each VPD
probably due to the temperature incubation.
Next lets look at a consolidated graph with all the experimental reps together.
b_v[,percent := (conidia/sum(conidia))*100,
by = .(vpd)]
Figure_5 <-
b_v |>
ggplot(aes( x = as.factor(vpd),
y = percent ,
fill = as.factor(germtubes))) +
geom_col()+
scale_fill_manual(values = c("0" = "grey90",
"1" = "grey70",
"2" = "grey55",
"3" = "grey40"),
name = "Germ tubes")+
xlab("Vapour pressure deficit (kPa)")+
ylab("Percent germinated (%)")+
labs(colour = "Germ tubes")+
theme_classic()+
scale_fill_grey(start = 0.95,end = 0.25)+
theme(axis.text.x = element_text(angle = 70, vjust = 1.7, hjust = 2))+
geom_text(aes(x = as.factor(vpd), y = 105, label = paste0(Tm,"°C")),size = 2.5)+
scale_y_continuous(limits = c(0,111), breaks = c(0,25,50,75,100))+
annotate("text", x = 1, y = 111, label = "(a)", size =5)
Figure_5
Analysis of conidial germination and germtube proliferation
In order to analyse the data we need to set the data classes.
b_v[, c("exp",
"rep",
"germtubes",
"cham",
"vpd") := list(as.factor(exp),
as.integer(rep),
as.ordered(germtubes),
paste0(Tm, RH),
as.numeric(vpd))]
First we start with a simple additive model with no interaction or smooths.
Call:
polr(formula = germtubes ~ vpd + Tm + exp, data = b_v, weights = conidia,
Hess = TRUE)
Coefficients:
Value Std. Error t value
vpd -1.81748 0.039186 -46.381
Tm 0.05095 0.007146 7.129
exp2 -0.78489 0.044774 -17.530
exp3 -0.59512 0.040611 -14.654
Intercepts:
Value Std. Error t value
0|1 0.9503 0.1794 5.2981
1|2 1.5189 0.1795 8.4605
2|3 2.0327 0.1801 11.2889
Residual Deviance: 28328.77
AIC: 28342.77 Lets try predicting from this model so we can plot how it performs. First we need to make some test data.
# create a data.frame with predicted values and reference them to the respective temperature
m1_pred <-
as.data.table(cbind(test_dat,
round(predict(m1,
test_dat,
type = "p"),
3)))
# format to long
m1_pred <-
melt(
m1_pred,
id.vars = c("Tm","vpd","exp"),
measure.vars = c("0", "1", "2", "3"),
variable.name = "germtubes",
value.name = "prob"
)
# Plot the predicted model
m1_pred |>
ggplot(aes(x = vpd, y = prob, colour = germtubes))+
facet_grid(rows = vars(Tm), cols = vars(exp))+
geom_line(size = 1)+
theme_minimal()
There is a lot going on in this graph but as the m1
summary shows experimental rep, temperature and VPD were all
significant. While there might be a difference in the predicted values
between experiments, the trends are the same. Nothing drastically
different is happening between experimental rep which justifies keeping
them separate. For ease in interpretation, lets look at a consolidated
plot, taking the mean of each experiment.
m1_pred[,.(prob = mean(prob)),by= .(vpd,germtubes)] |>
ggplot(aes(x = vpd, y = prob, colour = germtubes))+
# facet_wrap(facets = vars(Tm))+
geom_line(size = 1)+
theme_minimal()
Model one tells us a bit about the influence of the VPD and temperature
variables but the model in the graph above does look like it will fit
our data very well. Model two we will try fit the model to our
observations by using a basis spline on our VPD variable.
m2 <- polr(germtubes ~ bs(vpd,knots = 0.5) + Tm + exp,
weights = conidia,
data = b_v, Hess = TRUE)
summary(m2)
Call:
polr(formula = germtubes ~ bs(vpd, knots = 0.5) + Tm + exp, data = b_v,
weights = conidia, Hess = TRUE)
Coefficients:
Value Std. Error t value
bs(vpd, knots = 0.5)1 -4.80829 0.095654 -50.27
bs(vpd, knots = 0.5)2 -2.23386 0.160078 -13.95
bs(vpd, knots = 0.5)3 -5.39222 0.194772 -27.68
bs(vpd, knots = 0.5)4 -4.44845 0.113469 -39.20
Tm 0.07002 0.007992 8.76
exp2 -0.69700 0.048324 -14.42
exp3 -0.77686 0.045167 -17.20
Intercepts:
Value Std. Error t value
0|1 -0.4810 0.2013 -2.3896
1|2 0.3077 0.2011 1.5304
2|3 1.0436 0.2014 5.1818
Residual Deviance: 24437.32
AIC: 24457.32 Let’s predict from this model to consider it’s performance.
# create a data.frame with predicted values and reference them to the respective temperature
m2_pred <- marginaleffects::predictions(m2,
newdata = test_dat,
type = "probs")
# set data.table format
setDT(m2_pred)
# calculate approximate confidence intervals
m2_pred[, CI_low := .(predicted - 1.96*std.error)]
m2_pred[, CI_high := .(predicted + 1.96*std.error)]
# change the column names
setnames(m2_pred,
old = c("group","Tm.x","vpd.x"),
new = c("germtubes", "Tm", "vpd"))
# summarise the 3 experiments into one
m2_plot <- m2_pred[,.(prob = mean(predicted),
CI_l = mean(CI_low),
CI_h = mean(CI_high)),
by= .(vpd,Tm, germtubes)]
m2_plot |>
ggplot(aes(x = vpd, y = prob, colour = germtubes))+
geom_line(size = 0.7)+
facet_wrap(facets = vars(Tm), ncol = 2)+
geom_line(data = m2_plot, aes(x = vpd, y = CI_l, #colour = group,
), linetype = "dotted", size = 0.6)+
geom_line(data = m2_plot, aes(x = vpd, y = CI_h, #colour = germtubes,
), linetype = "dotted", size = 0.6)+
theme_minimal()
There is a consistent effect of temperature that increases the number of germinated and branching germ-tubes. Low VPD has a strong effect on the number and proportion of branching germtubes. While the number of germinated conidia and germ tube branching decreases rapidly as VPD increases, however after 0.5 kPa, the proportion of germinating and branching conidia stabilise all the way up to 2 kPa.
Lets look at just VPD and drop the temperature effect.
# re-summarise the three experiments into one without looking at temperatures
m2_plot <- m2_pred[,.(prob = mean(predicted),
CI_l = mean(CI_low),
CI_h = mean(CI_high)),
by= .(vpd, germtubes)]
# plot with confidence intervals
Figure_5b <-
m2_plot |>
ggplot(aes(x = vpd, y = prob, colour = germtubes))+
geom_line(size = 0.7)+
geom_line(data = m2_plot, aes(x = vpd, y = CI_l, colour = germtubes,
), linetype = "dotted", size = 0.6)+
geom_line(data = m2_plot, aes(x = vpd, y = CI_h, colour = germtubes,
), linetype = "dotted", size = 0.6)+
theme_minimal() +
annotate("text", x = 0.01, y = 0.98, label = "(b)", size =5)+
scale_color_grey(start = 0.8, end = 0.3)+
xlab("Vapour pressure deficit (kPa)")+
ylab("Proportional probability")+
labs(colour = "Germ tubes")
Figure_5b
# Save table as a word document
ggsave(
filename = "FiguresTables/Fig5b_VPDGermtubeMod.tiff",
plot = Figure_5b,
device = "tiff",
height = 5,
width = 7,
dpi = 300,
units = "in"
)
Given the plots show a decent similarity to the raw data we will
choose m2. Lets now extract the coefficients
# get table of coefficients
coef_table <- data.table(coef(summary(m2)), keep.rownames = TRUE)
# calculate p-values
coef_table[, p_value := pnorm(abs(`t value`),lower.tail = FALSE) *2][, sig := p_star(p_value)]
coef_table
rn Value Std. Error t value
1: bs(vpd, knots = 0.5)1 -4.80829167 0.095653746 -50.267677
2: bs(vpd, knots = 0.5)2 -2.23386481 0.160077935 -13.954858
3: bs(vpd, knots = 0.5)3 -5.39221887 0.194771549 -27.684838
4: bs(vpd, knots = 0.5)4 -4.44845006 0.113469356 -39.203977
5: Tm 0.07001695 0.007992462 8.760374
6: exp2 -0.69700190 0.048324305 -14.423423
7: exp3 -0.77686230 0.045166989 -17.199781
8: 0|1 -0.48104065 0.201309824 -2.389554
9: 1|2 0.30770150 0.201058233 1.530410
10: 2|3 1.04364994 0.201407804 5.181775
p_value sig
1: 0.000000e+00 ***
2: 2.938875e-44 ***
3: 1.063089e-168 ***
4: 0.000000e+00 ***
5: 1.946036e-18 ***
6: 3.685972e-47 ***
7: 2.665388e-66 ***
8: 1.686885e-02 *
9: 1.259153e-01
10: 2.197842e-07 ***