This box continues with the forest bird data set used in Boxes 8.2 and 8.10
We calculated common measures of relative importance based on the fit of the following full model relating bird abundance to five predictors for the dataset from Loyn (1997; see Box 8.2): (bird abundance)i = β0 + β1(log patch area) + β2(log distance to nearest patch) +β3(grazing) +β4(altitude) + β5(year isolated) + εi
First, load the required packages (relaimpo, car, hier.part, MuMIn, lm.beta, sjstats)
Import loyn data file (data/loyn.csv); from previous example, we needed to transform area and distance
loyn <- read.csv("../data/loyn.csv")
head(loyn,10)loyn$logarea <- log10(loyn$area)
loyn$logdist <- log10(loyn$dist)For some analyses, we also need two dataframes that are subsets of loyn containing the response variable and the predictors
loyn_abund<-loyn$abund
loyn_pred<-subset(loyn, select = c("logarea","logdist","alt","yearisol","graze"))loyn.lm <- lm(abund~ logarea+logdist+graze+alt+yearisol, data=loyn)
tidy(loyn.lm)lm.beta.loyn <- lm.beta(loyn.lm)
tidy(lm.beta.loyn, conf.int=TRUE)std.coef(loyn.lm, partial.sd=FALSE) Estimate* Std. Error* df
(Intercept) 0.000000 0.000000 50
logarea 0.552111 0.101112 50
logdist -0.050237 0.089388 50
graze -0.229600 0.126203 50
alt 0.086970 0.092826 50
yearisol 0.182493 0.104669 50std.coef(loyn.lm, partial.sd=TRUE) Estimate* Std. Error* df
(Intercept) 0.00000 0.00000 50
logarea 4.61565 0.84529 50
logdist -0.47507 0.84529 50
graze -1.53783 0.84529 50
alt 0.79197 0.84529 50
yearisol 1.47379 0.84529 50calc.relimp(loyn.lm, type = c("lmg", "pmvd", "last", "first", "betasq", "pratt"), rela=FALSE)Response variable: abund
Total response variance: 115.2351
Analysis based on 56 observations
5 Regressors:
logarea logdist graze alt yearisol
Proportion of variance explained by model: 68.43%
Metrics are not normalized (rela=FALSE).
Relative importance metrics:
lmg pmvd last first betasq pratt
logarea 0.32413897 0.441235188 0.188237420 0.5476530 0.304826673 0.408581982
logdist 0.01039951 0.003176187 0.001994139 0.0160588 0.002523798 -0.006366253
graze 0.19421053 0.162925775 0.020895833 0.4658218 0.052716286 0.156704799
alt 0.05484119 0.013555615 0.005541867 0.1488696 0.007563750 0.033556104
yearisol 0.10074578 0.063443199 0.019191723 0.2533690 0.033303745 0.091859333
Average coefficients for different model sizes:
1X 2Xs 3Xs 4Xs 5Xs
logarea 9.77842156 8.68715779 7.93631294 7.49569857 7.29528918
logdist 3.28742481 1.82790369 0.48827974 -0.54907795 -1.30324642
graze -4.98129971 -4.21006477 -3.35571467 -2.47849706 -1.67573437
alt 0.09515173 0.06443007 0.04304481 0.02976254 0.02144778
yearisol 0.21122456 0.14501519 0.10324962 0.08419284 0.07657980Note: Available metrics from relaimpo package: lmg, pmvd (non-US version only), last, first, betasq, pratt, genizi and car. - from package help
Can only get pmvd to run if we obtain the non-US version from author’s website and install manually - copy into relevant R library folder.
Author web site is http://prof.beuth-hochschule.de/groemping/software/relaimpo/
Note: Use stand coefficents based on partial sd above instead of betasq
loyn.boot <- boot.relimp(loyn.lm, b=1000, type = c("lmg", "pmvd"))
booteval.relimp(loyn.boot)Response variable: abund
Total response variance: 115.2351
Analysis based on 56 observations
5 Regressors:
logarea logdist graze alt yearisol
Proportion of variance explained by model: 68.43%
Metrics are not normalized (rela=FALSE).
Relative importance metrics:
lmg pmvd
logarea 0.32413897 0.441235188
logdist 0.01039951 0.003176187
graze 0.19421053 0.162925775
alt 0.05484119 0.013555615
yearisol 0.10074578 0.063443199
Average coefficients for different model sizes:
1X 2Xs 3Xs 4Xs 5Xs
logarea 9.77842156 8.68715779 7.93631294 7.49569857 7.29528918
logdist 3.28742481 1.82790369 0.48827974 -0.54907795 -1.30324642
graze -4.98129971 -4.21006477 -3.35571467 -2.47849706 -1.67573437
alt 0.09515173 0.06443007 0.04304481 0.02976254 0.02144778
yearisol 0.21122456 0.14501519 0.10324962 0.08419284 0.07657980
Confidence interval information ( 1000 bootstrap replicates, bty= perc ):
Relative Contributions with confidence intervals:
Lower Upper
percentage 0.95 0.95 0.95
logarea.lmg 0.3241 AB___ 0.2152 0.4561
logdist.lmg 0.0104 ___DE 0.0049 0.0679
graze.lmg 0.1942 ABC__ 0.1112 0.2823
alt.lmg 0.0548 __CDE 0.0110 0.1338
yearisol.lmg 0.1007 _BCD_ 0.0300 0.2089
logarea.pmvd 0.4412 AB___ 0.1975 0.6373
logdist.pmvd 0.0032 __CDE 0.0000 0.0678
graze.pmvd 0.1629 ABCDE 0.0028 0.4705
alt.pmvd 0.0136 _BCDE 0.0000 0.1045
yearisol.pmvd 0.0634 _BCDE 0.0002 0.2738
Letters indicate the ranks covered by bootstrap CIs.
(Rank bootstrap confidence intervals always obtained by percentile method)
CAUTION: Bootstrap confidence intervals can be somewhat liberal.
Differences between Relative Contributions:
Lower Upper
difference 0.95 0.95 0.95
logarea-logdist.lmg 0.3137 * 0.1855 0.4432
logarea-graze.lmg 0.1299 -0.0421 0.3187
logarea-alt.lmg 0.2693 * 0.1113 0.4290
logarea-yearisol.lmg 0.2234 * 0.0389 0.3893
logdist-graze.lmg -0.1838 * -0.2640 -0.0779
logdist-alt.lmg -0.0444 -0.1162 0.0306
logdist-yearisol.lmg -0.0903 -0.1940 0.0051
graze-alt.lmg 0.1394 * 0.0419 0.2176
graze-yearisol.lmg 0.0935 -0.0224 0.1911
alt-yearisol.lmg -0.0459 -0.1648 0.0665
logarea-logdist.pmvd 0.4381 * 0.1889 0.6289
logarea-graze.pmvd 0.2783 -0.2508 0.6001
logarea-alt.pmvd 0.4277 * 0.1722 0.6175
logarea-yearisol.pmvd 0.3778 * 0.0860 0.5930
logdist-graze.pmvd -0.1597 -0.4490 0.0112
logdist-alt.pmvd -0.0104 -0.1010 0.0635
logdist-yearisol.pmvd -0.0603 -0.2617 0.0328
graze-alt.pmvd 0.1494 -0.0343 0.4325
graze-yearisol.pmvd 0.0995 -0.2123 0.4295
alt-yearisol.pmvd -0.0499 -0.2576 0.0813
* indicates that CI for difference does not include 0.
CAUTION: Bootstrap confidence intervals can be somewhat liberal. hier.part(loyn_abund, loyn_pred, family="gaussian", gof="Rsqu")Error in hier.part(loyn_abund, loyn_pred, family = "gaussian", gof = "Rsqu") :
could not find function "hier.part"The package hier.part was removed from CRAN in March 2023. The code above will work if you have hier.part installed already. An alternative is to use the package glmm.hp, which is done in the next code chunk.
Hier.part can also be installed from Github, though there may be issues with M1/M2 Macs. The quick way from Github is using devtools: devtools::install_github(“cjbwalsh/hier.part”)
library (glmm.hp)Loading required package: vegan
Loading required package: permute
This is vegan 2.6-4
Attaching package: ‘vegan’
The following object is masked from ‘package:survey’:
calibrateglmm.hp(loyn.lm, type="R2")Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.$Total.R2
[1] 0.684336
$hierarchical.partitioning
Unique Average.share Individual I.perc(%)
logarea 0.1882 0.1359 0.3241 47.37
logdist 0.0020 0.0084 0.0104 1.52
graze 0.0209 0.1733 0.1942 28.38
alt 0.0055 0.0493 0.0548 8.01
yearisol 0.0192 0.0815 0.1007 14.72
$variables
[1] "logarea" "logdist" "graze" "alt" "yearisol"
$type
[1] "hierarchical.partitioning"
attr(,"class")
[1] "glmmhp"options(na.action = "na.fail")
loyn.dredge <-dredge(loyn.lm, beta="none", evaluate=TRUE)Fixed term is "(Intercept)"loyn.dredgeGlobal model call: lm(formula = abund ~ logarea + logdist + graze + alt + yearisol,
data = loyn)
---
Model selection table
(Intrc) alt graze logar lgdst yersl df logLik AICc delta weight
23 -134.300 -1.902 7.166 0.07835 5 -180.555 372.3 0.00 0.211
24 -141.900 0.02586 -1.601 7.076 0.07991 6 -179.761 373.2 0.93 0.133
7 21.600 -2.854 6.890 4 -182.257 373.3 0.99 0.129
31 -120.500 -1.939 7.487 -2.0480 0.07354 6 -180.072 373.9 1.55 0.097
22 -236.700 0.03623 8.150 0.12480 5 -181.428 374.1 1.75 0.088
8 17.280 0.02468 -2.584 6.799 5 -181.577 374.4 2.05 0.076
15 26.630 -2.827 7.298 -2.4750 5 -181.582 374.4 2.06 0.076
21 -252.200 8.593 0.13520 4 -182.993 374.8 2.46 0.062
32 -131.800 0.02145 -1.676 7.295 -1.3030 0.07658 7 -179.585 375.5 3.19 0.043
16 22.110 0.01850 -2.632 7.126 -1.8480 6 -181.238 376.2 3.88 0.030
30 -233.700 0.03413 8.294 -0.6964 0.12420 6 -181.380 376.5 4.16 0.026
29 -241.900 8.904 -1.8350 0.13190 5 -182.638 376.5 4.17 0.026
6 3.959 0.04862 9.062 4 -187.342 383.5 11.16 0.001
14 6.929 0.04490 9.300 -1.1950 5 -187.226 385.7 13.34 0.000
5 10.400 9.778 3 -189.659 385.8 13.47 0.000
13 16.130 10.200 -2.7660 4 -189.012 386.8 14.50 0.000
3 34.370 -4.981 3 -194.315 395.1 22.78 0.000
4 28.560 0.03191 -4.597 4 -193.573 395.9 23.62 0.000
19 -62.750 -4.440 0.04898 4 -193.886 396.6 24.25 0.000
11 32.550 -4.950 0.7777 4 -194.268 397.3 25.01 0.000
20 -73.580 0.03285 -4.017 0.05143 5 -193.087 397.4 25.07 0.000
12 23.010 0.03793 -4.448 1.9060 5 -193.314 397.8 25.52 0.000
27 -72.190 -4.356 1.1420 0.05240 5 -193.786 398.8 26.46 0.000
28 -95.930 0.04056 -3.742 2.3960 0.05917 6 -192.681 399.1 26.77 0.000
26 -356.500 0.08167 5.3870 0.18060 5 -198.902 409.0 36.70 0.000
18 -348.500 0.07006 0.18350 4 -200.670 410.1 37.82 0.000
17 -392.300 0.21120 3 -203.690 413.8 41.53 0.000
25 -402.400 3.5440 0.21230 4 -202.981 414.7 42.44 0.000
10 -8.908 0.10710 5.7560 4 -205.772 420.3 48.02 0.000
2 5.598 0.09515 3 -207.358 421.2 48.87 0.000
1 19.510 2 -211.871 428.0 55.66 0.000
9 12.230 3.2870 3 -211.418 429.3 56.99 0.000
Models ranked by AICc(x) options(na.action = "na.fail")
loyn.dredge1 <-dredge(loyn.lm, beta="partial.sd", evaluate=TRUE)Fixed term is "(Intercept)"loyn.dredge1Global model call: lm(formula = abund ~ logarea + logdist + graze + alt + yearisol,
data = loyn)
---
Model selection table
(Intrc) alt graze logar lgdst yersl df logLik AICc delta weight
23 0 -1.834 4.749 1.5220 5 -180.555 372.3 0.00 0.211
24 0 1.0160 -1.485 4.681 1.5520 6 -179.761 373.2 0.93 0.133
7 0 -3.449 4.599 4 -182.257 373.3 0.99 0.129
31 0 -1.868 4.793 -0.7946 1.4190 6 -180.072 373.9 1.55 0.097
22 0 1.4810 6.141 2.9950 5 -181.428 374.1 1.75 0.088
8 0 0.9707 -2.961 4.529 5 -181.577 374.4 2.05 0.076
15 0 -3.415 4.689 -0.9671 5 -181.582 374.4 2.06 0.076
21 0 6.645 3.2930 4 -182.993 374.8 2.46 0.062
32 0 0.7920 -1.538 4.616 -0.4751 1.4740 7 -179.585 375.5 3.19 0.043
16 0 0.6850 -3.006 4.521 -0.6799 6 -181.238 376.2 3.88 0.030
30 0 1.3230 5.756 -0.2565 2.9740 6 -181.380 376.5 4.16 0.026
29 0 6.523 -0.7129 3.1910 5 -182.638 376.5 4.17 0.026
6 0 2.0160 7.014 4 -187.342 383.5 11.16 0.001
14 0 1.7580 6.640 -0.4412 5 -187.226 385.7 13.34 0.000
5 0 7.873 3 -189.659 385.8 13.47 0.000
13 0 7.832 -1.0810 4 -189.012 386.8 14.50 0.000
3 0 -7.261 3 -194.315 395.1 22.78 0.000
4 0 1.2570 -6.120 4 -193.573 395.9 23.62 0.000
19 0 -4.996 0.9587 4 -193.886 396.6 24.25 0.000
11 0 -7.142 0.3157 4 -194.268 397.3 25.01 0.000
20 0 1.2940 -4.242 1.0060 5 -193.087 397.4 25.07 0.000
12 0 1.4220 -5.718 0.7365 5 -193.314 397.8 25.52 0.000
27 0 -4.802 0.4588 1.0150 5 -193.786 398.8 26.46 0.000
28 0 1.5160 -3.767 0.9133 1.1420 6 -192.681 399.1 26.77 0.000
26 0 3.3420 2.1550 4.4500 5 -198.902 409.0 36.70 0.000
18 0 2.9390 4.5240 4 -200.670 410.1 37.82 0.000
17 0 5.3550 3 -203.690 413.8 41.53 0.000
25 0 1.4530 5.3820 4 -202.981 414.7 42.44 0.000
10 0 4.5090 2.3030 4 -205.772 420.3 48.02 0.000
2 0 4.1050 3 -207.358 421.2 48.87 0.000
1 0 2 -211.871 428.0 55.66 0.000
9 0 1.3480 3 -211.418 429.3 56.99 0.000
Models ranked by AICc(x) Note confint only provides CIs for conditional (natural) averages, based on z=1.96. Modified SE CIs in Box 9.1 calculated by hand for full averaging
loyn.ma<-model.avg(loyn.dredge)
summary(loyn.ma)
Call:
model.avg(object = loyn.dredge)
Component model call:
lm(formula = abund ~ <32 unique rhs>, data = loyn)
Component models:
df logLik AICc delta weight
235 5 -180.55 372.31 0.00 0.21
1235 6 -179.76 373.24 0.93 0.13
23 4 -182.26 373.30 0.99 0.13
2345 6 -180.07 373.86 1.55 0.10
135 5 -181.43 374.06 1.75 0.09
123 5 -181.58 374.35 2.05 0.08
234 5 -181.58 374.36 2.06 0.08
35 4 -182.99 374.77 2.46 0.06
12345 7 -179.59 375.50 3.19 0.04
1234 6 -181.24 376.19 3.88 0.03
1345 6 -181.38 376.47 4.16 0.03
345 5 -182.64 376.48 4.17 0.03
13 4 -187.34 383.47 11.16 0.00
134 5 -187.23 385.65 13.34 0.00
3 3 -189.66 385.78 13.47 0.00
34 4 -189.01 386.81 14.50 0.00
2 3 -194.31 395.09 22.78 0.00
12 4 -193.57 395.93 23.62 0.00
25 4 -193.89 396.56 24.25 0.00
24 4 -194.27 397.32 25.01 0.00
125 5 -193.09 397.37 25.07 0.00
124 5 -193.31 397.83 25.52 0.00
245 5 -193.79 398.77 26.46 0.00
1245 6 -192.68 399.08 26.77 0.00
145 5 -198.90 409.00 36.70 0.00
15 4 -200.67 410.13 37.82 0.00
5 3 -203.69 413.84 41.53 0.00
45 4 -202.98 414.75 42.44 0.00
14 4 -205.77 420.33 48.02 0.00
1 3 -207.36 421.18 48.87 0.00
(Null) 2 -211.87 427.97 55.66 0.00
4 3 -211.42 429.30 56.99 0.00
Term codes:
alt graze logarea logdist yearisol
1 2 3 4 5
Model-averaged coefficients:
(full average)
Estimate Std. Error Adjusted SE z value Pr(>|z|)
(Intercept) -106.85561 117.81334 118.76833 0.900 0.368
graze -1.73436 1.22630 1.23706 1.402 0.161
logarea 7.38974 1.38693 1.41472 5.223 2e-07 ***
yearisol 0.06376 0.05837 0.05885 1.083 0.279
alt 0.01094 0.01957 0.01980 0.553 0.580
logdist -0.56618 1.52862 1.55257 0.365 0.715
(conditional average)
Estimate Std. Error Adjusted SE z value Pr(>|z|)
(Intercept) -106.85561 117.81334 118.76833 0.900 0.3683
graze -2.17940 0.95905 0.97626 2.232 0.0256 *
logarea 7.38979 1.38678 1.41458 5.224 2e-07 ***
yearisol 0.09275 0.04761 0.04848 1.913 0.0557 .
alt 0.02750 0.02251 0.02302 1.195 0.2321
logdist -1.89107 2.30200 2.35492 0.803 0.4220
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1confint(loyn.ma) 2.5 % 97.5 %
(Intercept) -3.396373e+02 125.92604329
graze -4.092835e+00 -0.26596063
logarea 4.617275e+00 10.16231332
yearisol -2.262736e-03 0.18775937
alt -1.760897e-02 0.07261784
logdist -6.506628e+00 2.72449618loyn.ma1<-model.avg(loyn.dredge, beta="sd") #Recent update gives error with this averagingError in h(simpleError(msg, call)) :
error in evaluating the argument 'x' in selecting a method for function 't': incorrect number of dimensionssummary(loyn.ma1)Error in h(simpleError(msg, call)) :
error in evaluating the argument 'object' in selecting a method for function 'summary': object 'loyn.ma1' not foundconfint(loyn.ma1)Error: object 'loyn.ma1' not foundloyn.ma2<-model.avg(loyn.dredge, beta="partial.sd")
summary(loyn.ma2)
Call:
model.avg(object = get.models(object = loyn.dredge, subset = NA),
beta = "partial.sd")
Component model call:
lm(formula = abund ~ <32 unique rhs>, data = loyn)
Component models:
df logLik AICc delta weight
235 5 -180.55 372.31 0.00 0.21
1235 6 -179.76 373.24 0.93 0.13
23 4 -182.26 373.30 0.99 0.13
2345 6 -180.07 373.86 1.55 0.10
135 5 -181.43 374.06 1.75 0.09
123 5 -181.58 374.35 2.05 0.08
234 5 -181.58 374.36 2.06 0.08
35 4 -182.99 374.77 2.46 0.06
12345 7 -179.59 375.50 3.19 0.04
1234 6 -181.24 376.19 3.88 0.03
1345 6 -181.38 376.47 4.16 0.03
345 5 -182.64 376.48 4.17 0.03
13 4 -187.34 383.47 11.16 0.00
134 5 -187.23 385.65 13.34 0.00
3 3 -189.66 385.78 13.47 0.00
34 4 -189.01 386.81 14.50 0.00
2 3 -194.31 395.09 22.78 0.00
12 4 -193.57 395.93 23.62 0.00
25 4 -193.89 396.56 24.25 0.00
24 4 -194.27 397.32 25.01 0.00
125 5 -193.09 397.37 25.07 0.00
124 5 -193.31 397.83 25.52 0.00
245 5 -193.79 398.77 26.46 0.00
1245 6 -192.68 399.08 26.77 0.00
145 5 -198.90 409.00 36.70 0.00
15 4 -200.67 410.13 37.82 0.00
5 3 -203.69 413.84 41.53 0.00
45 4 -202.98 414.75 42.44 0.00
14 4 -205.77 420.33 48.02 0.00
1 3 -207.36 421.18 48.87 0.00
(Null) 2 -211.87 427.97 55.66 0.00
4 3 -211.42 429.30 56.99 0.00
Term codes:
alt graze logarea logdist yearisol
1 2 3 4 5
Model-averaged coefficients:
(full average)
Estimate Std. Error Adjusted SE z value Pr(>|z|)
(Intercept) 0.0000 0.0000 0.0000 NaN NaN
graze -1.8531 1.3946 1.4045 1.319 0.187
logarea 5.0078 1.0872 1.1032 4.539 5.6e-06 ***
yearisol 1.3600 1.3074 1.3166 1.033 0.302
alt 0.4314 0.7722 0.7813 0.552 0.581
logdist -0.2174 0.5855 0.5945 0.366 0.715
(conditional average)
Estimate Std. Error Adjusted SE z value Pr(>|z|)
(Intercept) 0.0000 0.0000 0.0000 NaN NaN
graze -2.3286 1.1561 1.1712 1.988 0.0468 *
logarea 5.0079 1.0871 1.1032 4.540 5.6e-06 ***
yearisol 1.9784 1.1239 1.1394 1.736 0.0825 .
alt 1.0841 0.8893 0.9090 1.193 0.2331
logdist -0.7261 0.8806 0.9007 0.806 0.4201
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1confint(loyn.ma2) 2.5 % 97.5 %
(Intercept) 0.0000000 0.00000000
graze -4.6239697 -0.03313715
logarea 2.8457176 7.16999469
yearisol -0.2547954 4.21154762
alt -0.6976252 2.86575899
logdist -2.4914401 1.03915465sw(loyn.dredge) logarea graze yearisol alt logdist
Sum of weights: 1.00 0.80 0.69 0.40 0.30
N containing models: 16 16 16 16 16 # importance(loyn.dredge) #this function now defunct, and sw recommendedcv_error(loyn, abund~ logarea+logdist+graze+alt+yearisol, k=5)Warning: `unnest()` has a new interface. See `?unnest` for details.
ℹ Try `df %>% unnest(c(predicted, residuals))`, with `mutate()` if needed.cv_error(loyn, abund~ logarea+logdist+graze+alt+yearisol, k=10)Warning: `unnest()` has a new interface. See `?unnest` for details.
ℹ Try `df %>% unnest(c(predicted, residuals))`, with `mutate()` if needed.cv_error(loyn, abund~ logarea+logdist+graze+alt+yearisol, k=56)Warning: `unnest()` has a new interface. See `?unnest` for details.
ℹ Try `df %>% unnest(c(predicted, residuals))`, with `mutate()` if needed.