I am working on a multipanel plot using the package ggplot2 in R. My dataset consists of 20 species, 6 individuals per species, two cores per individual, and each core is composed of several segments. Here's the dput output from a subset of my data frame:
all <- structure(list(Sp_ind_core = structure(c(82L, 150L, 127L, 103L,
30L, 25L, 112L, 159L, 170L, 27L), .Label = c("APEIME_1_a", "APEIME_1_b",
"APEIME_2_a", "APEIME_2_b", "APEIME_3_a", "APEIME_3_b", "APEIME_4_a",
"APEIME_4_b", "APEIME_5_a", "APEIME_5_b", "APEIME_6_a", "APEIME_6_b",
"ASTRGR_1_a", "ASTRGR_1_b", "ASTRGR_2_a", "ASTRGR_2_b", "ASTRGR_3_a",
"ASTRGR_3_b", "ASTRGR_4_a", "ASTRGR_4_b", "ASTRGR_5_a", "ASTRGR_5_b",
"ASTRGR_6_a", "ASTRGR_6_b", "CALOLO_1_a", "CALOLO_1_b", "CALOLO_2_a",
"CALOLO_2_b", "CALOLO_3_a", "CALOLO_3_b", "CALOLO_4_a", "CALOLO_4_b",
"CALOLO_5_a", "CALOLO_5_b", "CALOLO_6_a", "CALOLO_6_b", "CORDBI_1_a",
"CORDBI_1_b", "CORDBI_2_a", "CORDBI_2_b", "CORDBI_3_a", "CORDBI_3_b",
"CORDBI_4_a", "CORDBI_4_b", "CORDBI_6_a", "CORDBI_6_b", "CORDBI_7_a",
"CORDBI_7_b", "GUAZUL_1_a", "GUAZUL_1_b", "GUAZUL_2_a", "GUAZUL_2_b",
"GUAZUL_3_a", "GUAZUL_3_b", "GUAZUL_4_a", "GUAZUL_4_b", "GUAZUL_5_a",
"GUAZUL_5_b", "GUAZUL_6_a", "GUAZUL_6_b", "INGAMA_1_a", "INGAMA_1_b",
"INGAMA_2_a", "INGAMA_2_b", "INGAMA_3_a", "INGAMA_3_b", "INGAMA_4_a",
"INGAMA_4_b", "INGAMA_5_a", "INGAMA_5_b", "INGAMA_6_a", "INGAMA_6_b",
"JACACO_1_a", "JACACO_1_b", "JACACO_2_a", "JACACO_2_b", "JACACO_3_a",
"JACACO_3_b", "JACACO_4_a", "JACACO_4_b", "JACACO_6_a", "JACACO_6_b",
"JACACO_8_a", "JACACO_8_b", "LACMPA_1_a", "LACMPA_1_b", "LACMPA_2_a",
"LACMPA_2_b", "LACMPA_3_a", "LACMPA_3_b", "LACMPA_4_a", "LACMPA_4_b",
"LACMPA_5_a", "LACMPA_5_b", "LACMPA_6_a", "LACMPA_6_b", "LUEHSE_1_a",
"LUEHSE_1_b", "LUEHSE_2_a", "LUEHSE_2_b", "LUEHSE_3_a", "LUEHSE_3_b",
"LUEHSE_4_a", "LUEHSE_4_b", "LUEHSE_5_a", "LUEHSE_5_b", "LUEHSE_6_a",
"LUEHSE_6_b", "PLATEL_1_a", "PLATEL_1_b", "PLATEL_2_a", "PLATEL_2_b",
"PLATEL_3_a", "PLATEL_3_b", "PLATEL_4_a", "PLATEL_4_b", "PLATEL_5_a",
"PLATEL_5_b", "PLATEL_6_a", "PLATEL_6_b", "POULAR_1_a", "POULAR_1_b",
"POULAR_2_a", "POULAR_2_b", "POULAR_3_a", "POULAR_3_b", "POULAR_4_a",
"POULAR_4_b", "POULAR_5_a", "POULAR_5_b", "POULAR_6_a", "POULAR_6_b",
"PRIOCO_1_a", "PRIOCO_1_b", "PRIOCO_2_a", "PRIOCO_2_b", "PRIOCO_3_a",
"PRIOCO_3_b", "PRIOCO_4_a", "PRIOCO_4_b", "PRIOCO_5_a", "PRIOCO_5_b",
"PRIOCO_6_a", "PRIOCO_6_b", "SIMAAM_1_a", "SIMAAM_1_b", "SIMAAM_2_a",
"SIMAAM_2_b", "SIMAAM_5_a", "SIMAAM_5_d", "SIMAAM_6_a", "SIMAAM_6_b",
"SIMAAM_7_a", "SIMAAM_7_b", "SIMAAM_8_a", "SIMAAM_8_b", "SPONRA_1_a",
"SPONRA_1_b", "SPONRA_2_a", "SPONRA_2_b", "SPONRA_3_a", "SPONRA_3_b",
"SPONRA_4_a", "SPONRA_4_b", "SPONRA_5_a", "SPONRA_5_b", "SPONRA_7_a",
"SPONRA_7_b", "TABAGU_1_a", "TABAGU_1_b", "TABAGU_2_a", "TABAGU_2_b",
"TABAGU_3_a", "TABAGU_3_b", "TABAGU_4_a", "TABAGU_4_b", "TABAGU_5_a",
"TABAGU_5_b", "TABAGU_6_a", "TABAGU_6_b", "TABARO_1_a", "TABARO_1_b",
"TABARO_2_a", "TABARO_2_b", "TABARO_3_a", "TABARO_3_b", "TABARO_4_a",
"TABARO_4_b", "TABARO_5_a", "TABARO_5_b", "TABARO_6_a", "TABARO_6_b",
"TACHVE_1_a", "TACHVE_1_b", "TACHVE_3_a", "TACHVE_3_b", "TACHVE_4_a",
"TACHVE_4_b", "TACHVE_5_a", "TACHVE_5_b", "TACHVE_6_a", "TACHVE_6_b",
"TACHVE_7_a", "TACHVE_7_b", "TRIPCU_1_a", "TRIPCU_1_b", "TRIPCU_2_a",
"TRIPCU_2_b", "TRIPCU_3_a", "TRIPCU_3_b", "TRIPCU_4_a", "TRIPCU_4_b",
"TRIPCU_5_a", "TRIPCU_5_b", "TRIPCU_6_a", "TRIPCU_6_b", "VIROSE_1_a",
"VIROSE_1_b", "VIROSE_2_a", "VIROSE_2_b", "VIROSE_3_a", "VIROSE_3_b",
"VIROSE_4_a", "VIROSE_4_b", "VIROSE_5_a", "VIROSE_5_b", "VIROSE_6_a",
"VIROSE_6_b", "ZANTEK_1_a", "ZANTEK_1_b", "ZANTEK_2_a", "ZANTEK_2_b",
"ZANTEK_3_a", "ZANTEK_3_b", "ZANTEK_4_a", "ZANTEK_4_b", "ZANTEK_5_a",
"ZANTEK_5_b", "ZANTEK_6_a", "ZANTEK_6_b"), class = "factor"),
Species = structure(c(7L, 13L, 11L, 9L, 3L, 3L, 10L, 14L,
15L, 3L), .Label = c("APEIME", "ASTRGR", "CALOLO", "CORDBI",
"GUAZUL", "INGAMA", "JACACO", "LACMPA", "LUEHSE", "PLATEL",
"POULAR", "PRIOCO", "SIMAAM", "SPONRA", "TABAGU", "TABARO",
"TACHVE", "TRIPCU", "VIROSE", "ZANTEK"), class = "factor"),
Individual = c(6L, 5L, 4L, 4L, 3L, 1L, 2L, 2L, 1L, 2L), Core = structure(c(2L,
3L, 1L, 1L, 2L, 1L, 2L, 1L, 2L, 1L), .Label = c("a", "b",
"d"), class = "factor"), DP = c(3, 18, 5, 1, 9, 9, 1, 2,
1.6, 8), WD = c(0.469094923, 0.542043984, 0.218045113, 0.6192,
0.428431373, 0.542285345, 0.723011364, 0.282582938, 0.942528736,
0.58), DP.max = c(21, 18, 10, 15, 20, 22, 17, 18, 8.8, 8),
DPscale = c(0.142857143, 1, 0.5, 0.066666667, 0.45, 0.409090909,
0.058823529, 0.111111111, 0.181818182, 1), DP2 = c(0.020408163,
1, 0.25, 0.004444444, 0.2025, 0.167355372, 0.003460208, 0.012345679,
0.033057851, 1)), .Names = c("Sp_ind_core", "Species", "Individual",
"Core", "DP", "WD", "DP.max", "DPscale", "DP2"), row.names = c(1370L,
2480L, 2085L, 1730L, 474L, 390L, 1871L, 2619L, 2787L, 434L), class = "data.frame")
There are 20 panels in the multipanel plots, one for each species, with density on the y axis and pith distance on the x axis. Each panel plot shows 12 different lines for the 12 cores for that species. Here's the code for the multipanel plot:
library(ggplot2)
ggplot(all, aes(x=DPscale, y=WD, group=Sp_ind_core))+ geom_line()
+theme_bw()+facet_wrap(facets=~Species, ncol=5)
The challenge for me is to add a line of the fitted values from a linear mixed effects model using the lme function in the nlme package in R. Here's the code for mixed effects model for each species in the subset data:
library(nlme)
jacaco.qm<-lme(WD~DPscale+DP2, random=~1|Individual/Core, jacaco)
luehse.m<-lme(WD~1, random=~1|Individual/Core, luehse)
simaam.lm<-lme(WD~DPscale, random=~1|Individual/Core, simaam)
calolo.qm<-lme(WD~DPscale+DP2, random=~1|Individual/Core, calolo)
platel.m<-lme(WD~1, random=~1|Individual/Core, platel)
tabagu.m<-lme(WD~1, random=~1|Individual/Core, tabagu)
sponra.lm<-lme(WD~DPscale, random=~1|Individual/Core, sponra)
There is a relevant question here:
Plotting results of lme4 with ggplot2
which leads me to think that the solution is to create a second data frame containing the fixed effects coefficients for each facet in the multipanel plot, and then use the geom_abline function to draw the fitted lines based on the fixed effects coefficient (intercept and slope(s)).
However, I am unable to apply this solution because the relationship between density and pith distance varies for the 20 species. Some of the species have a random effects model, others have a linear model, and yet others have quadratic models, meaning there are 1, 2, and 3 fixed effects coefficients depending on the model. I tried to create a data frame containing the fixed effects for all species:
fixed.m1 <- data.frame(fixef(jacaco.qm), fixef(luehse.m), fixef(platel.m), fixef(poular.qm),
fixef(simaam.lm), fixef(calolo.qm), fixef(tabagu.m), fixef(sponra.lm))
I get the following error:
Error in data.frame(fixef(apeime.lm), fixef(astrgr.qm), fixef(calolo.qm), : arguments imply
differing number of rows: 2, 3
Lastly, I have also tried the stat_smooth function in ggplot2 before for the purpose of adding fitted lines, and although it worked by adding fitted lines to each panel in the plot, it was an incorrect solution because the stat_smooth function doesn't have "lme" as one of its methods, and lines fitted by lm are not the same as the lines based on the fitted values from the lme models.
In summary, my problem is that I am unable to add lines based on fitted values from lme models to multipanel plots drawn by ggplot because the species in each multipanel plot have different numbers of fixed effects coefficients depending on whether they have a mean, linear or quadratic model. Specifically, I want fitted lines for species best fit by linear and quadratic models, and no lines for species best fit by mean model.
Thanks in advance.
