library(tidyverse)
url <- "https://raw.githubusercontent.com/jlacasa/stat705_fall2024/main/classes/data/lotus_part2.csv"
dd <- read.csv(url)
dd %>%
filter(species %in% c("A")) %>%
ggplot(aes(doy, crown_g))+
geom_point(shape = 21, size =3.5, fill = "black")+
labs(x = "Day of the Year",
y = "Crown biomass (g)",
fill = "Treatment")+
theme_classic()+
theme(aspect.ratio = 1)
\[y_i = \beta_0 + \beta_1x_i + \varepsilon_i, \\ \varepsilon_i\sim N(0,\sigma^2) .\]
m1 <- lm(crown_g ~ doy, data = dd)
summary(m1)
Call:
lm(formula = crown_g ~ doy, data = dd)
Residuals:
Min 1Q Median 3Q Max
-0.4178 -0.1364 -0.0201 0.1081 0.7966
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.3722453 0.1363156 -10.07 5.51e-16 ***
doy 0.0091941 0.0006899 13.33 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.1961 on 82 degrees of freedom
Multiple R-squared: 0.6841, Adjusted R-squared: 0.6803
F-statistic: 177.6 on 1 and 82 DF, p-value: < 2.2e-16dd %>%
filter(species %in% c("A")) %>%
ggplot(aes(doy, crown_g))+
stat_function(fun = function(x){
coef(m1)[1] + coef(m1)[2]*x
})+
geom_point(shape = 21, size =3.5, fill = "black")+
labs(x = "Day of the Year",
y = "Crown biomass (g)",
fill = "Treatment")+
theme_classic()+
theme(aspect.ratio = 1)
\[y_{ij} = \beta_{0j} + \beta_1x_{ij} + \varepsilon_{ij}, \\ \varepsilon_{ij}\sim N(0,\sigma^2) .\]
m2 <- lm(crown_g ~ trt + doy, data = dd)
summary(m2)
Call:
lm(formula = crown_g ~ trt + doy, data = dd)
Residuals:
Min 1Q Median 3Q Max
-0.40020 -0.13449 -0.00894 0.09488 0.78319
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.359492 0.137887 -9.859 1.6e-15 ***
trtflood -0.030954 0.043359 -0.714 0.477
doy 0.009197 0.000692 13.291 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.1967 on 81 degrees of freedom
Multiple R-squared: 0.6861, Adjusted R-squared: 0.6784
F-statistic: 88.53 on 2 and 81 DF, p-value: < 2.2e-16dd %>%
ggplot(aes(doy, crown_g))+
stat_function(fun = function(x){
coef(m2)[1] + coef(m2)[3]*x
})+
stat_function(fun = function(x){
(coef(m2)[1] + coef(m2)[2]) + coef(m2)[3]*x
},
linetype =2)+
geom_point(shape = 21, size =3.5, aes(fill = trt))+
scale_fill_manual(values = c("#DB504A", "#43AA8B"))+
labs(x = "Day of the Year",
y = "Crown biomass (g)",
fill = "Treatment")+
theme_classic()+
theme(aspect.ratio = 1)
\[y_{ij} = \beta_{0j} + \beta_{1j} x_{ij} + \varepsilon_{ij}, \\ \varepsilon_{ij}\sim N(0,\sigma^2) .\]
m3 <- lm(crown_g ~ trt * doy, data = dd)
summary(m3)
Call:
lm(formula = crown_g ~ trt * doy, data = dd)
Residuals:
Min 1Q Median 3Q Max
-0.35287 -0.13701 -0.01415 0.09433 0.73550
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.5809069 0.1908532 -8.283 2.25e-12 ***
trtflood 0.4123418 0.2706453 1.524 0.132
doy 0.0103322 0.0009681 10.673 < 2e-16 ***
trtflood:doy -0.0022714 0.0013692 -1.659 0.101
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.1946 on 80 degrees of freedom
Multiple R-squared: 0.6966, Adjusted R-squared: 0.6852
F-statistic: 61.21 on 3 and 80 DF, p-value: < 2.2e-16plot_data <- expand.grid(trt=c("ctrl", "flood"),
doy = min(dd$doy):max(dd$doy))
plot_data <- plot_data %>%
bind_cols(as.data.frame(predict(m2, newdata = plot_data, interval = "confidence")))dd %>%
filter(species %in% c("A")) %>%
ggplot(aes(doy, crown_g))+
stat_function(fun = function(x){
coef(m3)[1] + coef(m3)[3]*x
})+
stat_function(fun = function(x){
(coef(m3)[1] + coef(m3)[2]) + (coef(m3)[3] + coef(m3)[4])*x
},
linetype =2)+
geom_ribbon(aes(ymin = lwr, ymax = upr, x = doy, y = fit, fill = trt), data = plot_data, alpha = .2)+
geom_point(shape = 21, size =3.5, aes(fill = trt))+
scale_fill_manual(values = c("#DB504A", "#43AA8B"))+
labs(x = "Day of the Year",
y = "Crown biomass (g)",
fill = "Treatment")+
theme_classic()+
theme(aspect.ratio = 1)