Model Selection
Jocelyn Stalker
10/5/2021
Reptile Diversity in Agricultural Landscapes
pairs(div.data[,c(3,5:11)], lower.panel = NULL)
mod1 <- lm(diversity_index ~ ai + human_settlement, data= div.data)
mod2 <- lm(diversity_index ~ agriculture + ai, data= div.data)
mod3 <- lm(diversity_index ~ ai + natural_vegetation, data= div.data)
mod4 <- lm(diversity_index ~ natural_vegetation + seminatural_vegetation, data= div.data)
mod5 <- lm(diversity_index ~ agriculture + human_settlement, data= div.data)mod.sel <- model.sel(mod1,mod2,mod3,mod4,mod5)
mod.sel## Model selection table
## (Int) ai hmn_stt agr ntr_vgt smn_vgt
## mod5 1.501 -1.393e-06 1.972e-07
## mod1 1.443 0.002057 -3.153e-07
## mod3 1.462 0.001535 -2.839e-08
## mod2 1.441 0.001731 1.494e-08
## mod4 1.532 -5.522e-08 6.042e-08
## family df logLik AICc delta weight
## mod5 gaussian(identity) 4 -6.377 21.6 0.00 0.300
## mod1 gaussian(identity) 4 -6.757 22.4 0.76 0.205
## mod3 gaussian(identity) 4 -6.849 22.5 0.94 0.187
## mod2 gaussian(identity) 4 -6.914 22.7 1.08 0.175
## mod4 gaussian(identity) 4 -7.185 23.2 1.62 0.134
## Models ranked by AICc(x)subset(mod.sel, cumsum(mod.sel$weight) <= .95)## Model selection table
## (Int) ai hmn_stt agr ntr_vgt family df
## mod5 1.501 -1.393e-06 1.972e-07 gaussian(identity) 4
## mod1 1.443 0.002057 -3.153e-07 gaussian(identity) 4
## mod3 1.462 0.001535 -2.839e-08 gaussian(identity) 4
## mod2 1.441 0.001731 1.494e-08 gaussian(identity) 4
## logLik AICc delta weight
## mod5 -6.377 21.6 0.00 0.346
## mod1 -6.757 22.4 0.76 0.236
## mod3 -6.849 22.5 0.94 0.216
## mod2 -6.914 22.7 1.08 0.202
## Models ranked by AICc(x)# Consistent AIC with Fishers information matrix
model.sel(mod1,mod2,mod3,mod4,mod5, rank = CAICF) ## Model selection table
## (Int) ai hmn_stt agr ntr_vgt smn_vgt
## mod1 1.443 0.002057 -3.153e-07
## mod2 1.441 0.001731 1.494e-08
## mod3 1.462 0.001535 -2.839e-08
## mod5 1.501 -1.393e-06 1.972e-07
## mod4 1.532 -5.522e-08 6.042e-08
## family df logLik CAICF delta weight
## mod1 gaussian(identity) 4 -6.757 85.5 0.00 0.808
## mod2 gaussian(identity) 4 -6.914 89.8 4.25 0.097
## mod3 gaussian(identity) 4 -6.849 89.8 4.26 0.096
## mod5 gaussian(identity) 4 -6.377 103.7 18.15 0.000
## mod4 gaussian(identity) 4 -7.185 108.9 23.38 0.000
## Models ranked by CAICF(x)importance(mod.sel)## ai human_settlement agriculture natural_vegetation
## Sum of weights: 0.57 0.50 0.47 0.32
## N containing models: 3 2 2 2
## seminatural_vegetation
## Sum of weights: 0.13
## N containing models: 1# Model average using all candidate models, always use revised.var = TRUE
mod.avg.ests <- model.avg(mod.sel, revised.var = TRUE)
mod.avg.ests##
## Call:
## model.avg(object = mod.sel, revised.var = TRUE)
##
## Component models:
## '23' '12' '14' '13' '45'
##
## Coefficients:
## (Intercept) agriculture human_settlement ai natural_vegetation
## full 1.475424 6.171814e-08 -4.819682e-07 0.001011250 -1.268253e-08
## subset 1.475424 1.299910e-07 -9.554144e-07 0.001784363 -3.957664e-08
## seminatural_vegetation
## full 8.069682e-09
## subset 6.041603e-08anova(mod1)## Analysis of Variance Table
##
## Response: diversity_index
## Df Sum Sq Mean Sq F value Pr(>F)
## ai 1 0.1099 0.109862 1.3634 0.2486
## human_settlement 1 0.0270 0.027004 0.3351 0.5653
## Residuals 49 3.9482 0.080576f1 <- ggplot(div.data, aes(ai, diversity_index)) +
geom_point() +
geom_smooth(method="lm") + labs(title= "Reptile Diversity as a Function of \n Agricultural Intensity", x= "Agricultural Intensity", y= "Reptile Diversity Index")
f2 <- ggplot(div.data, aes(human_settlement, diversity_index)) +
geom_point() +
geom_smooth(method="lm") + labs(title= "Reptile Diversity as a Function of \n Nearby Human Settlements", x= "Human Settlements", y= "Reptile Diversity Index")
f1+f2## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
Drivers of Density of an Urban Lizard
Tropidurus hispidus
pairs(lagartixas[,9:13], lower.panel = NULL)
# arthropods and herb cover look suspicious, so let's not include them in the same modelm1 <- lm(encounter_rate_lagartixas_ind_m ~ trees_m + mean_herbcover, data= lagartixas)
m2 <- lm(encounter_rate_lagartixas_ind_m ~ mean_impervious + mean_shelters, data= lagartixas)
m3 <- lm(encounter_rate_lagartixas_ind_m ~ mean_arthropods + mean_shelters, data= lagartixas)
m4 <- lm(encounter_rate_lagartixas_ind_m ~ trees_m + mean_arthropods, data= lagartixas)
m5 <- lm(encounter_rate_lagartixas_ind_m ~ mean_impervious + mean_herbcover, data= lagartixas)mod.sel2 <- model.sel(m1,m2,m3,m4,m5)
mod.sel2## Model selection table
## (Int) men_hrb trs_m men_imp men_shl men_art family
## m3 0.005113 0.004326 -1.582e-04 gaussian(identity)
## m2 0.003786 9.156e-06 0.004311 gaussian(identity)
## m1 0.003289 4.597e-05 0.08305 gaussian(identity)
## m4 0.004543 0.08213 2.235e-05 gaussian(identity)
## m5 -0.002826 2.172e-04 1.247e-04 gaussian(identity)
## df logLik AICc delta weight
## m3 4 167.520 -326.1 0.00 0.462
## m2 4 167.420 -325.9 0.20 0.418
## m1 4 165.708 -322.5 3.62 0.075
## m4 4 165.136 -321.3 4.77 0.043
## m5 4 162.064 -315.2 10.91 0.002
## Models ranked by AICc(x)m6 <- lm(encounter_rate_lagartixas_ind_m ~ mean_shelters + mean_impervious + mean_arthropods, data= lagartixas)
m7 <- lm(encounter_rate_lagartixas_ind_m ~ mean_shelters + mean_impervious + mean_arthropods + trees_m, data= lagartixas)
mod.sel3 <- model.sel(m1,m2,m3,m4,m5,m6,m7)
mod.sel3## Model selection table
## (Int) men_hrb trs_m men_imp men_shl men_art
## m7 0.003083 0.05193 -5.397e-06 0.003360 -8.286e-05
## m3 0.005113 0.004326 -1.582e-04
## m2 0.003786 9.156e-06 0.004311
## m6 0.005265 -1.677e-06 0.004329 -1.713e-04
## m1 0.003289 4.597e-05 0.08305
## m4 0.004543 0.08213 2.235e-05
## m5 -0.002826 2.172e-04 1.247e-04
## family df logLik AICc delta weight
## m7 gaussian(identity) 6 170.128 -326.2 0.00 0.297
## m3 gaussian(identity) 4 167.520 -326.1 0.07 0.287
## m2 gaussian(identity) 4 167.420 -325.9 0.27 0.260
## m6 gaussian(identity) 5 167.522 -323.6 2.58 0.082
## m1 gaussian(identity) 4 165.708 -322.5 3.69 0.047
## m4 gaussian(identity) 4 165.136 -321.3 4.84 0.026
## m5 gaussian(identity) 4 162.064 -315.2 10.98 0.001
## Models ranked by AICc(x)subset(mod.sel3, delta <2)## Model selection table
## (Int) trs_m men_imp men_shl men_art family df
## m7 0.003083 0.05193 -5.397e-06 0.003360 -8.286e-05 gaussian(identity) 6
## m3 0.005113 0.004326 -1.582e-04 gaussian(identity) 4
## m2 0.003786 9.156e-06 0.004311 gaussian(identity) 4
## logLik AICc delta weight
## m7 170.128 -326.2 0.00 0.352
## m3 167.520 -326.1 0.07 0.340
## m2 167.420 -325.9 0.27 0.308
## Models ranked by AICc(x)anova(m3)## Analysis of Variance Table
##
## Response: encounter_rate_lagartixas_ind_m
## Df Sum Sq Mean Sq F value Pr(>F)
## mean_arthropods 1 0.00001788 0.00001788 0.3566 0.5535
## mean_shelters 1 0.00093566 0.00093566 18.6565 8.78e-05 ***
## Residuals 44 0.00220669 0.00005015
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m2)## Analysis of Variance Table
##
## Response: encounter_rate_lagartixas_ind_m
## Df Sum Sq Mean Sq F value Pr(>F)
## mean_impervious 1 0.00001777 0.00001777 0.3528 0.5556
## mean_shelters 1 0.00092632 0.00092632 18.3915 9.68e-05 ***
## Residuals 44 0.00221614 0.00005037
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m7)## Analysis of Variance Table
##
## Response: encounter_rate_lagartixas_ind_m
## Df Sum Sq Mean Sq F value Pr(>F)
## mean_shelters 1 0.00093741 0.00093741 19.9354 5.923e-05 ***
## mean_impervious 1 0.00000668 0.00000668 0.1421 0.70807
## mean_arthropods 1 0.00000956 0.00000956 0.2034 0.65431
## trees_m 1 0.00023164 0.00023164 4.9262 0.03191 *
## Residuals 42 0.00197494 0.00004702
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1f3 <- ggplot(lagartixas, aes(mean_shelters, encounter_rate_lagartixas_ind_m)) +
geom_point() +
geom_smooth(method="lm") + labs(title= "The Effect of Shelters on \nLagartixa Encounter Rate", x= "Mean # of Shelters", y= "Lagartixa Encounter Rate (individuals/meter)")
f4 <- ggplot(lagartixas, aes(trees_m, encounter_rate_lagartixas_ind_m)) +
geom_point() +
geom_smooth(method="lm") + labs(title= "The Effect of Trees on \nLagartixa Encounter Rate", x= "Trees per Meter", y= "Lagartixa Encounter Rate (individuals/meter)")
f3+f4## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 1 rows containing non-finite values (stat_smooth).## Warning: Removed 1 rows containing missing values (geom_point).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 1 rows containing non-finite values (stat_smooth).
## Warning: Removed 1 rows containing missing values (geom_point).