8 1992 full system analyses
8.1 Total leaf area in 1992.
This analysis uses all systems that have a non-NA value for total leaf area of the front and back in 1991, excluding those systems for which 1990 leaf area was NA.
8.1.1 Data & visualization
Read in data
# set na.strings="." to recode '.' to NA
mayappleDatasetFull<-readxl::read_excel("Datasets/mayapple_1990_updated_version2.xlsx", na=".")Select the total leaf area variable, and the independent variables. Also filtering out rows that don’t have a leaf area record in 1990. We need to do this in order to be able to compare models from model set A (which don’t include leaf area, and so wouldn’t naturally be affected by the missing data) with models from model set B (which do include leaf area and would be affected by the missing data). In other words, removing these rows means that model set A and model set B will run using the exact same dataset, which is necessary in order for us to compare the models.
# select variables for analysis and filter out rows that don't have a leaf area recorded in 1990
mayappleData <- mayappleDatasetFull %>%
dplyr::select(COLONY,Lf_90,Sex_90,Sever,Time, Lftot_91, Lftot_92) %>%
mutate(Lftot_92=parse_double(Lftot_92), Lftot_91=parse_double(Lftot_91)) %>%
filter(!is.na(Lf_90)) %>%
filter(!is.na(Lftot_91)) %>%
filter(!is.na(Lftot_92))
ggplot(data=mayappleData, aes(x=Sever, y=Lftot_92, fill = Time)) +
stat_summary(fun = mean, geom = "bar",position="dodge") +
stat_summary(fun.data = mean_se, geom = "errorbar", position="dodge")+
facet_wrap(~Sex_90) +
ylab("Mean total leaf area \n in front and back shoots in 1992") +
scale_fill_brewer(palette="Set1")
8.1.2 Model selection with leaf area in 1990 and 1991
d_mod92la.0l<-lmer(Lftot_92~Sex_90*Sever*Time +Lftot_91+Lf_90+(1|COLONY), data = mayappleData, na.action = "na.fail", REML=FALSE)
#this is the global model
#MuMIn::getAllTerms(d_mod92la.0l) #this lets us see how to specify terms in the model
dr<-dredge(d_mod92la.0l, fixed=~Lf_90+Lftot_91)
knitr::kable(head(dr, n=10)) %>%
kable_styling() %>%
scroll_box(width = "100%", box_css = "border: 0px;")| (Intercept) | Lf_90 | Lftot_91 | Sever | Sex_90 | Time | Sever:Sex_90 | Sever:Time | Sex_90:Time | Sever:Sex_90:Time | df | logLik | AICc | delta | weight | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 6 | -18.98414 | 0.2778212 | 0.7753865 |
|
NA |
|
NA | NA | NA | NA | 10 | -1455.489 | 2932.022 | 0.0000000 | 0.18486575 |
| 2 | -11.40856 | 0.3133500 | 0.7830574 |
|
NA | NA | NA | NA | NA | NA | 8 | -1457.689 | 2932.055 | 0.0332159 | 0.18182086 |
| 1 | -14.75243 | 0.3049134 | 0.7980890 | NA | NA | NA | NA | NA | NA | NA | 5 | -1461.208 | 2932.694 | 0.6728440 | 0.13205360 |
| 5 | -24.40721 | 0.2699304 | 0.7914756 | NA | NA |
|
NA | NA | NA | NA | 7 | -1459.182 | 2932.888 | 0.8663366 | 0.11987644 |
| 4 | -44.03431 | 0.3534655 | 0.7845174 |
|
|
NA | NA | NA | NA | NA | 9 | -1457.382 | 2933.613 | 1.5908845 | 0.08344499 |
| 40 | -44.46533 | 0.2972305 | 0.7771511 |
|
|
|
NA | NA |
|
NA | 13 | -1453.031 | 2933.813 | 1.7913502 | 0.07548657 |
| 8 | -29.18421 | 0.2918686 | 0.7761747 |
|
|
|
NA | NA | NA | NA | 11 | -1455.459 | 2934.176 | 2.1543934 | 0.06295582 |
| 39 | -45.87488 | 0.2827734 | 0.7913642 | NA |
|
|
NA | NA |
|
NA | 10 | -1456.590 | 2934.223 | 2.2009571 | 0.06150702 |
| 3 | -40.25139 | 0.3382077 | 0.7985259 | NA |
|
NA | NA | NA | NA | NA | 6 | -1460.998 | 2934.387 | 2.3654655 | 0.05665029 |
| 7 | -30.35397 | 0.2786637 | 0.7917473 | NA |
|
|
NA | NA | NA | NA | 8 | -1459.171 | 2935.017 | 2.9956628 | 0.04133867 |
From this round of model selection, when Lf_90 and LfFtot_91 are required to be in the models that are considered, there are six models within two AICc units: (1) Lf_90, Lft_91, Sever, and Time (AICc 2932.0) (2) Lf_90, Lft_91, and Sever (AICc 2932.1) (3) Lf_90 and Lft_91 (AICc 2932.7) (4) Lf_90, Lft_91, and Time (AICc 2932.9) (5) Lf_90, Lft_91, Sever, and Sex_90 (AICc 2933.6) (6) Lf_90, Lft_91, Sever, Time, Sex_90, and the Sex_90 x Time interaction (AICc 2933.8)
8.1.3 Model selection with leaf area in 1991
dr<-dredge(d_mod92la.0l, fixed=~Lftot_91)
knitr::kable(head(dr, n=10)) %>%
kable_styling() %>%
scroll_box(width = "100%", box_css = "border: 0px;")| (Intercept) | Lf_90 | Lftot_91 | Sever | Sex_90 | Time | Sever:Sex_90 | Sever:Time | Sex_90:Time | Sever:Sex_90:Time | df | logLik | AICc | delta | weight | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 12 | -18.98414 | 0.2778212 | 0.7753865 |
|
NA |
|
NA | NA | NA | NA | 10 | -1455.489 | 2932.022 | 0.0000000 | 0.18486575 |
| 4 | -11.40856 | 0.3133500 | 0.7830574 |
|
NA | NA | NA | NA | NA | NA | 8 | -1457.689 | 2932.055 | 0.0332159 | 0.18182086 |
| 2 | -14.75243 | 0.3049134 | 0.7980890 | NA | NA | NA | NA | NA | NA | NA | 5 | -1461.208 | 2932.694 | 0.6728440 | 0.13205360 |
| 10 | -24.40721 | 0.2699304 | 0.7914756 | NA | NA |
|
NA | NA | NA | NA | 7 | -1459.182 | 2932.888 | 0.8663366 | 0.11987644 |
| 8 | -44.03431 | 0.3534655 | 0.7845174 |
|
|
NA | NA | NA | NA | NA | 9 | -1457.382 | 2933.613 | 1.5908845 | 0.08344499 |
| 80 | -44.46533 | 0.2972305 | 0.7771511 |
|
|
|
NA | NA |
|
NA | 13 | -1453.031 | 2933.813 | 1.7913502 | 0.07548657 |
| 16 | -29.18421 | 0.2918686 | 0.7761747 |
|
|
|
NA | NA | NA | NA | 11 | -1455.459 | 2934.176 | 2.1543934 | 0.06295582 |
| 78 | -45.87488 | 0.2827734 | 0.7913642 | NA |
|
|
NA | NA |
|
NA | 10 | -1456.590 | 2934.223 | 2.2009571 | 0.06150702 |
| 6 | -40.25139 | 0.3382077 | 0.7985259 | NA |
|
NA | NA | NA | NA | NA | 6 | -1460.998 | 2934.387 | 2.3654655 | 0.05665029 |
| 14 | -30.35397 | 0.2786637 | 0.7917473 | NA |
|
|
NA | NA | NA | NA | 8 | -1459.171 | 2935.017 | 2.9956628 | 0.04133867 |
From this round of model selection, when LfFtot_91 only is required to be in the models that are considered, the same six best models are returned as above.
8.1.4 Model selection with leaf area in 1990
dr<-dredge(d_mod92la.0l, fixed=~Lf_90)
knitr::kable(head(dr, n=10)) %>%
kable_styling() %>%
scroll_box(width = "100%", box_css = "border: 0px;")| (Intercept) | Lf_90 | Lftot_91 | Sever | Sex_90 | Time | Sever:Sex_90 | Sever:Time | Sex_90:Time | Sever:Sex_90:Time | df | logLik | AICc | delta | weight | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 12 | -18.98414 | 0.2778212 | 0.7753865 |
|
NA |
|
NA | NA | NA | NA | 10 | -1455.489 | 2932.022 | 0.0000000 | 0.18486575 |
| 4 | -11.40856 | 0.3133500 | 0.7830574 |
|
NA | NA | NA | NA | NA | NA | 8 | -1457.689 | 2932.055 | 0.0332159 | 0.18182086 |
| 2 | -14.75243 | 0.3049134 | 0.7980890 | NA | NA | NA | NA | NA | NA | NA | 5 | -1461.208 | 2932.694 | 0.6728440 | 0.13205360 |
| 10 | -24.40721 | 0.2699304 | 0.7914756 | NA | NA |
|
NA | NA | NA | NA | 7 | -1459.182 | 2932.888 | 0.8663366 | 0.11987644 |
| 8 | -44.03431 | 0.3534655 | 0.7845174 |
|
|
NA | NA | NA | NA | NA | 9 | -1457.382 | 2933.613 | 1.5908845 | 0.08344499 |
| 80 | -44.46533 | 0.2972305 | 0.7771511 |
|
|
|
NA | NA |
|
NA | 13 | -1453.031 | 2933.813 | 1.7913502 | 0.07548657 |
| 16 | -29.18421 | 0.2918686 | 0.7761747 |
|
|
|
NA | NA | NA | NA | 11 | -1455.459 | 2934.176 | 2.1543934 | 0.06295582 |
| 78 | -45.87488 | 0.2827734 | 0.7913642 | NA |
|
|
NA | NA |
|
NA | 10 | -1456.590 | 2934.223 | 2.2009571 | 0.06150702 |
| 6 | -40.25139 | 0.3382077 | 0.7985259 | NA |
|
NA | NA | NA | NA | NA | 6 | -1460.998 | 2934.387 | 2.3654655 | 0.05665029 |
| 14 | -30.35397 | 0.2786637 | 0.7917473 | NA |
|
|
NA | NA | NA | NA | 8 | -1459.171 | 2935.017 | 2.9956628 | 0.04133867 |
From this round of model selection, when Lf_90 is required to be in the models that are considered, the same six best models are returned as above.
8.1.5 Model selection without leaf area
## Fixed term is "(Intercept)"knitr::kable(head(dr, n=10)) %>%
kable_styling() %>%
scroll_box(width = "100%", box_css = "border: 0px;")| (Intercept) | Lf_90 | Lftot_91 | Sever | Sex_90 | Time | Sever:Sex_90 | Sever:Time | Sex_90:Time | Sever:Sex_90:Time | df | logLik | AICc | delta | weight | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 24 | -18.98414 | 0.2778212 | 0.7753865 |
|
NA |
|
NA | NA | NA | NA | 10 | -1455.489 | 2932.022 | 0.0000000 | 0.18486575 |
| 8 | -11.40856 | 0.3133500 | 0.7830574 |
|
NA | NA | NA | NA | NA | NA | 8 | -1457.689 | 2932.055 | 0.0332159 | 0.18182086 |
| 4 | -14.75243 | 0.3049134 | 0.7980890 | NA | NA | NA | NA | NA | NA | NA | 5 | -1461.208 | 2932.694 | 0.6728440 | 0.13205360 |
| 20 | -24.40721 | 0.2699304 | 0.7914756 | NA | NA |
|
NA | NA | NA | NA | 7 | -1459.182 | 2932.888 | 0.8663366 | 0.11987644 |
| 16 | -44.03431 | 0.3534655 | 0.7845174 |
|
|
NA | NA | NA | NA | NA | 9 | -1457.382 | 2933.613 | 1.5908845 | 0.08344499 |
| 160 | -44.46533 | 0.2972305 | 0.7771511 |
|
|
|
NA | NA |
|
NA | 13 | -1453.031 | 2933.813 | 1.7913502 | 0.07548657 |
| 32 | -29.18421 | 0.2918686 | 0.7761747 |
|
|
|
NA | NA | NA | NA | 11 | -1455.459 | 2934.176 | 2.1543934 | 0.06295582 |
| 156 | -45.87488 | 0.2827734 | 0.7913642 | NA |
|
|
NA | NA |
|
NA | 10 | -1456.590 | 2934.223 | 2.2009571 | 0.06150702 |
| 12 | -40.25139 | 0.3382077 | 0.7985259 | NA |
|
NA | NA | NA | NA | NA | 6 | -1460.998 | 2934.387 | 2.3654655 | 0.05665029 |
| 28 | -30.35397 | 0.2786637 | 0.7917473 | NA |
|
|
NA | NA | NA | NA | 8 | -1459.171 | 2935.017 | 2.9956628 | 0.04133867 |
From this round of model selection, when the leaf area effects are not required to be in the models that are considered, the same six best models are returned as above.
8.1.6 Conclusion
The simplest and best model here includes fixed effects of total leaf area in 1991 and leaf area in 1990.
8.1.7 Estimated marginal means
the main effect of leaf area in 1990
#generated based on the min and max of Lf_90
mylist <- list(Lf_90=seq(150,1250,by=100))
g1 <- ggplot(data=mayappleData, aes(x=Lf_90, y=Lftot_92, color=Sever)) +
geom_point(alpha=0.5) +
ylab("Mean total leaf area \n in front and back shoots in 1992")
g2 <- emmip(d_mod92la.bm, ~ Lf_90 , at=mylist, CI=TRUE)
gridExtra::grid.arrange(g1,g2,nrow=1)
While there are a few systems that had leaf area in 1990 but zero leaf area in 1992, there is an overall positive relationship.
the main effect of leaf area in 1991
mylist <- list(Lftot_91=seq(0,1850,by=100))
g1 <- ggplot(data=mayappleData, aes(x=Lftot_91, y=Lftot_92, color=Sever)) +
geom_point(alpha=0.5) +
ylab("Mean total leaf area \n in front and back shoots in 1992")
g2 <- emmip(d_mod92la.bm, ~ Lftot_91 , at=mylist, CI=TRUE)
gridExtra::grid.arrange(g1,g2,nrow=1)
Aside from a small number of systems that had leaf area in 1991 but zero leaf area in 1992, there is an overall positive relationship.