Using emmeans to Estimate Marginal Means of spmodel Objects
Michael Dumelle, Matt Higham, and Jay M. Ver Hoef
Source:vignettes/articles/emmeans.Rmd
emmeans.RmdIntroduction
spmodel is an R package used to fit,
summarize, and predict for a variety of spatial statistical models. This
vignette demonstrates estimating marginal means (i.e., least-squares
means) of spmodel objects using the emmeans
R package (Lenth 2024).
Before proceeding, we load spmodel and emmeans
by running
If using spmodel in a formal publication or report,
please cite it. Citing spmodel lets us devote more
resources to the package in the future. We view the spmodel
citation by running
citation(package = "spmodel")#> To cite spmodel in publications use:
#>
#> Dumelle M, Higham M, Ver Hoef JM (2023). spmodel: Spatial statistical
#> modeling and prediction in R. PLOS ONE 18(3): e0282524.
#> https://doi.org/10.1371/journal.pone.0282524
#>
#> A BibTeX entry for LaTeX users is
#>
#> @Article{,
#> title = {{spmodel}: Spatial statistical modeling and prediction in {R}},
#> author = {Michael Dumelle and Matt Higham and Jay M. {Ver Hoef}},
#> journal = {PLOS ONE},
#> year = {2023},
#> volume = {18},
#> number = {3},
#> pages = {1--32},
#> doi = {10.1371/journal.pone.0282524},
#> url = {https://doi.org/10.1371/journal.pone.0282524},
#> }Applying emmeans to spmodel Objects
In this section, we use the point-referenced lake data,
an sf object that contains data on lake conductivty for
some southwestern states (Arizona, Colorado, Nevada, Utah) in the United
States. We view the first few rows of lake by running
lake#> Simple feature collection with 102 features and 9 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: -2004016 ymin: 1031593 xmax: -753669.4 ymax: 2338804
#> Projected CRS: NAD83 / Conus Albers
#> # A tibble: 102 × 10
#> comid log_cond cond state temp precip elev origin year
#> * <chr> <dbl> <dbl> <chr> <dbl> <dbl> <dbl> <chr> <fct>
#> 1 20451100 6.32 554 AZ 12.7 49.4 1567 HUMAN_MADE 2012
#> 2 20476542 7.02 1121 AZ 21.8 57.8 459 HUMAN_MADE 2012
#> 3 10001770 7.13 1246 AZ 23.2 9.12 69.1 HUMAN_MADE 2012
#> 4 20584396 6.17 477 AZ 11.2 44.4 1822 HUMAN_MADE 2012
#> 5 20524727 5.48 239 AZ 8.31 61.2 2168 HUMAN_MADE 2012
#> 6 20479908 7.00 1096 AZ 22.4 23.1 366. HUMAN_MADE 2012
#> 7 10001834 7.85 2570 AZ 23.5 9.89 57.7 NATURAL 2012
#> 8 20695686 4.77 118. AZ 9.13 59.1 2072. HUMAN_MADE 2012
#> 9 21327603 5.30 201 AZ 17.9 28.3 1006. HUMAN_MADE 2012
#> 10 20449310 4.33 76 AZ 9.79 61.4 1998. HUMAN_MADE 2012
#> # ℹ 92 more rows
#> # ℹ 1 more variable: geometry <POINT [m]>We can learn more about lake by running
help("lake", "spmodel"), and we can visualize the
distribution of log conductivity in lake by state and year
by running
ggplot(lake, aes(color = log_cond)) +
geom_sf() +
scale_color_viridis_c() +
theme_gray(base_size = 14)
Distribution of log conductivity in the lake data.
A Single-Factor Model
First we explore a single-factor model that characterizes the response variable, log conductivity, by each state (AZ, CO, NV, UT). We fit and summarize this model by running:
spmod1 <- splm(
formula = log_cond ~ state,
data = lake,
spcov_type = "exponential"
)
summary(spmod1)#>
#> Call:
#> splm(formula = log_cond ~ state, data = lake, spcov_type = "exponential")
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -3.1485 -0.9733 0.0490 0.8168 3.2232
#>
#> Coefficients (fixed):
#> Estimate Std. Error df t value Pr(>|t|)
#> (Intercept) 5.82803 0.41633 16.16319 13.999 1.86e-10 ***
#> stateCO -0.93746 0.60188 14.06127 -1.558 0.142
#> stateNV 0.08518 0.57571 17.80264 0.148 0.884
#> stateUT -0.02632 0.54839 12.49787 -0.048 0.962
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Pseudo R-squared: 0.03776
#>
#> Coefficients (exponential spatial covariance):
#> de ie range
#> 1.704e+00 1.704e-04 4.123e+04The summary() output provides mean estimates for each
state relative to the difference from a reference group (here,
AZ). Often, however, the question “What is the mean in each group?” is
of interest, and this is not straightforward to obtain from
summary(). Fortunately, emmeans makes this
information readily available via the emmeans function:
em11 <- emmeans(spmod1, ~ state)which, when printed, returns the mean estimates, standard errors, and confidence intervals for each factor level:
em11#> state emmean SE df lower.CL upper.CL
#> AZ 5.83 0.416 16.16 4.95 6.71
#> CO 4.89 0.435 12.44 3.95 5.83
#> NV 5.91 0.398 19.69 5.08 6.74
#> UT 5.80 0.357 9.07 4.99 6.61
#>
#> Degrees-of-freedom method: satterthwaite
#> Confidence level used: 0.95The emmeans object em11 is an
emmGrid object, but it can easily be converted into a data
frame using data.frame() or
tibble::as_tibble():
data.frame(em11)#> state emmean SE df lower.CL upper.CL
#> 1 AZ 5.828033 0.4163274 16.16319 4.946182 6.709884
#> 2 CO 4.890569 0.4346870 12.43790 3.947153 5.833986
#> 3 NV 5.913214 0.3977708 19.69159 5.082644 6.743783
#> 4 UT 5.801714 0.3571189 9.07055 4.994812 6.608616We then visualize the means and confidence intervals:
plot(em11)
Recall that summary() provides mean estimates for each
state relative to the difference from a reference group (i.e.,
contrasts with a reference group). Contrasts between mean estimates that
are not reference groups, however, is again not straightforward to
obtain from summary(). Fortunately, pairs()
provides a simple solution, creating contrasts for comparisons of each
factor level to all other factor levels that are easily visualized:
pairs(em11)#> contrast estimate SE df t.ratio p.value
#> AZ - CO 0.9375 0.602 14.1 1.558 0.4320
#> AZ - NV -0.0852 0.576 17.8 -0.148 0.9988
#> AZ - UT 0.0263 0.548 12.5 0.048 1.0000
#> CO - NV -1.0226 0.589 15.2 -1.736 0.3403
#> CO - UT -0.9111 0.562 11.1 -1.621 0.4065
#> NV - UT 0.1115 0.530 15.1 0.210 0.9966
#>
#> Degrees-of-freedom method: satterthwaite
#> P value adjustment: tukey method for comparing a family of 4 estimates
By default, the p-values and confidence intervals from the output and
plot above are adjusted according to the Tukey method. The model
suggests no significant evidence (p-values > 0.1) that the average
log conductivity is different among the states. Other p-value adjustment
methods can be passed via adjust. For example, we can use
the Bonferroni method instead of Tukey method
pairs(em11, adjust = "bonferroni")#> contrast estimate SE df t.ratio p.value
#> AZ - CO 0.9375 0.602 14.1 1.558 0.8493
#> AZ - NV -0.0852 0.576 17.8 -0.148 1.0000
#> AZ - UT 0.0263 0.548 12.5 0.048 1.0000
#> CO - NV -1.0226 0.589 15.2 -1.736 0.6172
#> CO - UT -0.9111 0.562 11.1 -1.621 0.7981
#> NV - UT 0.1115 0.530 15.1 0.210 1.0000
#>
#> Degrees-of-freedom method: satterthwaite
#> P value adjustment: bonferroni method for 6 testsor apply no adjustment method at all:
pairs(em11, adjust = "none")#> contrast estimate SE df t.ratio p.value
#> AZ - CO 0.9375 0.602 14.1 1.558 0.1416
#> AZ - NV -0.0852 0.576 17.8 -0.148 0.8840
#> AZ - UT 0.0263 0.548 12.5 0.048 0.9625
#> CO - NV -1.0226 0.589 15.2 -1.736 0.1029
#> CO - UT -0.9111 0.562 11.1 -1.621 0.1330
#> NV - UT 0.1115 0.530 15.1 0.210 0.8362
#>
#> Degrees-of-freedom method: satterthwaiteA Multi-Factor Model
Now we explore a model that adds a second factor: year,
with two levels, 2012 and 2017:
spmod2 <- splm(
formula = log_cond ~ state + year,
data = lake,
spcov_type = "exponential"
)We can view the factors separately:
em21 <- emmeans(spmod2, ~ state)
em21#> state emmean SE df lower.CL upper.CL
#> AZ 5.77 0.415 16.8 4.90 6.65
#> CO 4.86 0.431 12.7 3.93 5.79
#> NV 5.87 0.396 20.6 5.04 6.69
#> UT 5.69 0.366 10.3 4.88 6.51
#>
#> Results are averaged over the levels of: year
#> Degrees-of-freedom method: satterthwaite
#> Confidence level used: 0.95
em22 <- emmeans(spmod2, ~ year)
em22#> year emmean SE df lower.CL upper.CL
#> 2012 5.67 0.208 15.7 5.23 6.11
#> 2017 5.43 0.261 33.7 4.90 5.96
#>
#> Results are averaged over the levels of: state
#> Degrees-of-freedom method: satterthwaite
#> Confidence level used: 0.95Or, we can view the factors simultaneously (by providing both
variables separated by +):
em23 <- emmeans(spmod2, ~ state + year)
em23#> state year emmean SE df lower.CL upper.CL
#> AZ 2012 5.89 0.417 17.36 5.02 6.77
#> CO 2012 4.98 0.438 13.50 4.04 5.92
#> NV 2012 5.99 0.401 20.95 5.15 6.82
#> UT 2012 5.81 0.353 9.17 5.02 6.61
#> AZ 2017 5.66 0.443 21.16 4.73 6.58
#> CO 2017 4.74 0.453 15.55 3.78 5.70
#> NV 2017 5.75 0.423 27.08 4.88 6.62
#> UT 2017 5.57 0.412 15.75 4.70 6.45
#>
#> Degrees-of-freedom method: satterthwaite
#> Confidence level used: 0.95
plot(em23)
A Multi-Factor Model With an Interaction
We can supplement the model with an interaction, which lets the
effect of state to vary by year. Recall that shorthand for
state + year + state:year is state * year:
spmod3 <- splm(
formula = log_cond ~ state * year,
data = lake,
spcov_type = "exponential"
)Because the effect of state varies by year, single-variable summaries
of emmeans can be misleading, which emmeans
warns users about:
emmeans(spmod3, ~ state)#> NOTE: Results may be misleading due to involvement in interactions#> state emmean SE df lower.CL upper.CL
#> AZ 5.70 0.435 17.3 4.79 6.62
#> CO 4.88 0.441 12.5 3.93 5.84
#> NV 5.82 0.413 22.6 4.97 6.68
#> UT 5.76 0.398 12.8 4.90 6.62
#>
#> Results are averaged over the levels of: year
#> Degrees-of-freedom method: satterthwaite
#> Confidence level used: 0.95Instead, it is helpful to quantify the effect of state separately for each year:
em31 <- emmeans(spmod2, ~ state, by = "year")
em31#> year = 2012:
#> state emmean SE df lower.CL upper.CL
#> AZ 5.89 0.417 17.36 5.02 6.77
#> CO 4.98 0.438 13.50 4.04 5.92
#> NV 5.99 0.401 20.95 5.15 6.82
#> UT 5.81 0.353 9.17 5.02 6.61
#>
#> year = 2017:
#> state emmean SE df lower.CL upper.CL
#> AZ 5.66 0.443 21.16 4.73 6.58
#> CO 4.74 0.453 15.55 3.78 5.70
#> NV 5.75 0.423 27.08 4.88 6.62
#> UT 5.57 0.412 15.75 4.70 6.45
#>
#> Degrees-of-freedom method: satterthwaite
#> Confidence level used: 0.95
pairs(em31)#> year = 2012:
#> contrast estimate SE df t.ratio p.value
#> AZ - CO 0.9168 0.596 14.3 1.539 0.4416
#> AZ - NV -0.0944 0.570 18.1 -0.165 0.9983
#> AZ - UT 0.0826 0.545 12.8 0.152 0.9987
#> CO - NV -1.0112 0.583 15.4 -1.734 0.3406
#> CO - UT -0.8342 0.560 11.4 -1.490 0.4738
#> NV - UT 0.1770 0.528 15.3 0.335 0.9865
#>
#> year = 2017:
#> contrast estimate SE df t.ratio p.value
#> AZ - CO 0.9168 0.596 14.3 1.539 0.4416
#> AZ - NV -0.0944 0.570 18.1 -0.165 0.9983
#> AZ - UT 0.0826 0.545 12.8 0.152 0.9987
#> CO - NV -1.0112 0.583 15.4 -1.734 0.3406
#> CO - UT -0.8342 0.560 11.4 -1.490 0.4738
#> NV - UT 0.1770 0.528 15.3 0.335 0.9865
#>
#> Degrees-of-freedom method: satterthwaite
#> P value adjustment: tukey method for comparing a family of 4 estimates
plot(em31)
And similarly, we can quantify the effect of year separately for each state
#> state = AZ:
#> contrast estimate SE df t.ratio p.value
#> year2012 - year2017 0.239 0.227 63.6 1.053 0.2965
#>
#> state = CO:
#> contrast estimate SE df t.ratio p.value
#> year2012 - year2017 0.239 0.227 63.6 1.053 0.2965
#>
#> state = NV:
#> contrast estimate SE df t.ratio p.value
#> year2012 - year2017 0.239 0.227 63.6 1.053 0.2965
#>
#> state = UT:
#> contrast estimate SE df t.ratio p.value
#> year2012 - year2017 0.239 0.227 63.6 1.053 0.2965
#>
#> Degrees-of-freedom method: satterthwaiteAnd we can quantify the effect of each combination of
state and year:
em33 <- emmeans(spmod2, ~ state + year)
em33#> state year emmean SE df lower.CL upper.CL
#> AZ 2012 5.89 0.417 17.36 5.02 6.77
#> CO 2012 4.98 0.438 13.50 4.04 5.92
#> NV 2012 5.99 0.401 20.95 5.15 6.82
#> UT 2012 5.81 0.353 9.17 5.02 6.61
#> AZ 2017 5.66 0.443 21.16 4.73 6.58
#> CO 2017 4.74 0.453 15.55 3.78 5.70
#> NV 2017 5.75 0.423 27.08 4.88 6.62
#> UT 2017 5.57 0.412 15.75 4.70 6.45
#>
#> Degrees-of-freedom method: satterthwaite
#> Confidence level used: 0.95
pairs(em33)#> contrast estimate SE df t.ratio p.value
#> AZ year2012 - CO year2012 0.9168 0.596 14.3 1.539 0.7761
#> AZ year2012 - NV year2012 -0.0944 0.570 18.1 -0.165 1.0000
#> AZ year2012 - UT year2012 0.0826 0.545 12.8 0.152 1.0000
#> AZ year2012 - AZ year2017 0.2389 0.227 63.6 1.053 0.9640
#> AZ year2012 - CO year2017 1.1557 0.630 18.1 1.833 0.6077
#> AZ year2012 - NV year2017 0.1445 0.610 24.4 0.237 1.0000
#> AZ year2012 - UT year2017 0.3215 0.609 19.4 0.528 0.9993
#> CO year2012 - NV year2012 -1.0112 0.583 15.4 -1.734 0.6680
#> CO year2012 - UT year2012 -0.8342 0.560 11.4 -1.490 0.7986
#> CO year2012 - AZ year2017 -0.6779 0.645 18.9 -1.052 0.9594
#> CO year2012 - CO year2017 0.2389 0.227 63.6 1.053 0.9640
#> CO year2012 - NV year2017 -0.7723 0.630 21.1 -1.226 0.9148
#> CO year2012 - UT year2017 -0.5953 0.630 17.4 -0.945 0.9767
#> NV year2012 - UT year2012 0.1770 0.528 15.3 0.335 1.0000
#> NV year2012 - AZ year2017 0.3332 0.617 23.2 0.540 0.9993
#> NV year2012 - CO year2017 1.2501 0.622 19.4 2.010 0.5004
#> NV year2012 - NV year2017 0.2389 0.227 63.6 1.053 0.9640
#> NV year2012 - UT year2017 0.4159 0.598 22.8 0.695 0.9963
#> UT year2012 - AZ year2017 0.1562 0.570 15.6 0.274 1.0000
#> UT year2012 - CO year2017 1.0731 0.577 13.3 1.861 0.5951
#> UT year2012 - NV year2017 0.0619 0.550 19.3 0.112 1.0000
#> UT year2012 - UT year2017 0.2389 0.227 63.6 1.053 0.9640
#> AZ year2017 - CO year2017 0.9168 0.596 14.3 1.539 0.7761
#> AZ year2017 - NV year2017 -0.0944 0.570 18.1 -0.165 1.0000
#> AZ year2017 - UT year2017 0.0826 0.545 12.8 0.152 1.0000
#> CO year2017 - NV year2017 -1.0112 0.583 15.4 -1.734 0.6680
#> CO year2017 - UT year2017 -0.8342 0.560 11.4 -1.490 0.7986
#> NV year2017 - UT year2017 0.1770 0.528 15.3 0.335 1.0000
#>
#> Degrees-of-freedom method: satterthwaite
#> P value adjustment: tukey method for comparing a family of 8 estimatesA Single-Factor Numeric Model With an Interaction
Suppose it is of interest to supplement the state model
(spmod1) with a continuous temperature variable:
spmod4 <- splm(
formula = log_cond ~ state * temp,
data = lake,
spcov_type = "exponential"
)Because our model has a state-by-year interaction, the effect of
temperature varies by state. Supplying the by argument lets
us quantify the effect of state at the average temperature value:
em41 <- emmeans(spmod4, ~ state, by = "temp")
em41#> temp = 7.63:
#> state emmean SE df lower.CL upper.CL
#> AZ 4.67 0.304 58.3 4.06 5.27
#> CO 5.70 0.209 14.2 5.25 6.15
#> NV 5.64 0.203 62.6 5.24 6.05
#> UT 6.05 0.143 39.5 5.76 6.34
#>
#> Degrees-of-freedom method: satterthwaite
#> Confidence level used: 0.95
pairs(em41)#> temp = 7.63:
#> contrast estimate SE df t.ratio p.value
#> AZ - CO -1.0345 0.368 35.2 -2.808 0.0386
#> AZ - NV -0.9771 0.365 59.6 -2.677 0.0460
#> AZ - UT -1.3829 0.336 54.4 -4.122 0.0007
#> CO - NV 0.0573 0.291 27.3 0.197 0.9972
#> CO - UT -0.3484 0.253 18.9 -1.377 0.5280
#> NV - UT -0.4058 0.248 54.0 -1.637 0.3672
#>
#> Degrees-of-freedom method: satterthwaite
#> P value adjustment: tukey method for comparing a family of 4 estimatesIf we want to quantify the effect of state at specific temperature
values, we can supply them via the at argument:
#> temp = 2:
#> contrast estimate SE df t.ratio p.value
#> AZ - CO -0.3803 0.505 56.9 -0.753 0.8751
#> AZ - NV -0.8775 0.590 61.7 -1.487 0.4513
#> AZ - UT -0.4471 0.496 50.9 -0.901 0.8042
#> CO - NV -0.4972 0.435 63.5 -1.144 0.6638
#> CO - UT -0.0668 0.295 38.1 -0.227 0.9958
#> NV - UT 0.4304 0.424 55.1 1.016 0.7411
#>
#> temp = 8:
#> contrast estimate SE df t.ratio p.value
#> AZ - CO -1.0769 0.366 33.1 -2.940 0.0289
#> AZ - NV -0.9836 0.356 60.1 -2.762 0.0371
#> AZ - UT -1.4436 0.330 55.0 -4.371 0.0003
#> CO - NV 0.0933 0.295 25.9 0.316 0.9888
#> CO - UT -0.3667 0.263 19.0 -1.392 0.5193
#> NV - UT -0.4600 0.249 56.1 -1.847 0.2627
#>
#> Degrees-of-freedom method: satterthwaite
#> P value adjustment: tukey method for comparing a family of 4 estimatesWe use emmip() to visualize the change in the effect of
state at varying temperature values:
And emtrends to quantify the effect of temperature
separately for each state:
em43 <- emtrends(spmod4, ~ state, var = "temp")
em43#> state temp.trend SE df lower.CL upper.CL
#> AZ 0.159 0.0317 56.7 0.0951 0.222
#> CO 0.275 0.0424 24.4 0.1872 0.362
#> NV 0.176 0.0491 81.5 0.0786 0.274
#> UT 0.325 0.0365 55.2 0.2515 0.398
#>
#> Degrees-of-freedom method: satterthwaite
#> Confidence level used: 0.95A Multi-Factor Numeric Model With Interactions
We can extend the single-factor numeric model to include multiple factors:
spmod5 <- splm(log_cond ~ state * year * temp, data = lake, spcov_type = "exponential")
em51 <- emmeans(spmod5, ~ state, by = c("temp", "year"))
em51#> temp = 7.63, year = 2012:
#> state emmean SE df lower.CL upper.CL
#> AZ 4.93 0.341 59.0 4.24 5.61
#> CO 5.71 0.246 14.9 5.18 6.23
#> NV 5.76 0.250 74.7 5.26 6.26
#> UT 6.06 0.145 38.1 5.76 6.35
#>
#> temp = 7.63, year = 2017:
#> state emmean SE df lower.CL upper.CL
#> AZ 3.71 0.647 85.6 2.43 5.00
#> CO 5.69 0.314 50.4 5.06 6.32
#> NV 5.39 0.352 83.2 4.69 6.09
#> UT 6.17 1.445 56.8 3.27 9.06
#>
#> Degrees-of-freedom method: satterthwaite
#> Confidence level used: 0.95#> temp = 7.63, state = AZ:
#> year emmean SE df lower.CL upper.CL
#> 2012 4.93 0.341 59.0 4.24 5.61
#> 2017 3.71 0.647 85.6 2.43 5.00
#>
#> temp = 7.63, state = CO:
#> year emmean SE df lower.CL upper.CL
#> 2012 5.71 0.246 14.9 5.18 6.23
#> 2017 5.69 0.314 50.4 5.06 6.32
#>
#> temp = 7.63, state = NV:
#> year emmean SE df lower.CL upper.CL
#> 2012 5.76 0.250 74.7 5.26 6.26
#> 2017 5.39 0.352 83.2 4.69 6.09
#>
#> temp = 7.63, state = UT:
#> year emmean SE df lower.CL upper.CL
#> 2012 6.06 0.145 38.1 5.76 6.35
#> 2017 6.17 1.445 56.8 3.27 9.06
#>
#> Degrees-of-freedom method: satterthwaite
#> Confidence level used: 0.95#> temp = 7.63:
#> contrast estimate SE df t.ratio p.value
#> AZ year2012 - CO year2012 -0.7827 0.420 35.5 -1.862 0.5840
#> AZ year2012 - NV year2012 -0.8365 0.423 64.7 -1.979 0.5037
#> AZ year2012 - UT year2012 -1.1311 0.370 55.4 -3.055 0.0635
#> AZ year2012 - AZ year2017 1.2113 0.715 84.3 1.694 0.6910
#> AZ year2012 - CO year2017 -0.7609 0.463 55.0 -1.644 0.7222
#> AZ year2012 - NV year2017 -0.4668 0.490 73.9 -0.953 0.9795
#> AZ year2012 - UT year2017 -1.2399 1.484 60.7 -0.835 0.9903
#> CO year2012 - NV year2012 -0.0537 0.351 33.0 -0.153 1.0000
#> CO year2012 - UT year2012 -0.3484 0.286 18.5 -1.218 0.9162
#> CO year2012 - AZ year2017 1.9940 0.692 81.1 2.881 0.0902
#> CO year2012 - CO year2017 0.0219 0.368 83.1 0.059 1.0000
#> CO year2012 - NV year2017 0.3160 0.430 51.9 0.735 0.9955
#> CO year2012 - UT year2017 -0.4572 1.466 61.1 -0.312 1.0000
#> NV year2012 - UT year2012 -0.2946 0.289 64.8 -1.018 0.9701
#> NV year2012 - AZ year2017 2.0477 0.694 85.7 2.952 0.0749
#> NV year2012 - CO year2017 0.0756 0.401 59.8 0.188 1.0000
#> NV year2012 - NV year2017 0.3697 0.419 79.8 0.882 0.9869
#> NV year2012 - UT year2017 -0.4034 1.466 58.5 -0.275 1.0000
#> UT year2012 - AZ year2017 2.3424 0.663 85.7 3.533 0.0147
#> UT year2012 - CO year2017 0.3702 0.346 48.0 1.071 0.9597
#> UT year2012 - NV year2017 0.6644 0.381 79.0 1.743 0.6593
#> UT year2012 - UT year2017 -0.1088 1.447 55.6 -0.075 1.0000
#> AZ year2017 - CO year2017 -1.9721 0.719 84.5 -2.743 0.1242
#> AZ year2017 - NV year2017 -1.6780 0.737 85.6 -2.278 0.3178
#> AZ year2017 - UT year2017 -2.4512 1.583 63.2 -1.549 0.7781
#> CO year2017 - NV year2017 0.2941 0.472 71.4 0.624 0.9984
#> CO year2017 - UT year2017 -0.4790 1.478 60.5 -0.324 1.0000
#> NV year2017 - UT year2017 -0.7731 1.487 59.5 -0.520 0.9995
#>
#> Degrees-of-freedom method: satterthwaite
#> P value adjustment: tukey method for comparing a family of 8 estimates
em54 <- emtrends(spmod5, ~ state, by = "year", var = "temp")
em54#> year = 2012:
#> state temp.trend SE df lower.CL upper.CL
#> AZ 0.145 0.0359 67.3 0.0733 0.217
#> CO 0.269 0.0487 27.7 0.1694 0.369
#> NV 0.168 0.1030 83.4 -0.0370 0.373
#> UT 0.324 0.0391 70.1 0.2458 0.402
#>
#> year = 2017:
#> state temp.trend SE df lower.CL upper.CL
#> AZ 0.225 0.0657 85.3 0.0943 0.355
#> CO 0.292 0.0707 64.8 0.1512 0.434
#> NV 0.178 0.0563 80.6 0.0664 0.290
#> UT 0.359 0.2426 59.3 -0.1268 0.844
#>
#> Degrees-of-freedom method: satterthwaite
#> Confidence level used: 0.95
em55 <- emtrends(spmod5, ~ year, by = "state", var = "temp")
em55#> state = AZ:
#> year temp.trend SE df lower.CL upper.CL
#> 2012 0.145 0.0359 67.3 0.0733 0.217
#> 2017 0.225 0.0657 85.3 0.0943 0.355
#>
#> state = CO:
#> year temp.trend SE df lower.CL upper.CL
#> 2012 0.269 0.0487 27.7 0.1694 0.369
#> 2017 0.292 0.0707 64.8 0.1512 0.434
#>
#> state = NV:
#> year temp.trend SE df lower.CL upper.CL
#> 2012 0.168 0.1030 83.4 -0.0370 0.373
#> 2017 0.178 0.0563 80.6 0.0664 0.290
#>
#> state = UT:
#> year temp.trend SE df lower.CL upper.CL
#> 2012 0.324 0.0391 70.1 0.2458 0.402
#> 2017 0.359 0.2426 59.3 -0.1268 0.844
#>
#> Degrees-of-freedom method: satterthwaite
#> Confidence level used: 0.95
em56 <- emtrends(spmod5, ~ state + year, var = "temp")
em56#> state year temp.trend SE df lower.CL upper.CL
#> AZ 2012 0.145 0.0359 67.3 0.0733 0.217
#> CO 2012 0.269 0.0487 27.7 0.1694 0.369
#> NV 2012 0.168 0.1030 83.4 -0.0370 0.373
#> UT 2012 0.324 0.0391 70.1 0.2458 0.402
#> AZ 2017 0.225 0.0657 85.3 0.0943 0.355
#> CO 2017 0.292 0.0707 64.8 0.1512 0.434
#> NV 2017 0.178 0.0563 80.6 0.0664 0.290
#> UT 2017 0.359 0.2426 59.3 -0.1268 0.844
#>
#> Degrees-of-freedom method: satterthwaite
#> Confidence level used: 0.95Pairing with anova() from spmodel
The anova() function from spmodel is
especially helpful to further contextualize emmeans output.
The emmeans functions are very helpful for contrasting
factor levels, but anova() is built to answer the question
“Are any of these factor levels significantly related to the response
variable?”. Recall spmod5, which quantifies the effects of
state, year, and temp (and their
interactions) on log conductivity:
summary(spmod5)#>
#> Call:
#> splm(formula = log_cond ~ state * year * temp, data = lake, spcov_type = "exponential")
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -2.354623 -0.404763 -0.001426 0.444013 2.801426
#>
#> Coefficients (fixed):
#> Estimate Std. Error df t value Pr(>|t|)
#> (Intercept) 3.81882 0.57406 61.09623 6.652 9.08e-09 ***
#> stateCO -0.16530 0.65566 61.32568 -0.252 0.801794
#> stateNV 0.66205 1.06415 78.21506 0.622 0.535657
#> stateUT -0.23392 0.63944 58.59585 -0.366 0.715810
#> year2017 -1.82141 1.21707 83.96986 -1.497 0.138258
#> temp 0.14498 0.03592 67.29314 4.036 0.000142 ***
#> stateCO:year2017 1.62235 1.33773 83.22223 1.213 0.228652
#> stateNV:year2017 1.37143 1.60001 82.40432 0.857 0.393853
#> stateUT:year2017 1.66483 1.34976 78.42032 1.233 0.221102
#> stateCO:temp 0.12418 0.06051 37.92851 2.052 0.047110 *
#> stateNV:temp 0.02285 0.10906 82.58084 0.209 0.834582
#> stateUT:temp 0.17879 0.05308 68.82831 3.368 0.001242 **
#> year2017:temp 0.07992 0.07309 81.99717 1.093 0.277391
#> stateCO:year2017:temp -0.05671 0.10825 84.18029 -0.524 0.601721
#> stateNV:year2017:temp -0.06941 0.13715 84.60862 -0.506 0.614127
#> stateUT:year2017:temp -0.04516 0.25289 53.38844 -0.179 0.858949
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Pseudo R-squared: 0.6694
#>
#> Coefficients (exponential spatial covariance):
#> de ie range
#> 6.197e-01 6.197e-05 9.249e+03We perform an analysis of variance by running:
anova(spmod5)#> Analysis of Variance Table
#>
#> Response: log_cond
#> NumDF DenDF F value Pr(>F)
#> (Intercept) 1 61.096 44.2537 9.082e-09 ***
#> state 3 65.978 0.3261 0.8064792
#> year 1 83.970 2.2397 0.1382576
#> temp 1 67.293 16.2874 0.0001417 ***
#> state:year 3 70.710 0.5515 0.6487916
#> state:temp 3 57.614 4.1171 0.0102788 *
#> year:temp 1 81.997 1.1957 0.2773906
#> state:year:temp 3 71.574 0.1308 0.9414474
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The analysis of variance suggests a significant intercept, temperature effect, and state-by-temperature interaction (p-values < 0.01) but no other significant effects (p-values > 0.1).
Uncertainty Quantification
Standard errors and confidence intervals returned by
emmeans for spmodel objects use the fixed
effect covariance matrix (returned by vcov()), accounting
for spatial dependence. For splm() and
spautor() model objects fit with
ddf = "satterthwaite" (the default when the sample size is
500 or fewer), emmeans also uses Satterthwaite denominator
degrees of freedom. Otherwise, asymptotic (Inf) degrees of
freedom are used. For spglm() and spgautor(),
asymptotic degrees of freedom are always used. See the Technical Details
vignette for more.
Additional Tools
We only showed a small subset of all possible tools that
emmeans can apply to spmodel objects.
Additional tools provide support for enhanced visualizations and
printing, contrasts, joint hypothesis testing, p-value adjustments, and
variable transformations (e.g., when using a spglm() or
spgautor() model object), among many others. To learn more
about emmeans, visit the website here.