How you quantify your model effects varies with your model structure (e.g. error distribution assumption).
If it is your first time reading this, read through the examples that present models assuming a normal error distribution assumption. This will help you understand why we are reporting modelled effects in this way. Then, (if relevant) look at the section on the error distribution assumption you are interested in for your model.
In this section, you will learn how to report your modelled effects directly - the change in your response given a unit change in your predictor. These effects are quantified in the estimates of your coefficients (i.e. the magnitude, direction and uncertainty of your coefficient estimates).
How you report your modelled effects will depend on:
whether you are reporting effects for a categorical or numeric predictor, and
your error distribution assumption
We will consider each of these now:
Effects are the change in your response due to a unit change in your predictor.
For categorical predictors, the effect is the change in the response when you change the level (category) of your categorical predictor.
You capture the effect of your categorical predictor on your response by reporting the modelled value of the response at each level (category) of your predictor.1
You also report the uncertainty around this modelled estimate. Two common uncertainty estimates are
standard errors (SE) which are a measure of uncertainty of the average modelled effect, and
the 95% confidence intervals around your modelled effects, which tell you the range your modelled effects would be expected to fall into (with a 95% probability) if you redid your experiment.
We will focus on the 95% confidence intervals as they tend to be easier to interpret for many different types of models.
When you have a categorical predictor with more than two categories (levels), you will also want to test to see where the effects on the response differ among the categories. This is done using multiple comparison testing, i.e. you will compare the modelled effect on your response of each level of your categorical predictor vs. every other level of your categorical predictor to determine which effects are different from each other (called pairwise testing).
This is done by testing the null hypothesis that the modelled effects of each pair of category levels are similar to one another, typically rejecting the null hypothesis when P < 0.05. Remember that the P-value is the probability that you observe a value at least as big as the one you observed even though your null hypothesis is true. In this case, you are looking at the value you are observing as the difference between coefficients estimated for the two levels of your predictor. The P-value is the probability of getting a difference at least as big as the one you observed even though there is actually no difference between the effects (i.e. the null hypothesis is true).
A few of things to note about multiple comparison testing:
Multiple comparison testing is only needed for categorical predictors.
Multiple comparison testing is only needed for categorical predictors with more than two levels (or categories).
Multiple comparison testing should only be done after your hypothesis testing has given you evidence that the predictor has an effect on your response. That is, you should only do a multiple comparison test on a predictor if that predictor is in your best-specified model. As this is a test done after your hypothesis testing, it is called a post-hoc2 test.
Multiple comparison testing can be a problem because you are essentially repeating a hypothesis test many times on the same data (i.e. are the effects of Sp1 different than those of Sp2? are the effects of Sp1 different than those of Sp3? are the effects of Sp2 different than those of Sp3?…). These repeated tests mean there is a high chance that you will find a difference in one test purely due to random chance, vs. due to there being an actual difference. The multiple comparison tests you will perform have been formulated to correct for this increased error.
In summary: reporting effects of categorical predictors on your response means reporting:
the modelled value of the response at each predictor level,
the uncertainty of this modelled value (95% confidence interval), and
the evidence that the effects among levels are different from one another.
Effects are the change in your response due to a unit change in your predictor.
For numeric predictors, the effect is the change in the response when you increase your numeric predictor by one unit.
You capture the effect of your numeric predictor on your response by reporting this rate of change as a slope. When the slope (coefficient) is positive, you have an increase in your response with a unit increase in your predictor. When the slope (coefficient) is negative, you have a decrease in your response with a unit increase in your predictor. When the slope (coefficient) is zero, you have no change in your response for a unit increase in your predictor, and you would conclude that there is no evidence that your predictor explains variability in your response.
Your best-specified model might also include an interaction term. For example, you might have an interaction between categorical and numeric predictors. To report the effect of an interaction, you report the effect of the numeric predictor for each level in your categorical predictor. See Example 4 below for how to do this.
When you use a GLM to test your hypothesis, the coefficients are estimated in the “linked world”. You need to convert these back to the response (real) world in order to communicate your effects.
This is because of the difference between the “link” scale and “response” scale for generalized linear models.
The link scale is where the model fitting happens, and where the evidence for your modelled effects is gathered.
The response scale is back in the real world - these are the actual effects in your response you would expect to see with a change in your predictor.
How you do the conversion to interpret your coefficients depends on your specific error distribution assumption as you need to consider the link function in your model.
For a normal error distribution assumption, the effect of a numeric predictor on the response is the absolute change in the response for a unit change in the predictor.
If you have a normal error distribution assumption (and are using link = "identity"), your coefficient on the response scale is identical to your coefficient on the link scale.
Here’s a model equation for a generalized linear model:
\[ \begin{align} g(\mu_i|Pred_i) &= \beta_1 \cdot Pred_i + \beta_0 \\ Resp_i &\sim F(\mu_i) \end{align} \]
where * \(g()\) is the link function
\(\mu_i\) is the mean fitted value \(\mu_i\) dependent on \(Pred_i\)
\(\beta_1\) is the coefficient of \(Pred_i\)
\(Pred_i\) is the value of your predictor
\(\beta_0\) is the intercept
\(Resp_i\) is the value of the response variable
\(F\) represents some error distribution assumption
As above, the effect of a numeric predictor on your response is the change in your fitted value (\(g(\mu_i)\)) for a unit change in your predictor (\(Pred_i\)):
\[
\begin{align}
g(\mu|Pred+1) - g(\mu|Pred) &= (\beta_1 \cdot (Pred+1) + \beta_0) - (\beta_1 \cdot Pred_i + \beta_0) \\
&=\beta_1 \cdot Pred + \beta_1 + \beta_0 - \beta_1 \cdot Pred_i - \beta_0 \\
&=\beta_1
\end{align}
\] For linear models with a normal error distribution assumption (i.e. when you are using link = "identity"), the link function \(g()\) is not there. For this reason, you can read the effects of your model directly from the model coefficients (\(\beta_1\), \(\beta_0\)) when you have a normal error distribution assumption. And for this reason, you will need to convert your model coefficients to quantify your modelled effects when you have an error distribution assumption that is not normal/uses a link function that is not “identity” (see next section).
Recall that Example 1 contains only one predictor and it is categorical:
Resp1 ~ Cat1 + 1
What are your modelled effects (with uncertainty)?
In the chapter on your Starting Model, you found your coefficients in the “Estimate” column of the summary() output of your model:
coef(summary(bestMod1)) # extract the coefficients from summary() Estimate Std. Error t value Pr(>|t|)
(Intercept) 20.954839 1.316231 15.920341 8.313153e-29
Cat1Sp2 5.876740 1.773637 3.313384 1.296521e-03
Cat1Sp3 -4.509677 1.861431 -2.422694 1.726174e-02Interpreting the coefficients from this output takes practice - especially for categorical predictors because of the way that R treats categorical predictors in regression. With R’s “dummy coding”, one level of the predictor (here Sp1) is incorporated into the intercept estimate (21) as the reference level. The other coefficients in the Estimate column represent the change in modelled response when you move from the reference level (Sp1) to another level.
The modelled response when Cat1 = Sp2 is 5.9 units higher than this reference level, or 21 + 5.9 = 26.9 units.
The modelled Resp1 when Cat1 = Sp3 is -4.5 units lower than the reference level, or 21 -4.5 = 16.5 units.3
So all the information we need is in this summary() output, but not easy to see immediately. An easier way is to use the emmeans4 package which helps you by reporting the coefficients directly5. In our case, you use the emmeans package to get the mean value of the response at each level of the predictor. For categorical predictors, you do this with the emmeans() function:
library(emmeans) # load the emmeans package
emmOut <- emmeans(object = bestMod1, # your model
specs = ~ Cat1, # your predictors
type = "response") # report coefficients on the response scale
emmOut Cat1 emmean SE df lower.CL upper.CL
Sp1 21.0 1.32 97 18.3 23.6
Sp2 26.8 1.19 97 24.5 29.2
Sp3 16.4 1.32 97 13.8 19.1
Confidence level used: 0.95 So now you have a modelled value of your response for each level of your categorical predictor - this captures the effect of your categorical predictor on your response. You also need to report uncertainty around this effect, and you can do this by reporting the 95% confidence level (CL) also reported by the emmeans() function.
When Cat1 is Sp1, Resp1 is estimated to be 21 units (95% confidence level: 18.3 to 23.6 units).
When Cat1 is Sp2, Resp1 is estimated to be 26.8 units (95% confidence level: 24.5 to 29.2 units).
When Cat1 is Sp3, Resp1 is estimated to be 16.4 units (95% confidence level: 13.8 to 19.1 units.
Note that this is the same information in the summary() output just easier to read.6
Note also that two types of uncertainty are measured here. SE stands for the standard error around the prediction, and is a measure of uncertainty of the average modelled effect. The lower.CL and upper.CL represent the 95% confidence limits of the prediction - so if I observed a new Resp1 at a particular Cat1, there would be a 95% chance it would lie between the bounds of the confidence limits.
Finally, note that you can also get a quick plot of the effects by handing the emmeans() output to plot().
Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)
Multiple comparison testing is quick to do with the emmeans package. It just requires you to hand the output from emmeans() to a function called pairs():
forComp <- pairs(emmOut, # your emmeans object
adjust = 'fdr') # multiple comparison adjustment
forComp contrast estimate SE df t.ratio p.value
Sp1 - Sp2 -5.88 1.77 97 -3.313 0.0019
Sp1 - Sp3 4.51 1.86 97 2.423 0.0173
Sp2 - Sp3 10.39 1.77 97 5.856 <0.0001
P value adjustment: fdr method for 3 tests The output shows the results of the multiple comparison (pairwise) testing. The values in the p.value column tell you the results of the hypothesis test comparing the coefficients between the two levels of your categorical predictor. For example, for Sp1 vs. Sp2, there is a 0.19% (p = 0.0019) probability of getting a difference in coefficients at least as big as 26.8 - 21 = 5.8 even though the null hypothesis (no difference) is true. This value 0.19% (P = 0.0019) is small enough that you can say you have evidence that the effect of Cat1 on Resp1 is different for these two different levels (Sp1 and Sp2).
Based on a threshold p-value of 0.05, you can see that:
There is evidence that the value of Resp1 when Cat1 is Sp1 is different (lower) than that when Cat1 is Sp2 as P = 0.0019 is less than P = 0.05.
There is evidence that the value of Resp1 when Cat1 is Sp1 is different (higher) than that when Cat1 is Sp3 as P = 0.0173 is less than P = 0.05.
There is evidence that the value of Resp1 when Cat1 is Sp2 is different (higher) than that when Cat1 is Sp3 as P < 0.00017 is smaller than P = 0.05.8.
Note that the p-values have been adjusted via the “false discovery rate” (fdr) method (Verhoeven, Simonsen, and McIntyre 2005) which adjusts the difference that the two coefficients need to have to allow for the fact that we are making multiple comparisons.9
Note that you can also get the results from the pairwise testing visually by handing the output of pairs() to plot().
Recall that Example 2 contains two predictors and both are categorical:
Resp2 ~ Cat2 + Cat3 + Cat2:Cat3 + 1
What are your modelled effects (with uncertainty)?
Again, for categorical predictors, there is a coefficient estimated for each level of the predictor. If there is one or more interactions among predictors, there will be one coefficient for each combination of levels among predictors. Let’s look at the summary() output of your model to understand better:
coef(summary(bestMod2)) # extract the coefficients from summary() Estimate Std. Error t value Pr(>|t|)
(Intercept) 364.90879 24.73850 14.7506420 1.629062e-25
Cat2TypeB 37.37325 32.98467 1.1330491 2.602712e-01
Cat2TypeC -75.94350 33.87459 -2.2419017 2.748208e-02
Cat2TypeD -141.96558 38.32472 -3.7042819 3.695537e-04
Cat3Treat2 -43.62974 36.41409 -1.1981554 2.340735e-01
Cat3Control 64.11573 30.29835 2.1161457 3.715672e-02
Cat2TypeB:Cat3Treat2 42.70879 53.60009 0.7968045 4.277091e-01
Cat2TypeC:Cat3Treat2 73.60837 54.15228 1.3592849 1.775295e-01
Cat2TypeD:Cat3Treat2 66.02110 55.13228 1.1975037 2.343260e-01
Cat2TypeB:Cat3Control 81.25030 43.24327 1.8789121 6.356678e-02
Cat2TypeC:Cat3Control -19.07730 41.29748 -0.4619483 6.452584e-01
Cat2TypeD:Cat3Control -71.72444 46.17117 -1.5534463 1.239057e-01Comparing this output to that of Example 1 shows many more estimated coefficients for Example 2. This is because we have one coefficient for each level of each predictor, as well as coefficients for each combination (the interaction term) of levels of the predictors.
Again, it takes a little practice to read the coefficients from the summary() output. For Example 2:
The modelled prediction for Resp2 when Cat2 is TypeA and Cat3 is Treat1 is 365 units (the intercept).
The modelled prediction for Resp2 when Cat2 is TypeB and Cat3 is Treat1 is 365 + 37 = 402 units.
The modelled prediction for Resp2 when Cat2 is TypeA and Cat3 is Treat2 is 365 - 44 = 321 units.
The modelled prediction for Resp2 when Cat2 is TypeC and Cat3 is Treat2 is 365 - 76 - 44 + 74 = 319 units.
As above, we can use the emmeans package to more easily see these coefficients:
emmOut <- emmeans(object = bestMod2, # your model
specs = ~ Cat2 + Cat3 + Cat2:Cat3, # your predictors
type = "response") # report coefficients on the response scale
emmOut Cat2 Cat3 emmean SE df lower.CL upper.CL
TypeA Treat1 365 24.7 88 316 414
TypeB Treat1 402 21.8 88 359 446
TypeC Treat1 289 23.1 88 243 335
TypeD Treat1 223 29.3 88 165 281
TypeA Treat2 321 26.7 88 268 374
TypeB Treat2 401 32.7 88 336 466
TypeC Treat2 319 32.7 88 254 384
TypeD Treat2 245 29.3 88 187 304
TypeA Control 429 17.5 88 394 464
TypeB Control 548 21.8 88 504 591
TypeC Control 334 15.9 88 302 366
TypeD Control 215 18.9 88 178 253
Confidence level used: 0.95 So now we have a modelled value of our response for each level of our categorical predictor. For example:
The modelled prediction for Resp2 when Cat2 is TypeA and Cat3 is Treat1 is 365 units (95% confidence limit: 316-414).
The modelled prediction for Resp2 when Cat2 is TypeB and Cat3 is Treat1 is 402 units (95% confidence limit: 359-446).
The modelled prediction for Resp2 when Cat2 is TypeA and Cat3 is Treat2 is 321 units (95% confidence limit: 268-374).
The modelled prediction for Resp2 when Cat2 is TypeC and Cat3 is Treat2 is 319 units (95% confidence limit: 254-384).
Finally, note that you can also get a quick plot of the effects by handing the emmeans() output to plot().
Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)
As with Example 1, you can also find out which combinations of predictor levels are leading to statistically different model predictions in Resp2:
forComp <- pairs(emmOut, # your emmeans object
adjust = 'fdr') # multiple comparison adjustment
forComp contrast estimate SE df t.ratio p.value
TypeA Treat1 - TypeB Treat1 -37.373 33.0 88 -1.133 0.3303
TypeA Treat1 - TypeC Treat1 75.944 33.9 88 2.242 0.0465
TypeA Treat1 - TypeD Treat1 141.966 38.3 88 3.704 0.0010
TypeA Treat1 - TypeA Treat2 43.630 36.4 88 1.198 0.3090
TypeA Treat1 - TypeB Treat2 -36.452 41.0 88 -0.889 0.4361
TypeA Treat1 - TypeC Treat2 45.965 41.0 88 1.120 0.3307
TypeA Treat1 - TypeD Treat2 119.574 38.3 88 3.120 0.0054
TypeA Treat1 - TypeA Control -64.116 30.3 88 -2.116 0.0584
TypeA Treat1 - TypeB Control -182.739 33.0 88 -5.540 <0.0001
TypeA Treat1 - TypeC Control 30.905 29.4 88 1.051 0.3617
TypeA Treat1 - TypeD Control 149.574 31.1 88 4.805 <0.0001
TypeB Treat1 - TypeC Treat1 113.317 31.8 88 3.563 0.0016
TypeB Treat1 - TypeD Treat1 179.339 36.5 88 4.912 <0.0001
TypeB Treat1 - TypeA Treat2 81.003 34.5 88 2.348 0.0367
TypeB Treat1 - TypeB Treat2 0.921 39.3 88 0.023 0.9814
TypeB Treat1 - TypeC Treat2 83.338 39.3 88 2.119 0.0584
TypeB Treat1 - TypeD Treat2 156.947 36.5 88 4.299 0.0002
TypeB Treat1 - TypeA Control -26.742 28.0 88 -0.956 0.4098
TypeB Treat1 - TypeB Control -145.366 30.9 88 -4.711 <0.0001
TypeB Treat1 - TypeC Control 68.278 27.0 88 2.531 0.0248
TypeB Treat1 - TypeD Control 186.947 28.9 88 6.477 <0.0001
TypeC Treat1 - TypeD Treat1 66.022 37.3 88 1.769 0.1128
TypeC Treat1 - TypeA Treat2 -32.314 35.3 88 -0.914 0.4280
TypeC Treat1 - TypeB Treat2 -112.396 40.1 88 -2.804 0.0128
TypeC Treat1 - TypeC Treat2 -29.979 40.1 88 -0.748 0.5038
TypeC Treat1 - TypeD Treat2 43.631 37.3 88 1.169 0.3176
TypeC Treat1 - TypeA Control -140.059 29.0 88 -4.828 <0.0001
TypeC Treat1 - TypeB Control -258.683 31.8 88 -8.134 <0.0001
TypeC Treat1 - TypeC Control -45.038 28.1 88 -1.605 0.1510
TypeC Treat1 - TypeD Control 73.631 29.9 88 2.465 0.0279
TypeD Treat1 - TypeA Treat2 -98.336 39.6 88 -2.481 0.0275
TypeD Treat1 - TypeB Treat2 -178.418 43.9 88 -4.064 0.0003
TypeD Treat1 - TypeC Treat2 -96.001 43.9 88 -2.186 0.0519
TypeD Treat1 - TypeD Treat2 -22.391 41.4 88 -0.541 0.6383
TypeD Treat1 - TypeA Control -206.081 34.1 88 -6.043 <0.0001
TypeD Treat1 - TypeB Control -324.705 36.5 88 -8.894 <0.0001
TypeD Treat1 - TypeC Control -111.061 33.3 88 -3.335 0.0029
TypeD Treat1 - TypeD Control 7.609 34.8 88 0.218 0.8535
TypeA Treat2 - TypeB Treat2 -80.082 42.2 88 -1.895 0.0920
TypeA Treat2 - TypeC Treat2 2.335 42.2 88 0.055 0.9708
TypeA Treat2 - TypeD Treat2 75.945 39.6 88 1.916 0.0899
TypeA Treat2 - TypeA Control -107.746 31.9 88 -3.374 0.0027
TypeA Treat2 - TypeB Control -226.369 34.5 88 -6.562 <0.0001
TypeA Treat2 - TypeC Control -12.725 31.1 88 -0.409 0.7158
TypeA Treat2 - TypeD Control 105.945 32.7 88 3.237 0.0039
TypeB Treat2 - TypeC Treat2 82.417 46.3 88 1.781 0.1125
TypeB Treat2 - TypeD Treat2 156.026 43.9 88 3.554 0.0016
TypeB Treat2 - TypeA Control -27.663 37.1 88 -0.745 0.5038
TypeB Treat2 - TypeB Control -146.287 39.3 88 -3.719 0.0010
TypeB Treat2 - TypeC Control 67.357 36.4 88 1.852 0.0989
TypeB Treat2 - TypeD Control 186.027 37.8 88 4.923 <0.0001
TypeC Treat2 - TypeD Treat2 73.609 43.9 88 1.677 0.1336
TypeC Treat2 - TypeA Control -110.081 37.1 88 -2.967 0.0083
TypeC Treat2 - TypeB Control -228.704 39.3 88 -5.815 <0.0001
TypeC Treat2 - TypeC Control -15.060 36.4 88 -0.414 0.7158
TypeC Treat2 - TypeD Control 103.609 37.8 88 2.742 0.0148
TypeD Treat2 - TypeA Control -183.690 34.1 88 -5.387 <0.0001
TypeD Treat2 - TypeB Control -302.313 36.5 88 -8.281 <0.0001
TypeD Treat2 - TypeC Control -88.669 33.3 88 -2.663 0.0179
TypeD Treat2 - TypeD Control 30.000 34.8 88 0.861 0.4455
TypeA Control - TypeB Control -118.624 28.0 88 -4.242 0.0002
TypeA Control - TypeC Control 95.021 23.6 88 4.023 0.0004
TypeA Control - TypeD Control 213.690 25.7 88 8.299 <0.0001
TypeB Control - TypeC Control 213.644 27.0 88 7.918 <0.0001
TypeB Control - TypeD Control 332.314 28.9 88 11.514 <0.0001
TypeC Control - TypeD Control 118.669 24.7 88 4.809 <0.0001
P value adjustment: fdr method for 66 tests This is a hard table to navigate as every possible combination of the levels across predictors is compared. This may not be what you want. Instead, you can look for differences of one predictor based on a value of another. For example:
forComp <- pairs(emmOut, # emmeans output
adjust = 'fdr', # multiple comparison adjustment
simple = "Cat2") # contrasts within this categorical predictor
forCompCat3 = Treat1:
contrast estimate SE df t.ratio p.value
TypeA - TypeB -37.37 33.0 88 -1.133 0.2603
TypeA - TypeC 75.94 33.9 88 2.242 0.0412
TypeA - TypeD 141.97 38.3 88 3.704 0.0011
TypeB - TypeC 113.32 31.8 88 3.563 0.0012
TypeB - TypeD 179.34 36.5 88 4.912 <0.0001
TypeC - TypeD 66.02 37.3 88 1.769 0.0964
Cat3 = Treat2:
contrast estimate SE df t.ratio p.value
TypeA - TypeB -80.08 42.2 88 -1.895 0.1166
TypeA - TypeC 2.34 42.2 88 0.055 0.9560
TypeA - TypeD 75.94 39.6 88 1.916 0.1166
TypeB - TypeC 82.42 46.3 88 1.781 0.1166
TypeB - TypeD 156.03 43.9 88 3.554 0.0037
TypeC - TypeD 73.61 43.9 88 1.677 0.1166
Cat3 = Control:
contrast estimate SE df t.ratio p.value
TypeA - TypeB -118.62 28.0 88 -4.242 <0.0001
TypeA - TypeC 95.02 23.6 88 4.023 0.0001
TypeA - TypeD 213.69 25.7 88 8.299 <0.0001
TypeB - TypeC 213.64 27.0 88 7.918 <0.0001
TypeB - TypeD 332.31 28.9 88 11.514 <0.0001
TypeC - TypeD 118.67 24.7 88 4.809 <0.0001
P value adjustment: fdr method for 6 tests Here you can more easily see the contrasts, and the adjustment to the P-value will be limited to what you actually need. Based on a threshold P-value of 0.05, you can see that the effects at some combinations of your predictors are not statistically different from each other.
For example, the coefficient (predicted Resp2) when Cat3 = Treat1 and Cat2 = TypeA is 365 and this is 37 lower than the coefficient when Cat3 = Treat1 and Cat2 = TypeB, but there is no evidence that the difference has any meaning as P = 0.26 for this comparison.
On the other hand, some combinations of your predictors are statistically different from each other. For example, a comparison of modelled Resp2 when Cat3 = Control and Cat2 = TypeB (548) vs. Cat3 = Control and Cat2 = TypeD (215) shows that they differ by 332 and that there is strong evidence that this difference is meaningful as P < 0.0001.
Note that you could also organize this with contrasts by Cat3 instead:
forComp <- pairs(emmOut, # emmeans output
adjust = 'fdr', # multiple comparison adjustment
simple = "Cat3") # contrasts within this categorical predictor
forCompCat2 = TypeA:
contrast estimate SE df t.ratio p.value
Treat1 - Treat2 43.630 36.4 88 1.198 0.2341
Treat1 - Control -64.116 30.3 88 -2.116 0.0557
Treat2 - Control -107.745 31.9 88 -3.374 0.0033
Cat2 = TypeB:
contrast estimate SE df t.ratio p.value
Treat1 - Treat2 0.921 39.3 88 0.023 0.9814
Treat1 - Control -145.366 30.9 88 -4.711 <0.0001
Treat2 - Control -146.287 39.3 88 -3.719 0.0005
Cat2 = TypeC:
contrast estimate SE df t.ratio p.value
Treat1 - Treat2 -29.979 40.1 88 -0.748 0.6799
Treat1 - Control -45.038 28.1 88 -1.605 0.3363
Treat2 - Control -15.060 36.4 88 -0.414 0.6799
Cat2 = TypeD:
contrast estimate SE df t.ratio p.value
Treat1 - Treat2 -22.391 41.4 88 -0.541 0.8276
Treat1 - Control 7.609 34.8 88 0.218 0.8276
Treat2 - Control 30.000 34.8 88 0.861 0.8276
P value adjustment: fdr method for 3 tests Note that you can also get the results from the pairwise testing visually by handing the output of pairs() to plot().
Recall that Example 3 contains one predictor and the predictor is numeric:
Resp3 ~ Num4 + 1
For a numeric predictor, the effect is captured in the coefficient that describes the change in the response for a unit change in the numeric predictor.
For models with a normal error distribution assumption, this is communicated as the slope.
Let’s look at the summary() output for Example 3:
coef(summary(bestMod3)) # extract the coefficients from summary() Estimate Std. Error t value Pr(>|t|)
(Intercept) -226.9131 166.09360 -1.366176 1.750105e-01
Num4 260.0574 14.56593 17.853818 1.438227e-32This tells you that for a unit change in Num4, you get a 260 ± 1510 change in Resp3.
Note that this same information can be found using the emmeans package. Because Num4 is a numeric predictor, you use the emtrends() function from the emmeans package, which also provides the 95% confidence interval on the estimate:
library(emmeans) # load the emmeans package
trendsOut <- emtrends(bestMod3, # your model
specs = ~ Num4 + 1, # your predictors
var = "Num4") # your numeric predictors
trendsOut Num4 Num4.trend SE df lower.CL upper.CL
9.94 260 14.6 98 231 289
Confidence level used: 0.95 Recall that Example 4 contains two predictors, one is categorical and one is numeric:
Resp4 ~ Cat5 + Num6 + Cat5:Num6 + 1
What are your modelled effects (with uncertainty)?
The effects estimated for this model will include values of the response (Resp4) at each level of the categorical predictor (Cat5) as well as slopes describing the change in the response when the numeric predictor (Num6) changes, and a different slope will be estimated for each level in Cat5 (this is because of the interaction in the model).
As above, you can take a look at your modelled coefficients with the summary() output:
coef(summary(bestMod4)) Estimate Std. Error t value Pr(>|t|)
(Intercept) 98.34679337 1.868900753 52.62280151 1.536692e-71
Cat5Urban -0.23642403 2.895404981 -0.08165491 9.350947e-01
Cat5Wild -0.86733554 3.238562391 -0.26781498 7.894285e-01
Num6 -0.00293078 0.003740830 -0.78345704 4.353283e-01
Cat5Urban:Num6 0.01511745 0.005611511 2.69400763 8.361918e-03
Cat5Wild:Num6 0.03085821 0.006183238 4.99062262 2.756559e-06This shows:
The model prediction of Resp4 when Cat5 is Farm and Num6 is 0 is 98.34 units11 (the intercept).
The model prediction of Resp4 when Cat5 is Urban and Num6 is 0 is 98.34 - 0.24 = 98.1 units.
The model prediction of Resp4 when Cat5 is Wild and Num6 is 0 is 98.34 - 0.87 = 97.5 units.
The slope of the relationship between Num6 and Resp4 when Cat5 is Farm is -0.0029.12
The slope of the relationship between Num6 and Resp4 when Cat5 is Urban is -0.0029 + 0.015 = 0.0121.
The slope of the relationship between Num6 and Resp4 when Cat5 is Wild is -0.0029 + 0.031 = 0.0281.
Again, interpreting the coefficients from the summary() output is tedious and not necessary: You can use the emmeans package to give you the modelled response for each level of the categorical predictor (Cat5) directly:
emmOut <- emmeans(object = bestMod4, # your model
specs = ~ Cat5 + Num6 + Cat5:Num6 + 1, # your predictors
type = "response") # report coefficients on the response scale
emmOut Cat5 Num6 emmean SE df lower.CL upper.CL
Farm 510 96.9 0.458 94 95.9 97.8
Urban 510 104.3 0.421 94 103.5 105.2
Wild 510 111.7 0.511 94 110.7 112.7
Confidence level used: 0.95 Note that emmeans() sets our numeric predictor (Num6) to the mean value of Num6 (510 units). We can also control this with the at = argument:
emmOut <- emmeans(object = bestMod4, # your model
specs = ~ Cat5 + Num6 + Cat5:Num6 + 1, # your predictors
type = "response", # report coefficients on the response scale
at = list(Num6 = 0)) # control the value of your numeric predictor at which to make the coefficient estimates
emmOut Cat5 Num6 emmean SE df lower.CL upper.CL
Farm 0 98.3 1.87 94 94.6 102
Urban 0 98.1 2.21 94 93.7 103
Wild 0 97.5 2.64 94 92.2 103
Confidence level used: 0.95 By setting at = 0, you get the intercept - i.e. the modelled Resp4 when Num6 = 0 for each level of Cat5, and this is what is reported in the summary() output.
Similarly, you can get the emmeans package to give you the slope coefficients for the numeric predictor (Num6) using the emtrends() function, rather than interpreting them from the summary() output:
trendsOut <- emtrends(bestMod4, # your model
specs = ~ Cat5 + Num6 + Cat5:Num6 + 1, # your predictors
var = "Num6") # your numeric predictors
trendsOut Cat5 Num6 Num6.trend SE df lower.CL upper.CL
Farm 510 -0.00293 0.00374 94 -0.01036 0.0045
Urban 510 0.01219 0.00418 94 0.00388 0.0205
Wild 510 0.02793 0.00492 94 0.01815 0.0377
Confidence level used: 0.95 Note that the 95% confidence limits for the modelled effect for Num6 when Cat5 = Farm includes zero. We can also see this with the test() function.
test(trendsOut) # get test if coefficient is different than zero. Cat5 Num6 Num6.trend SE df t.ratio p.value
Farm 510 -0.00293 0.00374 94 -0.783 0.4353
Urban 510 0.01219 0.00418 94 2.914 0.0045
Wild 510 0.02793 0.00492 94 5.673 <0.0001combining the two output, we can report that:
there is no evidence of an effect on Resp4 of Num6 when Cat5 = Farm (P = 0.44).
there is evidence for a positive effect of Num6 on Resp4 when Cat5 = Urban (P = 0.0045). The effect is 0.01213 with a 95% confidence interval of 0.0039 - 0.020514.
there is evidence for a positive effect of Num6 on Resp4 when Cat5 = Wild (P < 0.0001). The effect is 0.02815 with a 95% confidence interval of 0.0182 - 0.037716.
Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)
Since you have a categorical predictor with more than two levels, you will want to report where the effects among levels differ from one another.
Note that the effects in Example 4 include slope coefficients (associated with the numeric predictor) and intercept coefficients (associated with the categorical predictor). Because the slope and intercept are dependent on one another (when the slope changes, the intercept will naturally change), you need to first check if the slope coefficients vary across your categorical levels. Only if they DO NOT vary (i.e. the modelled effects produce lines that are parallel across your categorical levels) would you go on to test if the intercepts vary across your categorical levels.
You can find out which slopes (i.e. the effect of Num6 on Resp4) are different across the levels of Cat5 using emtrends() with a request for pairwise testing:
forCompSlope <- pairs(trendsOut, # emmeans object
adjust = "fdr", # multiple comparison adjustment
simple = "Cat5") # contrasts within this categorical predictor
forCompSlopeNum6 = 510:
contrast estimate SE df t.ratio p.value
Farm - Urban -0.0151 0.00561 94 -2.694 0.0125
Farm - Wild -0.0309 0.00618 94 -4.991 <0.0001
Urban - Wild -0.0157 0.00646 94 -2.437 0.0167
P value adjustment: fdr method for 3 tests Here you see that the slope estimates (the effect of Num6 on Resp4) for all levels of Cat5 are significantly different from one another (0.0001 ≤ P ≤ 0.017).
Note that you can also get the results from the pairwise testing visually by handing the output of pairs() to plot().
Again: only if all the slopes were similar, would you want to test if the levels of your categorical predictor (Cat4) have significantly different coefficients (intercepts) from each other. You would do this with pairs() on the output from emmeans() (the emmOut object).
For a poisson error distribution assumption, the effect of a numeric predictor on the response is the percent change in the response for a unit change in the predictor.
The default (canonical) link function for models with a poisson error distribution assumption is the natural logarithm or link = "log".
If you are using link = "log" (as with a poisson error distribution, and sometimes also a Gamma error distribution - see below), you get your coefficient on the response scale by taking e to the coefficient on the link scale (with the R function exp()). This coefficient is called the rate ratio (RR)17 and it tells you the percent (%) change in the response for a unit change in your predictor.
But why do we convert coefficients from the link to response scale with \(e^x\)18 when link = "log"?
A reminder that we can present our linear model like this:
\[ \begin{align} g(\mu_i|Pred_i) &= \beta_1 \cdot Pred_i + \beta_0 \\ Resp_i &\sim F(\mu_i) \end{align} \]
where
\(g()\) is the link function
\(\mu_i\) is the mean fitted value \(\mu_i\) dependent on \(Pred_i\)
\(\beta_1\) is the coefficient of \(Pred_i\)
\(Pred_i\) is the value of your predictor
\(\beta_0\) is the intercept
\(Resp_i\) is the value of the response variable
\(F\) represents some error distribution assumption
If you have an error distribution assumption (\(F\)), the link function \(g()\) is \(log_e\)19:
\[ \begin{align} log_e(\mu_i|Pred_i) &= \beta_1 \cdot Pred_i + \beta_0 \\ Resp_i &\sim poisson(\mu_i) \end{align} \]
Then the rate ratio (RR), or % change in the response for a unit change in the predictor becomes,
\[ \begin{align} RR&=\frac{(\mu_i|Pred = a+1)}{(\mu_i|Pred = a)}\\[2em] log_e(RR)&=log_e(\frac{(\mu_i|Pred = a+1)}{(\mu_i|Pred = a)})\\[2em] log_e(RR)&=log_e(\mu_i|Pred = a+1)-log_e(\mu_i|Pred = a)\\[2em] log_e(RR)&=(\beta_1\cdot (a+1)+\beta_0)-(\beta_1\cdot (a)+\beta_0)\\[2em] log_e(RR)&=\beta_1\cdot a + \beta_1+\beta_0-\beta_1\cdot a-\beta_0\\[2em] log_e(RR)&=\beta_1\\[2em] RR &=e^{\beta_1} \end{align} \]
Here we have made new examples where the error distribution assumption is poisson. This changes the response variables, indicated by Resp1.p, Resp2.p, etc. to be integers (whole numbers) ≥ 0.
Recall that Example 1 contains only one predictor and it is categorical:
Resp1.p ~ Cat1 + 1
The best-specified model (now with poisson error distribution assumption) is:
bestMod1.p
Call: glm(formula = Resp1.p ~ Cat1 + 1, family = poisson(link = "log"),
data = Dat)
Coefficients:
(Intercept) Cat1StB Cat1StC
6.15524 -0.04869 0.09948
Degrees of Freedom: 49 Total (i.e. Null); 47 Residual
Null Deviance: 139
Residual Deviance: 40.01 AIC: 446that was fit to data in myDat1.p:
str(myDat1.p)'data.frame': 50 obs. of 2 variables:
$ Cat1 : Factor w/ 3 levels "StA","StB","StC": 2 2 2 1 3 2 3 1 3 1 ...
$ Resp1.p: int 434 476 443 473 516 449 529 462 512 429 ...and the dredge() table you used to pick your bestMod1.p is in dredgeOut1.p
dredgeOut1.pGlobal model call: glm(formula = Resp1.p ~ Cat1 + 1, family = poisson(link = "log"),
data = Dat)
---
Model selection table
(Intrc) Cat1 df logLik AICc delta weight
2 6.155 + 3 -219.983 446.5 0.00 1
1 6.164 1 -269.477 541.0 94.55 0
Models ranked by AICc(x) And here’s a visualization of the model effects:
# Set up your predictors for the visualized fit
forCat1<-unique(myDat1.p$Cat1) # every value of your categorical predictor
# create a data frame with your predictors
forVis<-expand.grid(Cat1=forCat1) # expand.grid() function makes sure you have all combinations of predictors
# Get your model fit estimates at each value of your predictors
modFit<-predict(bestMod1.p, # the model
newdata = forVis, # the predictor values
type = "link", # here, make the predictions on the link scale and we'll translate the back below
se.fit = TRUE) # include uncertainty estimate
ilink <- family(bestMod1.p)$linkinv # get the conversion (inverse link) from the model to translate back to the response scale
forVis$Fit<-ilink(modFit$fit) # add your fit to the data frame
forVis$Upper<-ilink(modFit$fit + 1.96*modFit$se.fit) # add your uncertainty to the data frame, 95% confidence intervals
forVis$Lower<-ilink(modFit$fit - 1.96*modFit$se.fit) # add your uncertainty to the data frame
library(ggplot2)
ggplot() + # start ggplot
geom_point(data = myDat1.p, # add observations to your plot
mapping = aes(x = Cat1, y = Resp1.p),
position=position_jitter(width=0.1)) + # control position of data points so they are easier to see on the plot
geom_errorbar(data = forVis, # add the uncertainty to your plot
mapping = aes(x = Cat1, y = Fit, ymin = Lower, ymax = Upper),
linewidth=1.2) + # control thickness of errorbar line
geom_point(data = forVis, # add the modelled fit to your plot
mapping = aes(x = Cat1, y = Fit),
shape = 22, # shape of point
size = 3, # size of point
fill = "white", # fill colour for plot
col = 'black') + # outline colour for plot
ylab("Resp1.p, (units)")+ # change y-axis label
xlab("Cat1, (units)")+ # change x-axis label
theme_bw() # change ggplot theme
What are your modelled effects (with uncertainty)?
For categorical predictors, you can get the coefficients on the response scale directly with the emmeans() function20 when you specify type = “response”:
library(emmeans) # load the emmeans package
emmOut <- emmeans(object = bestMod1.p, # your model
specs = ~ Cat1 + 1, # your predictors
type = "response") # report coefficients on the response scale
emmOut Cat1 rate SE df asymp.LCL asymp.UCL
StA 471 6.54 Inf 459 484
StB 449 4.32 Inf 440 457
StC 520 5.89 Inf 509 532
Confidence level used: 0.95
Intervals are back-transformed from the log scale Notice that the column reported is “rate”: these are the coefficients converted back onto the response scale.
Note also that two types of uncertainty are measured here. SE stands for the standard error around the prediction, and is a measure of uncertainty of the average modelled effect. The lower.CL and upper.CL represent the 95% confidence limits of the prediction - so if I observed a new Resp1.p at a particular Cat1, there would be a 95% chance it would lie between the bounds of the confidence limits.
From this output you see:
When Cat1 is SpA, Resp1.p is estimated to be 471 (95% confidence interval: 459 to 484 units).
When Cat1 is SpB, Resp1.p is estimated to be 449 (95% confidence interval: 440 to 457 units).
When Cat1 is SpC, Resp1.p is estimated to be 520 (95% confidence interval: 509 to 532 units).
Finally, note that you can also get a quick plot of the effects by handing the emmeans() output to plot().
Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)
Here you can compare the effects of the levels of Cat1 on Resp1.p to see if the effects differ between the levels.
forComp <- pairs(emmOut, #emmeans object
adjust = "fdr") # multiple testing adjustment
forComp contrast ratio SE df null z.ratio p.value
StA / StB 1.050 0.0177 Inf 1 2.880 0.0040
StA / StC 0.905 0.0162 Inf 1 -5.552 <0.0001
StB / StC 0.862 0.0128 Inf 1 -9.968 <0.0001
P value adjustment: fdr method for 3 tests
Tests are performed on the log scale The output shows the results of the multiple comparison (pairwise) testing.
The values in the p.value column tell you the results of the hypothesis test comparing the coefficients between the two levels of your categorical predictor.
For example:
there is evidence that Resp1.p is significantly larger when Cat 1 = SpA than when Cat1 = SpB (p = 0.004). Resp1.p is expected to be ~ 5 ± 2% bigger when Cat 1 = SpA than when Cat 1 = SpB.
there is evidence that Resp1.p is significantly smaller when Cat 1 = SpA than when Cat1 = SpC (p = 0). Resp1.p is expected to be ~ 9 ± 2% smaller when Cat 1 = SpA than when Cat 1 = SpC.
Note that you can also get the results from the pairwise testing visually by handing the output of pairs() to plot().
Recall that Example 2 contains two predictors and both are categorical:
Resp2.p ~ Cat2 + Cat3 + Cat2:Cat3 + 1
The best-specified model (now with poisson error distribution assumption) is:
bestMod2.p
Call: glm(formula = Resp2.p ~ Cat2 + Cat3 + Cat2:Cat3 + 1, family = poisson(link = "log"),
data = myDat2.p)
Coefficients:
(Intercept) Cat2TypeB Cat2TypeC
2.7213 -0.3365 -0.4526
Cat2TypeD Cat3Treat Cat2TypeB:Cat3Treat
-1.2744 -0.3810 -0.6176
Cat2TypeC:Cat3Treat Cat2TypeD:Cat3Treat
-0.5999 0.4381
Degrees of Freedom: 49 Total (i.e. Null); 42 Residual
Null Deviance: 132.5
Residual Deviance: 35.27 AIC: 239.5that was fit to data in myDat2.p:
str(myDat2.p)'data.frame': 50 obs. of 3 variables:
$ Cat2 : Factor w/ 4 levels "TypeA","TypeB",..: 1 2 4 1 2 4 1 1 2 2 ...
$ Cat3 : Factor w/ 2 levels "Control","Treat": 2 2 1 1 2 1 1 1 1 1 ...
$ Resp2.p: int 10 3 3 19 4 4 9 13 15 12 ...and the dredge() table you used to pick your bestMod2.p is in
dredgeOut2.pGlobal model call: glm(formula = Resp2.p ~ Cat2 + Cat3 + Cat2:Cat3 + 1, family = poisson(link = "log"),
data = myDat2.p)
---
Model selection table
(Int) Ct2 Ct3 Ct2:Ct3 df logLik AICc delta weight
8 2.721 + + + 8 -111.758 243.0 0.00 0.77
4 2.860 + + 5 -117.038 245.4 2.41 0.23
2 2.461 + 4 -134.275 277.4 34.41 0.00
3 2.344 + 2 -150.075 304.4 61.38 0.00
1 2.087 1 -160.376 322.8 79.81 0.00
Models ranked by AICc(x) And here’s a visualization of the model effects:
#### i) choosing the values of your predictors at which to make a prediction
# Set up your predictors for the visualized fit
forCat2<-unique(myDat2.p$Cat2) # every level of your categorical predictor
forCat3<-unique(myDat2.p$Cat3) # every level of your categorical predictor
# create a data frame with your predictors
forVis<-expand.grid(Cat2 = forCat2, Cat3 = forCat3) # expand.grid() function makes sure you have all combinations of predictors
#### ii) using `predict()` to use your model to estimate your response variable at those values of your predictor
# Get your model fit estimates at each value of your predictors
modFit<-predict(bestMod2.p, # the model
newdata = forVis, # the predictor values
type = "link", # here, make the predictions on the link scale and we'll translate the back below
se.fit = TRUE) # include uncertainty estimate
ilink <- family(bestMod2.p)$linkinv # get the conversion (inverse link) from the model to translate back to the response scale
forVis$Fit<-ilink(modFit$fit) # add your fit to the data frame
forVis$Upper<-ilink(modFit$fit + 1.96*modFit$se.fit) # add your uncertainty to the data frame, 95% confidence intervals
forVis$Lower<-ilink(modFit$fit - 1.96*modFit$se.fit) # add your uncertainty to the data frame
#### iii) use the model estimates to plot your model fit
library(ggplot2) # load ggplot2 library
ggplot() + # start ggplot
geom_point(data = myDat2.p, # add observations to your plot
mapping = aes(x = Cat2, y = Resp2.p, col = Cat3),
position=position_jitterdodge(jitter.width=0.75, dodge.width=0.75)) + # control position of data points so they are easier to see on the plot
geom_errorbar(data = forVis, # add the uncertainty to your plot
mapping = aes(x = Cat2, y = Fit, ymin = Lower, ymax = Upper, col = Cat3),
position=position_dodge(width=0.75), # control position of data points so they are easier to see on the plot
size=1.2) + # control thickness of errorbar line
geom_point(data = forVis, # add the modelled fit to your plot
mapping = aes(x = Cat2, y = Fit, fill = Cat3),
shape = 22, # shape of point
size=3, # size of point
col = 'black', # outline colour for point
position=position_dodge(width=0.75)) + # control position of data points so they are easier to see on the plot
ylab("Resp2.p, (units)")+ # change y-axis label
xlab("Cat2, (units)")+ # change x-axis label
labs(fill="Cat3, (units)", col="Cat3, (units)") + # change legend title
theme_bw() # change ggplot theme
What are your modelled effects (with uncertainty)?
As above, we can use the emmeans package to see the effects of the predictors on the scale of the response:
emmOut <- emmeans(object = bestMod2.p, # your model
specs = ~ Cat2 + Cat3 + Cat2:Cat3, # your predictors
type = "response") # report coefficients on the response scale
emmOut Cat2 Cat3 rate SE df asymp.LCL asymp.UCL
TypeA Control 15.20 1.740 Inf 12.14 19.03
TypeB Control 10.86 1.250 Inf 8.67 13.59
TypeC Control 9.67 1.800 Inf 6.72 13.91
TypeD Control 4.25 1.030 Inf 2.64 6.84
TypeA Treat 10.38 0.894 Inf 8.77 12.29
TypeB Treat 4.00 0.707 Inf 2.83 5.66
TypeC Treat 3.62 0.673 Inf 2.52 5.22
TypeD Treat 4.50 1.500 Inf 2.34 8.65
Confidence level used: 0.95
Intervals are back-transformed from the log scale So now we have a modelled value of our response for each level of our categorical predictors. For example:
The modelled prediction for Resp2.p when Cat2 is TypeA and Cat3 is Treat is 15.2 (95% confidence interval: 12.1 - 19.0 units).
The modelled prediction for Resp2.p when Cat2 is TypeC and Cat3 is Control is 9.7 units (95% confidence interval: 6.7 - 13.9 units).
Finally, note that you can also get a quick plot of the effects by handing the emmeans() output to plot().
Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)
As with Example 1, you can also find out which combinations of predictor levels are leading to statistically different model predictions in Resp2.p:
forComp <- pairs(emmOut, # emmeans output
adjust = "fdr", # multiple testing adjustment
simple = "Cat2") # contrasts within this categorical predictor
forCompCat3 = Control:
contrast ratio SE df null z.ratio p.value
TypeA / TypeB 1.400 0.227 Inf 1 2.074 0.0457
TypeA / TypeC 1.572 0.343 Inf 1 2.074 0.0457
TypeA / TypeD 3.576 0.960 Inf 1 4.750 <0.0001
TypeB / TypeC 1.123 0.245 Inf 1 0.532 0.5947
TypeB / TypeD 2.555 0.685 Inf 1 3.496 0.0014
TypeC / TypeD 2.275 0.695 Inf 1 2.690 0.0143
Cat3 = Treat:
contrast ratio SE df null z.ratio p.value
TypeA / TypeB 2.596 0.510 Inf 1 4.852 <0.0001
TypeA / TypeC 2.865 0.586 Inf 1 5.142 <0.0001
TypeA / TypeD 2.308 0.794 Inf 1 2.429 0.0303
TypeB / TypeC 1.103 0.283 Inf 1 0.384 0.7549
TypeB / TypeD 0.889 0.335 Inf 1 -0.312 0.7549
TypeC / TypeD 0.806 0.307 Inf 1 -0.567 0.7549
P value adjustment: fdr method for 6 tests
Tests are performed on the log scale Here you can more easily see the contrasts. Contrast ratios similar to 1 mean that the predictor levels lead to a similar value of the response.
Based on a threshold P-value of 0.05, you can see that the effects at some combinations of your predictors are not statistically different from each other. For example, the coefficient (predicted Resp2.p) when Cat3 = Control and Cat2 = TypeB (10.9) is 1.12 times (or 12%) more than the predicted Resp2.p when Cat3 = Control and Cat2 = TypeC (9.7), but there is no evidence that the difference has any meaning as P = 0.59 for this comparison.
On the other hand, some combinations of your predictors are statistically different from each other. For example, a comparison of modelled Resp2.p when Cat3 = Control and Cat2 = TypeB (10.9) is 2.6 times (or 160%) higher than when Cat3 = Control and Cat2 = TypeD (4.2) and there is strong evidence that this difference is meaningful as P < 0.0014.
Note that you can also get the results from the pairwise testing visually by handing the output of pairs() to plot().
Recall that Example 3 contains one predictor and the predictor is numeric:
Resp3 ~ Num4 + 1
The best-specified model (now with poisson error distribution assumption) is:
bestMod3.p
Call: glm(formula = Resp ~ Num4 + 1, family = poisson(link = "log"),
data = myDat3.p)
Coefficients:
(Intercept) Num4
1.577204 0.006674
Degrees of Freedom: 49 Total (i.e. Null); 48 Residual
Null Deviance: 74.74
Residual Deviance: 62.96 AIC: 256.7that was fit to data in myDat3.p:
str(myDat3.p)'data.frame': 50 obs. of 2 variables:
$ Num4 : num 67.5 61.9 71.5 116.6 55.9 ...
$ Resp3.p: int 12 8 4 12 0 2 4 10 8 7 ...and the dredge() table you used to pick your bestMod3.p is
dredgeOut3.pGlobal model call: glm(formula = Resp ~ Num4 + 1, family = poisson(link = "log"),
data = myDat3.p)
---
Model selection table
(Intrc) Num4 df logLik AICc delta weight
2 1.577 0.006674 2 -126.339 256.9 0.00 0.992
1 2.069 1 -132.230 266.5 9.61 0.008
Models ranked by AICc(x) And here’s a visualization of the model effects:
#### i) choosing the values of your predictors at which to make a prediction
# Set up your predictors for the visualized fit
forNum4<-seq(from = min(myDat3.p$Num4), to = max(myDat3.p$Num4), by = 1)# a sequence of numbers in your numeric predictor range
# create a data frame with your predictors
forVis<-expand.grid(Num4 = forNum4) # expand.grid() function makes sure you have all combinations of predictors.
#### ii) using `predict()` to use your model to estimate your response variable at those values of your predictor
# Get your model fit estimates at each value of your predictors
modFit<-predict(bestMod3.p, # the model
newdata = forVis, # the predictor values
type = "link", # here, make the predictions on the link scale and we'll translate the back below
se.fit = TRUE) # include uncertainty estimate
ilink <- family(bestMod3.p)$linkinv # get the conversion (inverse link) from the model to translate back to the response scale
forVis$Fit<-ilink(modFit$fit) # add your fit to the data frame
forVis$Upper<-ilink(modFit$fit + 1.96*modFit$se.fit) # add your uncertainty to the data frame, 95% confidence intervals
forVis$Lower<-ilink(modFit$fit - 1.96*modFit$se.fit) # add your uncertainty to the data frame
#### iii) use the model estimates to plot your model fit
library(ggplot2) # load ggplot2 library
ggplot() + # start ggplot
geom_point(data = myDat3.p, # add observations to your plot
mapping = aes(x = Num4, y = Resp3.p)) + # control position of data points so they are easier to see on the plot
geom_line(data = forVis, # add the modelled fit to your plot
mapping = aes(x = Num4, y = Fit),
linewidth = 1.2) + # control thickness of line
geom_line(data = forVis, # add uncertainty to your plot (upper line)
mapping = aes(x = Num4, y = Upper),
linewidth = 0.8, # control thickness of line
linetype = 2) + # control style of line
geom_line(data = forVis, # add uncertainty to your plot (lower line)
mapping = aes(x = Num4, y = Lower),
linewidth = 0.8, # control thickness of line
linetype = 2) + # control style of line
ylab("Resp3.p, (units)") + # change y-axis label
xlab("Num4, (units)") + # change x-axis label
theme_bw() # change ggplot theme
trendsOut <- emtrends(bestMod3.p,
specs = ~ Num4 + 1, # your predictors
var = "Num4") # your predictor of interest
trendsOut Num4 Num4.trend SE df asymp.LCL asymp.UCL
71.5 0.00667 0.00193 Inf 0.00288 0.0105
Confidence level used: 0.95 Note that these are the same estimates as we find in the summary() of your model object
coef(summary(bestMod3.p)) Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.577203770 0.155306248 10.155443 3.133832e-24
Num4 0.006674037 0.001934204 3.450534 5.594789e-04This is an estimate of the modelled effects in the linked world. We need to consider our log link function to report the modelled effects in the response world.
Since you have a poisson error distribution assumption, you can convert the estimate made by emtrends() on to the response scale with:
trendsOut.df <- data.frame(trendsOut) # convert trendsOut to a dataframe
RR <- exp(trendsOut.df$Num4.trend) # get the coefficient on the response scale
RR[1] 1.006696We can also use the same conversion on the confidence limits of the modelled effect, for the upper:
RR.up <- exp(trendsOut.df$asymp.UCL) # get the upper confidence interval on the response scale
RR.up[1] 1.01052and lower 95% confidence level:
RR.low <- exp(trendsOut.df$asymp.LCL) # get the lower confidence interval on the response scale
RR.low[1] 1.002887This tells you that for a unit change in Num4, Resp3 increases by 1.007x (95% confidence interval is 1.003 to 1.011) which is an increase of 0.67%.
Recall that Example 4 contains two predictors, one is categorical and one is numeric:
Resp4 ~ Cat5 + Num6 + Cat5:Num6 + 1
The best-specified model (now with poisson error distribution assumption) is:
bestMod4.p
Call: glm(formula = Resp4.p ~ Cat5 + Num6 + Cat5:Num6 + 1, family = poisson(link = "log"),
data = myDat4.p)
Coefficients:
(Intercept) Cat5StB Cat5StC Num6 Cat5StB:Num6
4.374118 0.074861 -0.262595 0.004075 -0.002618
Cat5StC:Num6
0.003898
Degrees of Freedom: 49 Total (i.e. Null); 44 Residual
Null Deviance: 183.1
Residual Deviance: 50.35 AIC: 386.8that was fit to data in myDat4.p:
str(myDat4.p)'data.frame': 50 obs. of 3 variables:
$ Cat5 : Factor w/ 3 levels "StA","StB","StC": 2 1 3 2 2 1 1 3 3 2 ...
$ Num6 : num 60.9 55.9 82.8 89.4 53.4 ...
$ Resp4.p: int 85 112 109 87 99 98 92 101 87 110 ...and the dredge() table you used to pick your bestMod4.p is
dredgeOut4.pGlobal model call: glm(formula = Resp4.p ~ Cat5 + Num6 + Cat5:Num6 + 1, family = poisson(link = "log"),
data = myDat4.p)
---
Model selection table
(Int) Ct5 Nm6 Ct5:Nm6 df logLik AICc delta weight
8 4.374 + 0.004075 + 6 -187.378 388.7 0.00 0.999
4 4.291 + 0.005194 4 -196.412 401.7 13.00 0.001
3 4.261 0.005378 2 -207.753 419.8 31.05 0.000
2 4.669 + 3 -239.402 485.3 96.62 0.000
1 4.665 1 -253.742 509.6 120.86 0.000
Models ranked by AICc(x) And here’s a visualization of the model effects:
#### i) choosing the values of your predictors at which to make a prediction
# Set up your predictors for the visualized fit
forCat5<-unique(myDat4.p$Cat5) # all levels of your categorical predictor
forNum6<-seq(from = min(myDat4.p$Num6), to = max(myDat4.p$Num6), by = 1)# a sequence of numbers in your numeric predictor range
# create a data frame with your predictors
forVis<-expand.grid(Cat5 = Cat5, Num6 = forNum6) # expand.grid() function makes sure you have all combinations of predictors.
#### ii) using `predict()` to use your model to estimate your response variable at those values of your predictor
# Get your model fit estimates at each value of your predictors
modFit<-predict(bestMod4.p, # the model
newdata = forVis, # the predictor values
type = "link", # here, make the predictions on the link scale and we'll translate the back below
se.fit = TRUE) # include uncertainty estimate
ilink <- family(bestMod4.p)$linkinv # get the conversion (inverse link) from the model to translate back to the response scale
forVis$Fit<-ilink(modFit$fit) # add your fit to the data frame
forVis$Upper<-ilink(modFit$fit + 1.96*modFit$se.fit) # add your uncertainty to the data frame, 95% confidence intervals
forVis$Lower<-ilink(modFit$fit - 1.96*modFit$se.fit) # add your uncertainty to the data frame
#### iii) use the model estimates to plot your model fit
library(ggplot2) # load ggplot2 library
ggplot() + # start ggplot
geom_point(data = myDat4.p, # add observations to your plot
mapping = aes(x = Num6, y = Resp4.p, col = Cat5)) + # control position of data points so they are easier to see on the plot
geom_line(data = forVis, # add the modelled fit to your plot
mapping = aes(x = Num6, y = Fit, col = Cat5),
linewidth = 1.2) + # control thickness of line
geom_line(data = forVis, # add uncertainty to your plot (upper line)
mapping = aes(x = Num6, y = Upper, col = Cat5),
linewidth = 0.4, # control thickness of line
linetype = 2) + # control style of line
geom_line(data = forVis, # add uncertainty to your plot (lower line)
mapping = aes(x = Num6, y = Lower, col = Cat5),
linewidth = 0.4, # control thickness of line
linetype = 2) + # control style of line
ylab("Resp4, (units)") + # change y-axis label
xlab("Num6, (units)") + # change x-axis label
labs(fill="Cat5, (units)", col="Cat5, (units)") + # change legend title
theme_bw() # change ggplot theme
What are your modelled effects (with uncertainty)?
The process for estimating effects of categorical vs. numeric predictors is different.
For categorical predictors (Cat5), you can use emmeans() to give you the effects on the response scale directly:
emmOut <- emmeans(object = bestMod4.p, # your model
specs = ~ Cat5 + Num6 + Cat5:Num6 + 1, # your predictors
type = "response") # report coefficients on the response scale
emmOut Cat5 Num6 rate SE df asymp.LCL asymp.UCL
StA 73.5 107.1 2.32 Inf 102.7 112
StB 73.5 95.2 2.61 Inf 90.2 100
StC 73.5 109.7 2.70 Inf 104.5 115
Confidence level used: 0.95
Intervals are back-transformed from the log scale When Cat5 is StA, Resp4.p is estimated to be 74 (95% confidence interval: to 103 units).
When Cat5 is StB, Resp4.p is estimated to be 74 (95% confidence interval: to 90 units). When Cat5 is StC, Resp4.p is estimated to be 74 (95% confidence interval: to 105 units).
For numeric predictors (Num6), you need to use the emtrends() function and then convert the coefficients to the response scale:
trendsOut <- emtrends(bestMod4.p,
specs = ~ Cat5 + Num6 + Cat5:Num6 + 1, # your predictors
var = "Num6") # your predictor of interest
trendsOut Cat5 Num6 Num6.trend SE df asymp.LCL asymp.UCL
StA 73.5 0.00408 0.000829 Inf 0.00245 0.00570
StB 73.5 0.00146 0.001450 Inf -0.00139 0.00430
StC 73.5 0.00797 0.000903 Inf 0.00620 0.00974
Confidence level used: 0.95 Since you have a poisson error distribution assumption, you can convert the estimate made by emtrends() on to the response scale with the exp() function:
trendsOut.df <- data.frame(trendsOut) # convert trendsOut to a dataframe
RR <- exp(trendsOut.df$Num6.trend) # get the coefficient on the response scale
RR[1] 1.004084 1.001458 1.008005You can also use the same conversion on the confidence limits of the modelled effect, for the upper:
RR.up <- exp(trendsOut.df$asymp.UCL) # get the upper confidence interval on the response scale
RR.up[1] 1.005716 1.004309 1.009791and lower 95% confidence level:
RR.low <- exp(trendsOut.df$asymp.LCL) # get the lower confidence interval on the response scale
RR.low[1] 1.0024542 0.9986158 1.0062224Let’s organize these numbers so we can read the effects easier:
RRAll <- data.frame(Cat5 = trendsOut.df$Cat5, # include info on the Cat5 level
RR = exp(trendsOut.df$Num6.trend), # the modelled effect as a rate ratio
RR.up = exp(trendsOut.df$asymp.UCL), # upper 95% confidence level of the modelled effect
RR.down = exp(trendsOut.df$asymp.LCL)) # lower 95% confidence level of the modelled effect
RRAll Cat5 RR RR.up RR.down
1 StA 1.004084 1.005716 1.0024542
2 StB 1.001458 1.004309 0.9986158
3 StC 1.008005 1.009791 1.0062224This tells you that for a unit change in Num6, Resp4.p
increases by 1.004 times (an increase of 0.41%) when Cat5 is StA (95% confidence interval: 1.002 to 1.006).
increases by 1.001 times (an increase of 0.15%) when Cat5 is StB (95% confidence interval: 0.999 to 1.004).
increases by 1.008 times (an increase of 0.8%) when Cat5 is StC (95% confidence interval: 1.006 to 1.01).
Note that the 95% confidence interval for the effect of Num6 on Resp4.p (the rate ratio) includes 1. When the rate ratio is 1, there is no effect of the numeric predictor on the response.
You can test for evidence that the effects of Num6 on Resp4.p are significant with:
test(trendsOut) # get test if coefficient is different than zero. Cat5 Num6 Num6.trend SE df z.ratio p.value
StA 73.5 0.00408 0.000829 Inf 4.918 <0.0001
StB 73.5 0.00146 0.001450 Inf 1.005 0.3150
StC 73.5 0.00797 0.000903 Inf 8.828 <0.0001Note that these tests are done on the link scale but can be used for reporting effects on the response scale.
From the results, you can see that:
there is evidence that there is an effect of Num6 on Resp4.p when Cat5 = StA. (i.e. the slope on the link scale is different than zero, and the rate ratio on the response scale is different than 1; P < 0.0001).
there is no evidence that there is an effect of Num6 on Resp4.p when Cat5 = StB. (i.e. the slope on the link scale is not different than zero, and the rate ratio on the response scale is not different than 1; P = 0.32).
there is evidence that there is an effect of Num6 on Resp4.p when Cat5 = StC. (i.e. the slope on the link scale is different than zero, and the rate ratio on the response scale is different than 1; P < 0.0001).
Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)
You can also find out which effects of Num6 on Resp4.p are different from one another using pairs() on the emtrends() output:
forCompRR <- pairs(trendsOut, # emmeans object
adjust = "fdr", # multiple comparison adjustment
simple = "Cat5") # contrasts within this categorical predictor
forCompRRNum6 = 73.5:
contrast estimate SE df z.ratio p.value
StA - StB 0.00262 0.00167 Inf 1.568 0.1170
StA - StC -0.00390 0.00123 Inf -3.180 0.0022
StB - StC -0.00652 0.00171 Inf -3.814 0.0004
P value adjustment: fdr method for 3 tests Remember that you need to convert any predictor coefficients to the response scale when reporting.
The table above tells you that:
There is no evidence the effect of Num6 on Resp4.p differs when Cat5 is StA vs. StB (P = 0.12).
There is evidence that the effect of Num6 on Resp4.p is different when Cat5 is StA vs. StC (P = 0.0022). The effect of Num6 on Resp4.p is exp(-0.0039) = 0.996 times as big when Cat5 is StA vs. StC.
There is evidence that the effect of Num6 on Resp4.p is different when Cat5 is StB vs. StC (P = 0.0004). The effect of Num6 on Resp4.p is exp(-0.00652) = 0.994 times as big when Cat5 is StB vs. StC.
Only if all the slopes were similar: you would want to test if the levels of your categorical predictor (Cat5) have significantly different coefficients (intercepts) from each other with pairs() on the output from emmeans() (the emmOut object).
For models with a Gamma error distribution assumption, the default link function is the inverse or link = "inverse". You need to take this link function into consideration when you report your modelled effects.
For a model where link = "inverse", the interpretation of the coefficients on the response scale are difficult to interpret.
Because of this, you can fit your GLM with a log link function (link = "log") similar to your models with a poisson error distribution assumption (and validate your model). It is much easier to interpret coefficients with a log (vs. inverse) link function - you can use the methods described above for the poisson error distribution assumption.
Alternatively, if you decide to use an inverse link function for your GLM21, you can report your modelled effects by visualizing your model and use your model to make predictions to give examples of your model’s effects.
For a binomial error distribution assumption, the effect of a numeric predictor on the response is the percent change in the odds for a unit change in the predictor.
Recall that models with a binomial error distribution assumption model the probability of an event happening. This might be presence (vs. absence), mature (vs. immature), death (vs. alive), etc., but in all cases this is represented as the probability of your response being 1 (vs. 0).
If \(p\) is the probability of an event happening (e.g. response being 1), the probability an event doesn’t happen is \(1-p\) when an event can only have two outcomes. The odds are the probability of the event happening divided by the probability it does not happen, or:
\(\frac{p}{1-p}\)
For models with a binomial error distribution assumption, you will typically use the logit function or link = "logit" for your link function22. The logit function is also called the “log-odds” function as it is equal to the logarithm of the odds:
\(log_e(\frac{p}{1-p})\)
Alternate link functions for binomially distributed data are probit and cloglog, but the logit link function is by far the one most used as the coefficients that result from the model when using a logit link function are the easiest to interpret. But, as always, you should choose the link function that leads to the best-behaved residuals for your model (See more in the Model Validation section).
You report your coefficients of a model with a binomial error distribution assumption (and link = “logit”) on the response scale by raising e to the power of the coefficients estimated on the link scale.
Here is the math that connects the coefficient to the odds ratio:
Given a binomial model, \[ \begin{align} log_e(\frac{p_i}{1-p_i}) &= \beta_1 \cdot Pred_i + \beta_0 \\ Resp_i &\sim binomial(p_i) \end{align} \] \(p_i\) is the probability of success and \(1-p_i\) is the probability of failure, and the odds are \(\frac{p_i}{1-p_i}\).
Then the odds ratio (OR), or % change in odds due to a change in predictor is:
\[ \begin{align} OR&=\frac{odds\;when\;Pred=a+1}{odds\;when \;Pred=a}\\[2em] OR&=\frac{(\frac{p}{1-p}|Pred=a+1)}{(\frac{p}{1-p}|Pred=a)}\\[2em] log_e(OR)&=(log_e\frac{p}{1-p}|Pred=a+1) - (log_e\frac{p}{1-p}|Pred=a)\\[2em] log_e(OR)&=(\beta_1\cdot (a+1)+\beta_0) - (\beta_1\cdot a+\beta_0)\\[2em] log_e(OR)&=\beta_1\cdot a + \beta_1-\beta_1\cdot a\\[2em] OR &=e^{\beta_1} \end{align} \]
We have made new examples where the error distribution assumption is binomial This changes the response variables, indicated by Resp1.b, Resp2.b, etc. to be values of either 0 or 1.
Recall that Example 1 contains only one predictor and it is categorical:
Resp1.b ~ Cat1 + 1
The best-specified model (now with a binomial error distribution assumption) is:
bestMod1.b
Call: glm(formula = Resp1.b ~ Cat1 + 1, family = binomial(link = "logit"),
data = myDat1)
Coefficients:
(Intercept) Cat1StB Cat1StC
-0.6931 -1.8718 1.6094
Degrees of Freedom: 49 Total (i.e. Null); 47 Residual
Null Deviance: 68.03
Residual Deviance: 51.43 AIC: 57.43that was fit to data in myDat1.b:
str(myDat1.b)'data.frame': 50 obs. of 2 variables:
$ Cat1 : Factor w/ 3 levels "StA","StB","StC": 3 3 2 3 3 1 1 2 3 3 ...
$ Resp1.b: int 1 0 0 1 1 0 0 0 1 0 ...and the dredge() table you used to pick your bestMod1.b is in dredgeOut1.b
dredgeOut1.bGlobal model call: glm(formula = Resp1.b ~ Cat1 + 1, family = binomial(link = "logit"),
data = myDat1)
---
Model selection table
(Intrc) Cat1 R^2 df logLik AICc delta weight
2 -0.6931 + 0.2825 3 -25.714 57.9 0.00 0.998
1 -0.3228 0.0000 1 -34.015 70.1 12.16 0.002
Models ranked by AICc(x) And here’s a visualization of the model effects:
# Set up your predictors for the visualized fit
forCat1<-unique(myDat1.b$Cat1) # every value of your categorical predictor
# create a data frame with your predictors
forVis<-expand.grid(Cat1=forCat1) # expand.grid() function makes sure you have all combinations of predictors
# Get your model fit estimates at each value of your predictors
modFit<-predict(bestMod1.b, # the model
newdata = forVis, # the predictor values
type = "link", # here, make the predictions on the link scale and we'll translate the back below
se.fit = TRUE) # include uncertainty estimate
ilink <- family(bestMod1.b)$linkinv # get the conversion (inverse link) from the model to translate back to the response scale
forVis$Fit<-ilink(modFit$fit) # add your fit to the data frame
forVis$Upper<-ilink(modFit$fit + 1.96*modFit$se.fit) # add your uncertainty to the data frame, 95% confidence intervals
forVis$Lower<-ilink(modFit$fit - 1.96*modFit$se.fit) # add your uncertainty to the data frame
library(ggplot2)
ggplot() + # start ggplot
geom_point(data = myDat1.b, # add observations to your plot
mapping = aes(x = Cat1, y = Resp1.b)) +
geom_errorbar(data = forVis, # add the uncertainty to your plot
mapping = aes(x = Cat1, y = Fit, ymin = Lower, ymax = Upper),
linewidth=1.2) + # control thickness of errorbar line
geom_point(data = forVis, # add the modelled fit to your plot
mapping = aes(x = Cat1, y = Fit),
shape = 22, # shape of point
size = 3, # size of point
fill = "white", # fill colour for plot
col = 'black') + # outline colour for plot
ylab("probability Resp1.b = 1")+ # change y-axis label
xlab("Cat1, (units)")+ # change x-axis label
theme_bw() # change ggplot theme
What are your modelled effects (with uncertainty)?
For categorical predictors, you can get the coefficients on the response scale directly with the emmeans() function23 when you specify type = “response”. These are the probabilities (not odds or odds ratio) that Resp1.b = 1 at each level of your categorical predictor:
library(emmeans) # load the emmeans package
emmOut <- emmeans(object = bestMod1.b, # your model
specs = ~ Cat1, # your predictors
type = "response") # report coefficients on the response scale
emmOut Cat1 prob SE df asymp.LCL asymp.UCL
StA 0.3333 0.1220 Inf 0.14596 0.594
StB 0.0714 0.0688 Inf 0.00996 0.370
StC 0.7143 0.0986 Inf 0.49239 0.866
Confidence level used: 0.95
Intervals are back-transformed from the logit scale Notice that the column reported is prob: these coefficients are already converted back onto the response scale (as probabilities).
Note also that two types of uncertainty are measured here. SE stands for the standard error around the prediction, and is a measure of uncertainty of the average modelled effect. The lower.CL and upper.CL represent the 95% confidence limits of the prediction - so if you observed a new Resp1.b at a particular Cat1, there is a 95% probability that the probability Resp1.b is 1 would lie within the bounds of the confidence limits.
From this output you see:
When Cat1 is SpA, the probability that Resp1.b is 1 is estimated to be 0.33 (95% confidence interval: 0.15 to 0.59 units). This is equivalent to odds (\(\frac{p}{(1-p)}\)) of 0.5.
When Cat1 is SpB, the probability that Resp1.b is 1 is estimated to be 0.07 (95% confidence interval: 0.01 to 0.37 units). This is equivalent to odds (\(\frac{p}{(1-p)}\)) of 0.08.
When Cat1 is SpC, the probability that Resp1.b is 1 is estimated to be 0.71 (95% confidence interval: 0.49 to 0.87 units). This is equivalent to odds (\(\frac{p}{(1-p)}\)) of 2.5.
Finally, note that you can also get a quick plot of the effects by handing the emmeans() output to plot().
The emmeans package provides an easy way of seeing the effects of Cat1 on Resp1.b the response scale, but you can also get this information from the summary() output:
coef(summary(bestMod1.b)) Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.6931472 0.5477226 -1.265508 0.20568935
Cat1StB -1.8718022 1.1734223 -1.595165 0.11067536
Cat1StC 1.6094379 0.7302967 2.203814 0.02753745You can convert the coefficients to odds with the exp() function, for example:
When Cat1 is SpA, the odds that Resp1.b is 1 is estimated to be exp(coef(summary(bestMod1.b))[1]) = 0.5. And you can get the probabilities directly with the plogis() function (plogis(coef(summary(bestMod1.b))[1]) = 0.333).
Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)
Here you can compare the effects of the levels of Cat1 on Resp1.b to see if the effects differ between the levels.
forComp <- pairs(emmOut, # emmeans object
adjust = "fdr") # multiple comparison adjustment
forComp contrast odds.ratio SE df null z.ratio p.value
StA / StB 6.5000 7.6300 Inf 1 1.595 0.1107
StA / StC 0.2000 0.1460 Inf 1 -2.204 0.0413
StB / StC 0.0308 0.0352 Inf 1 -3.041 0.0071
P value adjustment: fdr method for 3 tests
Tests are performed on the log odds ratio scale The output shows the results of the multiple comparison (pairwise) testing.
The values in the p.value column tell you the results of the hypothesis test comparing the coefficients between the two levels of your categorical predictor.
Note the odds.ratio column here. Comparing effects of Cat1 on the probability of Resp1.b is done by comparing odds of Resp1.b being 1 under each level of Cat1.
For example:
there is no evidence that the probability of Resp1.b being 1 is different when Cat 1 = SpA vs. when Cat1 = SpB (P = 0.1107).
there is little evidence that the probability of Resp1.b being 1 is different when Cat 1 = SpA than when Cat1 = SpC (P = 0.0413). The odds ratio is 0.2.
there is strong evidence that the probability of Resp1.b being 1 is different when Cat 1 = SpB than when Cat1 = SpC (P = 0.0071). The odds ratio is 0.0308.
This is a ratio of the odds of that Resp1.b is 1 when Cat 1 = SpA vs. the odds of that Resp1.b is 1 when Cat 1 = SpC:
For example, the odds that Resp1.b is 1 when Cat 1 = SpA are:
pCat1.SpA <- summary(emmOut)[1,2]
oddsCat1.SpA <- pCat1.SpA /(1-pCat1.SpA)
oddsCat1.SpA[1] 0.5and the odds that Resp1.b is 1 when Cat 1 = SpC are:
pCat1.SpC <- summary(emmOut)[3,2]
oddsCat1.SpC <- pCat1.SpC /(1-pCat1.SpC)
oddsCat1.SpC[1] 2.5and the odds ratio StA/StB is
oddsCat1.SpA/oddsCat1.SpC[1] 0.2Note that you can also get the results from the pairwise testing visually by handing the output of pairs() to plot().
Recall that Example 2 contains two predictors and both are categorical:
Resp2.b ~ Cat2 + Cat3 + Cat2:Cat3 + 1
The best-specified model (now with binomial error distribution assumption) is:
bestMod2.b
Call: glm(formula = Resp2.b ~ Cat2 + Cat3 + Cat2:Cat3 + 1, family = binomial(link = "logit"),
data = myDat2.b)
Coefficients:
(Intercept) Cat2StB Cat2StC Cat3Treat
-0.15415 0.23639 -0.08352 -0.10368
Cat2StB:Cat3Treat Cat2StC:Cat3Treat
-0.29222 1.99958
Degrees of Freedom: 499 Total (i.e. Null); 494 Residual
Null Deviance: 692.9
Residual Deviance: 649.7 AIC: 661.7that was fit to data in myDat2.b:
str(myDat2.b)'data.frame': 500 obs. of 3 variables:
$ Cat2 : Factor w/ 3 levels "StA","StB","StC": 2 2 3 2 2 3 2 1 1 1 ...
$ Cat3 : Factor w/ 2 levels "Control","Treat": 1 1 1 2 1 2 1 2 2 1 ...
$ Resp2.b: int 0 0 0 1 1 1 0 0 0 1 ...and the dredge() table you used to pick your bestMod2.b is in
dredgeOut2.bGlobal model call: glm(formula = Resp2.b ~ Cat2 + Cat3 + Cat2:Cat3 + 1, family = binomial(link = "logit"),
data = myDat2.b)
---
Model selection table
(Int) Ct2 Ct3 Ct2:Ct3 R^2 df logLik AICc delta weight
8 -0.15420 + + + 0.082720 6 -324.843 661.9 0.00 1
4 -0.38060 + + 0.031870 4 -338.331 684.7 22.89 0
2 -0.20190 + 0.023290 3 -340.538 687.1 25.27 0
3 -0.11690 + 0.007163 2 -344.632 693.3 31.43 0
1 0.04801 0.000000 1 -346.430 694.9 33.01 0
Models ranked by AICc(x) And here’s a visualization of the model effects:
#### i) choosing the values of your predictors at which to make a prediction
# Set up your predictors for the visualized fit
forCat2<-unique(myDat2.b$Cat2) # every level of your categorical predictor
forCat3<-unique(myDat2.b$Cat3) # every level of your categorical predictor
# create a data frame with your predictors
forVis<-expand.grid(Cat2 = forCat2, Cat3 = forCat3) # expand.grid() function makes sure you have all combinations of predictors
#### ii) using `predict()` to use your model to estimate your response variable at those values of your predictor
# Get your model fit estimates at each value of your predictors
modFit<-predict(bestMod2.b, # the model
newdata = forVis, # the predictor values
type = "link", # here, make the predictions on the link scale and we'll translate the back below
se.fit = TRUE) # include uncertainty estimate
ilink <- family(bestMod2.b)$linkinv # get the conversion (inverse link) from the model to translate back to the response scale
forVis$Fit<-ilink(modFit$fit) # add your fit to the data frame
forVis$Upper<-ilink(modFit$fit + 1.96*modFit$se.fit) # add your uncertainty to the data frame, 95% confidence intervals
forVis$Lower<-ilink(modFit$fit - 1.96*modFit$se.fit) # add your uncertainty to the data frame
#### iii) use the model estimates to plot your model fit
library(ggplot2) # load ggplot2 library
ggplot() + # start ggplot
geom_point(data = myDat2.b, # add observations to your plot
mapping = aes(x = Cat2, y = Resp2.b, col = Cat3),
position=position_jitterdodge(jitter.width=0.75, dodge.width=0.75)) + # control position of data points so they are easier to see on the plot
geom_errorbar(data = forVis, # add the uncertainty to your plot
mapping = aes(x = Cat2, y = Fit, ymin = Lower, ymax = Upper, col = Cat3),
position=position_dodge(width=0.75), # control position of data points so they are easier to see on the plot
size=1.2) + # control thickness of errorbar line
geom_point(data = forVis, # add the modelled fit to your plot
mapping = aes(x = Cat2, y = Fit, fill = Cat3),
shape = 22, # shape of point
size=3, # size of point
col = 'black', # outline colour for point
position=position_dodge(width=0.75)) + # control position of data points so they are easier to see on the plot
ylab("probability Resp2.b = 1")+ # change y-axis label
xlab("Cat2, (units)")+ # change x-axis label
labs(fill="Cat3, (units)", col="Cat3, (units)") + # change legend title
theme_bw() # change ggplot theme
What are your modelled effects (with uncertainty)?
As above, we can use the emmeans package to see the effects of the predictors on the scale of the response:
emmOut <- emmeans(object = bestMod2.b, # your model
specs = ~ Cat2 + Cat3 + Cat2:Cat3, # your predictors
type = "response") # report coefficients on the response scale
emmOut Cat2 Cat3 prob SE df asymp.LCL asymp.UCL
StA Control 0.462 0.0523 Inf 0.362 0.564
StB Control 0.521 0.0585 Inf 0.407 0.632
StC Control 0.441 0.0515 Inf 0.344 0.543
StA Treat 0.436 0.0561 Inf 0.331 0.547
StB Treat 0.422 0.0521 Inf 0.325 0.526
StC Treat 0.840 0.0423 Inf 0.739 0.907
Confidence level used: 0.95
Intervals are back-transformed from the logit scale So now you have a modelled value of your response for each level of your categorical predictors. For example:
When Cat2 is StA and Cat3 is Treat, the probability that Resp1.b is 1 is estimated to be 0.44 (95% confidence interval: 0.33 to 0.55 units).
When Cat2 is StB and Cat3 is Control, the probability that Resp1.b is 1 is estimated to be 0.52 (95% confidence interval: 0.41 to 0.63 units).
Finally, note that you can also get a quick plot of the effects by handing the emmeans() output to plot().
Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)
As with Example 1, you can also find out which combinations of predictor levels are leading to statistically different model predictions in Resp2.b:
forComp <- pairs(emmOut, # emmeans output
adjust = "fdr", # multiple comparison adjustment
simple = "Cat2") # contrasts within this categorical predictor
forCompCat3 = Control:
contrast odds.ratio SE df null z.ratio p.value
StA / StB 0.789 0.2490 Inf 1 -0.751 0.6791
StA / StC 1.087 0.3220 Inf 1 0.282 0.7781
StB / StC 1.377 0.4320 Inf 1 1.019 0.6791
Cat3 = Treat:
contrast odds.ratio SE df null z.ratio p.value
StA / StB 1.057 0.3300 Inf 1 0.179 0.8582
StA / StC 0.147 0.0573 Inf 1 -4.925 <0.0001
StB / StC 0.139 0.0530 Inf 1 -5.183 <0.0001
P value adjustment: fdr method for 3 tests
Tests are performed on the log odds ratio scale Based on a threshold P-value of 0.05, you can see that the effects at some combinations of your predictors are not statistically different from each other. For example, there is no evidence that the probability of Resp2.b being 1 when Cat3 = Control and Cat2 = StA (0.46)24 than when Cat3 = Control and Cat2 = StB (0.52) is different from each other (P = 0.68).
On the other hand, some combinations of your predictors are statistically different from each other. For example, there is strong evidence that the probability of Resp2.b being 1 when Cat3 = Treat and Cat2 = StA (0.46) is less than that when Cat3 = Treat and Cat2 = StC (0.84) is different from each other (P < 0.0001).
Note that you can also get the results from the pairwise testing visually by handing the output of pairs() to plot().
Recall that Example 3 contains one predictor and the predictor is numeric:
Resp3 ~ Num4 + 1
The best-specified model (now with a binomial error distribution assumption) is:
bestMod3.b
Call: glm(formula = Resp3.b ~ Num4 + 1, family = binomial(link = "logit"),
data = myDat3.b)
Coefficients:
(Intercept) Num4
-8.1519 0.1086
Degrees of Freedom: 299 Total (i.e. Null); 298 Residual
Null Deviance: 415.8
Residual Deviance: 199.1 AIC: 203.1that was fit to data in myDat3.b:
str(myDat3.b)'data.frame': 300 obs. of 2 variables:
$ Num4 : num 57.7 53.2 79.7 35.1 72.2 ...
$ Resp3.b: int 1 0 1 0 0 1 1 1 1 0 ...and the dredge() table you used to pick your bestMod3.b is
dredgeOut3.bGlobal model call: glm(formula = Resp3.b ~ Num4 + 1, family = binomial(link = "logit"),
data = myDat3.b)
---
Model selection table
(Intrc) Num4 R^2 df logLik AICc delta weight
2 -8.15200 0.1086 0.5144 2 -99.533 203.1 0.00 1
1 0.04001 0.0000 1 -207.884 417.8 214.68 0
Models ranked by AICc(x) And here’s a visualization of the model effects:
#### i) choosing the values of your predictors at which to make a prediction
# Set up your predictors for the visualized fit
forNum4<-seq(from = min(myDat3.b$Num4), to = max(myDat3.b$Num4), by = 1)# a sequence of numbers in your numeric predictor range
# create a data frame with your predictors
forVis<-expand.grid(Num4 = forNum4) # expand.grid() function makes sure you have all combinations of predictors.
#### ii) using `predict()` to use your model to estimate your response variable at those values of your predictor
# Get your model fit estimates at each value of your predictors
modFit<-predict(bestMod3.b, # the model
newdata = forVis, # the predictor values
type = "link", # here, make the predictions on the link scale and we'll translate the back below
se.fit = TRUE) # include uncertainty estimate
ilink <- family(bestMod3.b)$linkinv # get the conversion (inverse link) from the model to translate back to the response scale
forVis$Fit<-ilink(modFit$fit) # add your fit to the data frame
forVis$Upper<-ilink(modFit$fit + 1.96*modFit$se.fit) # add your uncertainty to the data frame, 95% confidence intervals
forVis$Lower<-ilink(modFit$fit - 1.96*modFit$se.fit) # add your uncertainty to the data frame
#### iii) use the model estimates to plot your model fit
library(ggplot2) # load ggplot2 library
ggplot() + # start ggplot
geom_point(data = myDat3.b, # add observations to your plot
mapping = aes(x = Num4, y = Resp3.b)) + # control position of data points so they are easier to see on the plot
geom_line(data = forVis, # add the modelled fit to your plot
mapping = aes(x = Num4, y = Fit),
linewidth = 1.2) + # control thickness of line
geom_line(data = forVis, # add uncertainty to your plot (upper line)
mapping = aes(x = Num4, y = Upper),
linewidth = 0.8, # control thickness of line
linetype = 2) + # control style of line
geom_line(data = forVis, # add uncertainty to your plot (lower line)
mapping = aes(x = Num4, y = Lower),
linetype = 2) + # control style of line
ylab("probability Resp2.b = 1")+ # change y-axis label
xlab("Num4, (units)") + # change x-axis label
theme_bw() # change ggplot theme
The model fit line may seem a little strange as the data are binary (can only be 0 or 1) but the model fit line goes through the space in between 0 and 1. This model fit line tells you how the probability of Resp3.b being 1 changes as Num4 changes.
What are your modelled effects (with uncertainty)?
trendsOut <- emtrends(bestMod3.b,
specs = ~ Num4 + 1, # your predictors
var = "Num4") # your predictor of interest
trendsOut Num4 Num4.trend SE df asymp.LCL asymp.UCL
75.8 0.109 0.0116 Inf 0.0859 0.131
Confidence level used: 0.95 Note that this is the same estimate as you find in the summary() output of your model object:
coef(summary(bestMod3.b)) Estimate Std. Error z value Pr(>|z|)
(Intercept) -8.1518996 0.89160247 -9.142976 6.075711e-20
Num4 0.1085946 0.01160302 9.359159 8.037377e-21This is because emtrends() returns coefficients on the link scale.
To report the coefficients on a response scale, you need to consider your link function. Since you have a binomial error distribution assumption, you can convert the estimate made by emtrends() to the response scale with:
trendsOut.df <- data.frame(trendsOut) # convert trendsOut to a dataframe
OR <- exp(trendsOut.df$Num4.trend) # get the coefficient on the response scale
OR[1] 1.11471You can also use the same conversion on the confidence limits of the modelled effect, for the upper:
OR.up <- exp(trendsOut.df$asymp.UCL) # get the upper confidence interval on the response scale
OR.up[1] 1.140351and lower 95% confidence level:
OR.low <- exp(trendsOut.df$asymp.LCL) # get the lower confidence interval on the response scale
OR.low[1] 1.089646This coefficient - called the odds ratio (?@sec-startBinomial) - tells you that for a unit change in Num4, the odds that Resp3.b is 1 increases by 1.115 times (95% confidence interval is 1.09 to 1.14) which is an increase in the odds of Resp3.b being 1 of 11.5%.
You can get a better sense of the odds ratio (vs. odds vs. probabilities) by estimating it directly from the modelled probability that Resp3.b is 1:
Let’s find the probability of Resp3.b = 1 when Num4 is 60:
pNum4.60 <- predict(bestMod3.b, # model
newdata = data.frame(Num4 = 60), # predictor values for the prediction
type = "response", # give prediction on the response scale
se.fit = TRUE) # include uncertainty
pNum4.60 <- pNum4.60$fit
pNum4.60 1
0.1629792 So there is a 16.3% probability of Resp3.b being 1 when Num4 is 60.
The odds associated with Resp3.b being 1 when Num4 is 60 is given by the probability of Num4 is divided by the probability of absence:
oddsNum4.60 <- pNum4.60/(1-pNum4.60)
oddsNum4.60 1
0.1947134 Now let’s find the probability of Resp3.b being 1 when Num4 is 61 (one unit more):
pNum4.61 <- predict(bestMod3.b, # model
newdata = data.frame(Num4 = 61), # predictor values for the prediction
type = "response", # give prediction on the response scale
se.fit = TRUE) # include uncertainty
pNum4.61 <- pNum4.61$fit
pNum4.61 1
0.1783404 So there is a 17.8% probability of Resp3.b being 1 when Num4 is 61.
The odds associated with Resp3.b being 1 when Num4 is 61 is given by the probability of presence divided by the probability of absence:
oddsNum4.61 <- pNum4.61/(1-pNum4.61)
oddsNum4.61 1
0.217049 Now the odds ratio is found by comparing the change in odds when your predictor (Num4) changes one unit
oddsNum4.61/oddsNum4.60 1
1.11471 Note that this is identical to the calculation of the odds ratio from the coefficient above:
OR[1] 1.11471and you interpret this as a OR - 1 or 11.5% change in the odds of Resp3.b being 1 with a unit change of Num4.
Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)
This question only applies to categorical predictors of which their are none in Example 3.
Recall that Example 4 contains two predictors, one is categorical and one is numeric:
Resp4.b ~ Cat5 + Num6 + Cat5:Num6 + 1
The best-specified model (now with a binomial error distribution assumption) is:
bestMod4.b
Call: glm(formula = Resp4.b ~ Cat5 + Num6 + Cat5:Num6 + 1, family = binomial(link = "logit"),
data = myDat4.b)
Coefficients:
(Intercept) Cat5StB Cat5StC Num6 Cat5StB:Num6
-28.55811 19.23667 -10.40094 0.15899 -0.09544
Cat5StC:Num6
0.03623
Degrees of Freedom: 499 Total (i.e. Null); 494 Residual
Null Deviance: 639.7
Residual Deviance: 169 AIC: 181that was fit to data in myDat4.b:
str(myDat4.b)'data.frame': 500 obs. of 3 variables:
$ Cat5 : Factor w/ 3 levels "StA","StB","StC": 2 1 1 1 2 1 2 1 3 2 ...
$ Num6 : num 206 191 312 223 210 ...
$ Resp4.b: int 1 1 1 1 1 1 1 0 1 1 ...and the dredge() table you used to pick your bestMod4.b is
dredgeOut4.bGlobal model call: glm(formula = Resp4.b ~ Cat5 + Num6 + Cat5:Num6 + 1, family = binomial(link = "logit"),
data = myDat4.b)
---
Model selection table
(Int) Ct5 Nm6 Ct5:Nm6 R^2 df logLik AICc delta weight
8 -28.5600 + 0.15900 + 0.60990 6 -84.513 181.2 0.00 1
4 -19.4300 + 0.10820 0.59390 4 -94.548 197.2 15.98 0
3 -12.4100 0.07056 0.49520 2 -148.970 302.0 120.77 0
2 0.5276 + 0.07603 3 -300.081 606.2 425.01 0
1 0.6722 0.00000 1 -319.850 641.7 460.51 0
Models ranked by AICc(x) And here’s a visualization of the model effects:
#### i) choosing the values of your predictors at which to make a prediction
# Set up your predictors for the visualized fit
forCat5<-unique(myDat4.b$Cat5) # all levels of your categorical predictor
forNum6<-seq(from = min(myDat4.b$Num6), to = max(myDat4.b$Num6), by = 1)# a sequence of numbers in your numeric predictor range
# create a data frame with your predictors
forVis<-expand.grid(Cat5 = forCat5, Num6 = forNum6) # expand.grid() function makes sure you have all combinations of predictors.
#### ii) using `predict()` to use your model to estimate your response variable at those values of your predictor
# Get your model fit estimates at each value of your predictors
modFit<-predict(bestMod4.b, # the model
newdata = forVis, # the predictor values
type = "link", # here, make the predictions on the link scale and we'll translate the back below
se.fit = TRUE) # include uncertainty estimate
ilink <- family(bestMod4.b)$linkinv # get the conversion (inverse link) from the model to translate back to the response scale
forVis$Fit<-ilink(modFit$fit) # add your fit to the data frame
forVis$Upper<-ilink(modFit$fit + 1.96*modFit$se.fit) # add your uncertainty to the data frame, 95% confidence intervals
forVis$Lower<-ilink(modFit$fit - 1.96*modFit$se.fit) # add your uncertainty to the data frame
#### iii) use the model estimates to plot your model fit
library(ggplot2) # load ggplot2 library
ggplot() + # start ggplot
geom_point(data = myDat4.b, # add observations to your plot
mapping = aes(x = Num6, y = Resp4.b, col = Cat5)) + # control position of data points so they are easier to see on the plot
geom_line(data = forVis, # add the modelled fit to your plot
mapping = aes(x = Num6, y = Fit, col = Cat5),
linewidth = 1.2) + # control thickness of line
geom_line(data = forVis, # add uncertainty to your plot (upper line)
mapping = aes(x = Num6, y = Upper, col = Cat5),
linewidth = 0.4, # control thickness of line
linetype = 2) + # control style of line
geom_line(data = forVis, # add uncertainty to your plot (lower line)
mapping = aes(x = Num6, y = Lower, col = Cat5),
linewidth = 0.4, # control thickness of line
linetype = 2) + # control style of line
ylab("probability Resp2.b = 1")+ # change y-axis label
xlab("Num6, (units)") + # change x-axis label
labs(fill="Cat5, (units)", col="Cat5, (units)") + # change legend title
theme_bw() # change ggplot theme
What are your modelled effects (with uncertainty)?
The process for estimating effects of categorical vs. numeric predictors is different.
For categorical predictors (Cat5), you can use emmeans() to give you the effects on the response scale directly:
emmOut <- emmeans(object = bestMod4.b, # your model
specs = ~ Cat5 + Num6 + Cat5:Num6 + 1, # your predictors
type = "response") # report coefficients on the response scale
emmOut Cat5 Num6 prob SE df asymp.LCL asymp.UCL
StA 198 0.951 0.0346 Inf 0.819 0.988
StB 198 0.964 0.0178 Inf 0.907 0.986
StC 198 0.438 0.0941 Inf 0.269 0.622
Confidence level used: 0.95
Intervals are back-transformed from the logit scale This output tells you that:
When Cat5 is StA and Num6 is 227 units25, there is a 95% probability that Resp4.b will be 1 (95% confidence interval: 0.82 to 0.99%)
When Cat5 is StB and Num6 is 227 units26, there is a 96% probability that Resp4.b will be 1 (95% confidence interval: 0.91 to 0.99 units).
When Cat5 is StB and Num6 is 227 units27, there is a 44% probability that Resp4.b will be 1 (95% confidence interval: 0.27 to 0.62 units).
Note that emmeans() sets your numeric predictor (Num6) to the mean value of Num6 (227 units). You can also control this with the at = argument in the emmeans() function:
emmOut <- emmeans(object = bestMod4.b, # your model
specs = ~ Cat5 + Num6 + Cat5:Num6 + 1, # your predictors
type = "response", # report coefficients on the response scale
at = list(Num6 = 150)) # control the value of your numeric predictor at which to make the coefficient estimates
emmOut Cat5 Num6 prob SE df asymp.LCL asymp.UCL
StA 150 8.93e-03 0.009230 Inf 1.17e-03 0.06500
StB 150 5.53e-01 0.078500 Inf 3.99e-01 0.69695
StC 150 6.27e-05 0.000142 Inf 7.00e-07 0.00524
Confidence level used: 0.95
Intervals are back-transformed from the logit scale For numeric predictors (Num6), you need to use the emtrends() function and then convert the coefficients to the response scale:
trendsOut <- emtrends(bestMod4.b,
specs = ~ Cat5 + Num6 + Cat5:Num6 + 1 , # your predictors
var = "Num6") # your predictor of interest
trendsOut Cat5 Num6 Num6.trend SE df asymp.LCL asymp.UCL
StA 198 0.1590 0.0330 Inf 0.0942 0.2237
StB 198 0.0635 0.0114 Inf 0.0413 0.0858
StC 198 0.1952 0.0444 Inf 0.1083 0.2821
Confidence level used: 0.95 Since you have a binomial error distribution assumption, you can convert the estimate made by emtrends() on to the response scale with the exp() function:
trendsOut.df <- data.frame(trendsOut) # convert trendsOut to a dataframe
OR <- exp(trendsOut.df$Num6.trend) # get the coefficient on the response scale
OR[1] 1.172328 1.065611 1.215576You can also use the same conversion on the confidence limits of the modelled effect, for the upper:
OR.up <- exp(trendsOut.df$asymp.UCL) # get the upper confidence interval on the response scale
OR.up[1] 1.250737 1.089631 1.325976and lower 95% confidence level:
OR.low <- exp(trendsOut.df$asymp.LCL) # get the lower confidence interval on the response scale
OR.low[1] 1.098834 1.042120 1.114368Let’s organize these numbers so we can read the effects easier:
ORAll <- data.frame(Cat5 = trendsOut.df$Cat5, # include info on the Cat5 level
OR = exp(trendsOut.df$Num6.trend), # the modelled effect as a rate ratio
OR.up = exp(trendsOut.df$asymp.UCL), # upper 95% confidence level of the modelled effect
OR.down = exp(trendsOut.df$asymp.LCL)) # lower 95% confidence level of the modelled effect
ORAll Cat5 OR OR.up OR.down
1 StA 1.172328 1.250737 1.098834
2 StB 1.065611 1.089631 1.042120
3 StC 1.215576 1.325976 1.114368This tells you that for a unit change in Num6, the odds that Resp4.b is 1:
increases by 17.2% with a unit change in Num6 when Cat5 is StA (95% confidence interval: 9.9 to 25.1).
increases by 6.6% with a unit change in Num6 when Cat5 is StB (95% confidence interval: 4.2 to 9).
increases by 21.6% with a unit change in Num6 when Cat5 is StC (95% confidence interval: 11.4 to 32.6).
You can test for evidence that the effects of Num6 on Resp4.b are significant (different than zero) with:
test(trendsOut) # get test if coefficient is different than zero. Cat5 Num6 Num6.trend SE df z.ratio p.value
StA 198 0.1590 0.0330 Inf 4.813 <0.0001
StB 198 0.0635 0.0114 Inf 5.588 <0.0001
StC 198 0.1952 0.0444 Inf 4.401 <0.0001Note that these tests are done on the link scale but can be used for reporting effects on the response scale.
From the results, you can see that the effects of Num6 on Resp4.b at all levels of Cat5 are significant (P < 0.0001).
Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)
You can also find out which effects of Num6 on Resp4.b are different from one another using pairs() on the emtrends() output:
forCompOR <- pairs(trendsOut, # emmeans object
adjust = "fdr", # multiple comparison adjustment
simple = "Cat5") # contrasts within this categorical predictor
forCompORNum6 = 198:
contrast estimate SE df z.ratio p.value
StA - StB 0.0954 0.0349 Inf 2.732 0.0094
StA - StC -0.0362 0.0553 Inf -0.655 0.5124
StB - StC -0.1317 0.0458 Inf -2.876 0.0094
P value adjustment: fdr method for 3 tests The table above tells you that:
There is evidence that the effect of Num6 on Resp4.b is different when Cat5 is StA vs. StB (P = 0.009). The effect of Num6 on Resp4.b is exp(0.0954) = 1.1 times as big when Cat5 is StA vs. StB (i.e. the effect of Num6 on Resp4.b is 10% more when Cat5 is StA vs. StB)
There is no evidence the effect of Num6 on Resp4.b differs when Cat5 is StA vs. StC (P = 0.512).
There is evidence that the effect of Num6 on Resp4.b is different when Cat5 is StB vs. StC (P = 0.009). The effect of Num6 on Resp4.b is exp(-0.1317) = 0.88 times as big when Cat5 is StB vs. StC (i.e. the effect of Num6 on Resp4.b is 12% less when Cat5 is StB vs. StC)
Only if all the slopes28 were similar, you would want to test if the levels of your categorical predictor (Cat5) have significantly different coefficients (intercepts) from each other with pairs() on the output from emmeans() (the emmOut object).
For categorical predictors, this is the same as the process in the section on giving examples of your modelled effects↩︎
“post hoc” is a Latin phrase meaning “after the event”↩︎
note that the p-values reported in the coefficient table are the result of a test of whether the coefficient associated with Sp2 is different than that of Sp1 (p = 0.001), and if the effect of Sp3 is different than that of Sp1 (p = 0.017)), but a comparison of effects of Sp2 vs. Sp3 is missing. We will discuss this more further along in these notes.↩︎
emmeans stands for “estimated marginal means”↩︎
another option is the marginaleffects package, where marginal effects represent the average effect - the average difference in your predicted response across a change in your categorical predictor level or a unit change in your numeric predictor(Arel-Bundock, Greifer, and Heiss 2024)↩︎
note that the summary() output becomes easier to read if you force the starting model not to have an intercept with: startMod<-glm(formula = Resp1 ~ Cat1 - 1, data = myDat, family = gaussian(link="identity")). This is fine to do here where you only have categorical predictors but causes problems when you start having numeric predictors in your model as well. So we won’t be practicing this in class and instead will use the very flexible and helpful emmeans package to help us report coefficients from our statistical models.↩︎
when the P-value is very low, R reports is as simple < 0.0001↩︎
note that now you get to assess the difference between coefficients for Sp2 and Sp3 which was missing in the summary() output above↩︎
The emmeans package is very flexible and has a lot of options as to how to make these corrections depending on your needs. Plenty more information is available here↩︎
* in units of Resp3 per units of Num4↩︎
units of Resp4↩︎
units will be the units of Resp4 divided by those of Num6↩︎
in units of Resp4 per units of Num6↩︎
in units of Resp4 per units of Num6↩︎
in units of Resp4 per units of Num6↩︎
in units of Resp4 per units of Num6↩︎
also called the incidence rate ratio, incidence density ratio or relative risk↩︎
a reminder that the number e that the function exp() uses as its base is Euler’s number. It (along with the natural logarithm, log()) is used to model exponential growth or decay, and can be used here when you want models where the response can not be below zero.↩︎
this is what R uses when you say link = "log"↩︎
from the emmeans package↩︎
e.g. because the model validation is better for your data↩︎
Jargon alert! Note that a GLM using the logit link function is also called logistic regression↩︎
from the emmeans package↩︎
taken from the emmOut output above↩︎
the mean of the range of Num6↩︎
the mean of the range of Num6↩︎
the mean of the range of Num6↩︎
i.e. effects associated with the numeric predictor↩︎