vif(lm(kg.ha ~ C3L + C3D + C4L + C4D + FbL + FbD +
WDY + LTR + LAI + VOR + UPRIGHT + LTRDEPTH +
intPAR, PARDat))
## Warning in summary.lm(object, ...): essentially perfect fit: summary may be
## unreliable
## C3L C3D C4L C4D FbL FbD WDY LTR
## 4.876125 3.959942 2.173508 4.037171 2.509100 1.926772 1.737838 3.592892
## LAI VOR UPRIGHT LTRDEPTH intPAR
## 9.882463 5.055673 1.678711 1.812821 5.154836
# Fit full combination of variables
m.null <- lm(kg.ha ~ 1, PARDat)
m.1a <- lm(kg.ha ~ intPAR, PARDat)
m.1b <- lm(kg.ha ~ LTR, PARDat)
m.1c <- lm(kg.ha ~ C3L, PARDat)
m.1d <- lm(kg.ha ~ C3D, PARDat)
m.1e <- lm(kg.ha ~ VOR, PARDat)
m.1f <- lm(kg.ha ~ LTRDEPTH, PARDat)
m.2a <- lm(kg.ha ~ intPAR + LTR, PARDat)
m.2b <- lm(kg.ha ~ intPAR + C3L, PARDat)
m.2c <- lm(kg.ha ~ intPAR + C3D, PARDat)
m.2d <- lm(kg.ha ~ intPAR + LTRDEPTH, PARDat)
m.3a <- lm(kg.ha ~ intPAR + C3L + C3D, PARDat)
m.3b <- lm(kg.ha ~ C3D + C3L + LTR, PARDat)
m.3c <- lm(kg.ha ~ intPAR + LTR + VOR, PARDat)
m.3d <- lm(kg.ha ~ intPAR + LTRDEPTH + VOR, PARDat)
m.4a <- lm(kg.ha ~ C3D + C3L + LTR +intPAR, PARDat)
m.5a <- lm(kg.ha ~ LAI + LTR, PARDat)
m.5b <- lm(kg.ha ~ LAI + C3L, PARDat)
m.5c <- lm(kg.ha ~ LAI + C3D, PARDat)
m.5d <- lm(kg.ha ~ LAI + LTRDEPTH, PARDat)
m.5e <- lm(kg.ha ~ LAI + C3L + C3D, PARDat)
m.5f <- lm(kg.ha ~ LAI + intPAR, PARDat)
cand.mod.names <- c("m.null", "m.1a", "m.1b", "m.1c",
"m.1d", "m.1e", "m.1f", "m.2a", "m.2b", "m.2c", "m.2d", "m.3a", "m.3b", "m.3c", "m.3d", "m.4a", "m.5a", "m.5b", "m.5c", "m.5d", "m.5e", "m.5f")
cand.mods <- list( )
for(i in 1:length(cand.mod.names)) {
cand.mods[[i]] <- get(cand.mod.names[i]) }
print(aictab(cand.set = cand.mods,
modnames = cand.mod.names))
##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## m.3c 5 1315.14 0.00 0.99 0.99 -652.16
## m.5a 4 1325.78 10.64 0.00 1.00 -658.62
## m.4a 6 1333.07 17.93 0.00 1.00 -659.95
## m.2a 4 1335.85 20.72 0.00 1.00 -663.66
## m.3b 5 1372.10 56.96 0.00 1.00 -680.64
## m.3d 5 1379.76 64.62 0.00 1.00 -684.47
## m.5d 4 1381.58 66.45 0.00 1.00 -686.52
## m.1b 3 1382.24 67.10 0.00 1.00 -687.96
## m.5b 4 1384.70 69.56 0.00 1.00 -688.08
## m.5e 5 1386.77 71.63 0.00 1.00 -687.97
## m.5f 4 1392.83 77.69 0.00 1.00 -692.14
## m.5c 4 1399.22 84.08 0.00 1.00 -695.34
## m.3a 5 1399.76 84.63 0.00 1.00 -694.47
## m.2b 4 1399.93 84.79 0.00 1.00 -695.69
## m.2d 4 1401.39 86.26 0.00 1.00 -696.43
## m.2c 4 1411.77 96.63 0.00 1.00 -701.61
## m.1a 3 1413.70 98.56 0.00 1.00 -703.69
## m.1e 3 1428.63 113.49 0.00 1.00 -711.15
## m.1f 3 1444.38 129.24 0.00 1.00 -719.03
## m.1c 3 1444.94 129.80 0.00 1.00 -719.31
## m.1d 3 1477.88 162.74 0.00 1.00 -735.78
## m.null 2 1492.05 176.91 0.00 1.00 -743.94
mod.tab <- aictab(cand.set = cand.mods,
modnames = cand.mod.names)
mod.tab <- as.data.frame(mod.tab)
anova(m.5a, m.3c)
## Analysis of Variance Table
##
## Model 1: kg.ha ~ LAI + LTR
## Model 2: kg.ha ~ intPAR + LTR + VOR
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 76 80640331
## 2 75 68468923 1 12171408 13.332 0.0004802 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(m.4a, m.2a)
## Analysis of Variance Table
##
## Model 1: kg.ha ~ C3D + C3L + LTR + intPAR
## Model 2: kg.ha ~ intPAR + LTR
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 74 83400885
## 2 76 91604671 -2 -8203786 3.6395 0.03107 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(m.5a.CIs <- as.data.frame(confint(m.5a) ) )
## 2.5 % 97.5 %
## (Intercept) 1149.66468 2068.29382
## LAI 868.20513 1351.30862
## LTR 21.20549 30.79355
(m.3c.CIs <- as.data.frame(confint(m.3c) ) )
## 2.5 % 97.5 %
## (Intercept) 330.5540133 1351.515978
## intPAR 0.9131481 2.402777
## LTR 22.4244672 31.026326
## VOR 243.9830829 563.531078
terms <- c("(Intercept)", "intPAR", "LAI", "LTR", "C3L", "C3D", "VOR", "LTRDEPTH")
av.params <- as.data.frame(array(NA,c(length(terms),4)))
colnames(av.params)<-c("term","estimate","ciL","ciU")
for(i in 1:length(terms)) {
av <- modavg(parm = paste(terms[i]),
cand.set = cand.mods,
modnames = cand.mod.names)
av.params[i,1] <- terms[i]
av.params[i,2] <- round(av$Mod.avg.beta, 2)
av.params[i,3] <- round(av$Lower.CL, 3)
av.params[i,4] <- round(av$Upper.CL, 3) }
av.params
## term estimate ciL ciU
## 1 (Intercept) 844.82 331.941 1357.693
## 2 intPAR 1.66 0.925 2.391
## 3 LAI 1109.76 872.051 1347.463
## 4 LTR 26.72 22.487 30.958
## 5 C3L -1.80 -10.082 6.485
## 6 C3D 7.04 1.925 12.150
## 7 VOR 403.76 246.561 560.953
## 8 LTRDEPTH 195.55 112.068 279.033
View(av.params)
av.params$term <- plyr::revalue(av.params$term,
c("intPAR" = "Intercepted\nPAR",
"LAI" = "Leaf Area Index",
"LTR" = "Actual Litter\nBiomass",
"C3L" = "Live C3\ngrass",
"C3D" = "Dead C3\ngrass",
"VOR" = "Visual\nObstruction\nReading",
"LTRDEPTH" = "Litter\nDepth"))
cept.CI <- ggplot(subset(av.params, terms != "(Intercept)")) +
coord_flip() + theme_bw(16) +
geom_hline(yintercept = 0) +
geom_errorbar(aes(x=term,
ymin=ciL, ymax=ciU),
width=0.2, size=2, color="#377eb8") +
geom_point(aes(x=term, y=estimate), size=7, pch=21,
bg="#377eb8", color="white", stroke=2) +
labs(y = "regression parameter estimate (95% CI)",
x = " ")+
theme(axis.text = element_text(color = "black",
size = 16))
cept.CI
