Introduction to Generalized Linear Models in the Psychological Sciences: R code
2024-05-01
Overview
This code is designed to walk through how to estimate, interpret, and plot results from generalized linear models (GLMs) using a marginal effects approach.
Setup Code
This chunk is for getting our data loaded in and setting up a few helper functions.
#Keep scientific notation to 3 places, digits to 3 places.
options(scipen=3,digits=3)
require(tidyverse)
require(RCurl)
require(marginaleffects)
require(knitr)
require(lm.beta)
require(packrat)
require(questionr)
require(MASS)
df<-read.table("ex_clean_QF.txt",header = T)%>%
mutate(aud=case_when(negpysum<=6~0,negpysum>6~1),
sex=as.factor(recode(gen,"1"="male","2"="female")),
female=recode(sex,"male"=0,"female"=1),
qf=QF,
ss=scale(ss)[,1])%>%
dplyr::select(-c(gen,purg, QF,pre,negpysum,sex))#This creates a quick descriptives table of the data.
kable(psych::describe(df)%>%
as.data.frame()%>%
dplyr::select(n, mean, sd,min,max))| n | mean | sd | min | max | |
|---|---|---|---|---|---|
| ss | 427 | 0.000 | 1.000 | -2.71 | 1.75 |
| aud | 423 | 0.416 | 0.493 | 0.00 | 1.00 |
| female | 427 | 0.520 | 0.500 | 0.00 | 1.00 |
| qf | 398 | 19.771 | 19.276 | 0.00 | 86.00 |
A few notes on the variables:
qf refers to alcohol quantity/frequency. It is count
distributed, meaning it has a lower bound of 0 and assumes integer
values. I will (hopefully, if time) use it as an example for count data
cases.
aud refers to alcohol use disorder (AUD) diagnosis. This
assumes only 2 values (0=“No”,1=“Yes). This is a made-up, dichotomized
variable derived from alcohol consequences.
female is a dummy-coded biological sex assigned at birth
variable (0=male, 1=female).
ss refers to sensation seeking, our primary continuous
variable predictor. This is a standardized variable (i.e., Mean = 0, SD
= 1).
Part 1: Model Estimation
To start, let’s estimate a logistic model regressing AUD on sex and sensation seeking.
logistic<-glm(aud~female + ss,data=df, family="binomial")
summary(logistic)##
## Call:
## glm(formula = aud ~ female + ss, family = "binomial", data = df)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.274 0.148 -1.86 0.063 .
## female -0.170 0.207 -0.82 0.412
## ss 0.445 0.110 4.05 0.000052 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 574.43 on 422 degrees of freedom
## Residual deviance: 553.57 on 420 degrees of freedom
## (4 observations deleted due to missingness)
## AIC: 559.6
##
## Number of Fisher Scoring iterations: 4Above, the value 0.445 refers to the effect of a 1-unit increase in
sensation seeking on the log-odds of AUD, above and beyond the effect of
female.
Note that these results are the transformed scale of our model - that is, the model used for estimating our outcome on the log-odds scale, as shown below:
\[log(\frac{\hat{\pi}}{1-\hat{\pi}}) = -0.274 + -0.17female + 0.445ss \]
We can generate a plot of these results using the
plot_predictions() function from
marginaleffects. For now, I will plot these effects from
the range of -3 to +3 for my sensation seeking variable,
and I will hold the female covariate at its sample mean
(0.52).
predplot.logistic<-plot_predictions(logistic,condition=list("ss"=-3:3,
"female"=0.52),
type="link")+
xlab("Sensation Seeking") +
ylab("log-odds of AUD") +
scale_fill_manual(values="black")+
scale_color_manual(values="black")+
theme(legend.position="none") +
theme_bw() +
theme(legend.position="none")
predplot.logistic
This plot shows the slope of sensation seeking on the log-odds of AUD. Because this is estimated on the linear (i.e., link) scale, this is simply the b value provided in our summary table.
However, a fundamental issue is that we may not be interested in log-odds, but rather odds or probabilities.
We can obtain these quantities by performing some transformations to our logistic regression model.
To start, we can exponentiate both sides of our equation to yield “odds” instead of “log-odds”:
\[\frac{\hat{\pi}}{1-\hat{\pi}} = exp(-0.274 + -0.17female + 0.445ss)\]
\[\frac{\hat{\pi}}{1-\hat{\pi}} = 0.76 + −0.844female + 1.56ss\]
These exponentiated values are the values typically reported in the results sections of manuscripts - known as odds ratios (ORs). ORs describe how the odds of an outcome change for a one-unit change in a predictor variable.
For instance, for our sensation seeking variable, the exponentiated
value for our coefficient is exp(0.445) = 1.56. This indicates that for
every one-unit increase in sensation seeking, the odds of AUD increases
by about 56% (95% CI = [1.262, 1.945]. The odds ratio
function from questionr provides a quick and easy summary
of this:
odds.ratio(logistic)## Waiting for profiling to be done...## OR 2.5 % 97.5 % p
## (Intercept) 0.760 0.567 1.01 0.063 .
## female 0.844 0.562 1.27 0.412
## ss 1.561 1.262 1.94 0.000052 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Although marginaleffects doesn’t give us a
straightforward way of plotting odds, we can do it by hand below:
dg<-datagrid(ss= seq(-3,3,.5),
female=0.52,
model=logistic)
orpreds<-predictions(model=logistic,
newdata=dg,
type="link")
orpreds$OR<-exp(orpreds$estimate)
orpreds$OR_low<-exp(orpreds$conf.low)
orpreds$OR_high<-exp(orpreds$conf.high)
ggplot(data=orpreds,aes(x=ss,y=OR))+
geom_line() +
ylab("Odds of AUD") +
geom_ribbon(aes(ymin=orpreds$OR_low,ymax=orpreds$OR_high),alpha=.1) +
theme_bw()
Above, we see our first major implication of GLMs: the relation between our predictors and our outcome is no longer linear. Namely, we see in this plot that the effect of sensation seeking on the odds of AUD seems to accelerate at higher levels of sensation seeking. This is the result of now having an exponential relation between sensation seeking and AUD odds, rather than a linear relation with AUD log-odds. We’ll come back to the implications of this later.
While most researchers stop with odds ratios when describing results, I would strongly advise going one step further to describe results in the metric of probabilities. This is primarily because ORs are not comparable to other study findings. Moreover, they are often conflated with the (more informative) metric of relative risks (or risk ratios), which are the ratios of probabilities, not odds.
We can do this by further transforming our model as follows:
\[\hat{\pi}=\frac{e^{-0.274 + -0.17female + 0.445ss}}{1+e^{-0.274 + -0.17female + 0.445ss}}\]
This is a full inversion of our logistic link function that describes the quantity we most care about: the probability of AUD.
I know this formula looks a little hairy, but the
marginaleffects package provides some useful tools for
putting this in lay terms.
For instance, we can generate predicted probabilities in the same
manner as above by changing type="link" to
type="response" in our plot_predictions
function:
predplot.probs<-plot_predictions(logistic,condition=list("ss"=seq(-3,3,.1),
"female"=0.52),
type="response")+
xlab("Sensation Seeking") +
ylab("AUD Probability") +
ylim(0,1)+
scale_fill_manual(values="black")+
scale_color_manual(values="black")+
theme(legend.position="none") +
theme_bw() +
theme(legend.position="none")
predplot.probs
Looking under the hood, the values on these plots are simply plugging in values of sensation seeking to the probability equation above to generate predicted probabilities at these sensation seeking levels. Note again that, for now, we are holding “female” constant at 0.52.
Note further that on the probability scale, we see more of an S-shape relationship:
Although beyond the scope of this workshop, the standard errors of these quantities are typically generated using the delta method (https://marginaleffects.com/vignettes/uncertainty.html). In short, the delta method can derive an asymptotic variance for any function of asymptotically normal variables by taking their derivatives.
To make this concrete, let’s take a look at these predicted values:
dg<-datagrid(ss= -3:3,
female=0.52,
model=logistic)
glmpreds<-predictions(model=logistic,
newdata=dg,
type="response") #this is what indicates that we want the "natural response" scale and not the "transformed"/logit scale.
kable(glmpreds)| rowid | estimate | std.error | statistic | p.value | s.value | conf.low | conf.high | aud | ss | female |
|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 0.155 | 0.046 | 3.33 | 0.001 | 10.2 | 0.064 | 0.246 | 0 | -3 | 0.52 |
| 2 | 0.222 | 0.043 | 5.11 | 0.000 | 21.6 | 0.137 | 0.307 | 0 | -2 | 0.52 |
| 3 | 0.308 | 0.034 | 9.19 | 0.000 | 64.5 | 0.243 | 0.374 | 0 | -1 | 0.52 |
| 4 | 0.410 | 0.025 | 16.68 | 0.000 | 205.1 | 0.362 | 0.459 | 0 | 0 | 0.52 |
| 5 | 0.521 | 0.035 | 14.70 | 0.000 | 160.0 | 0.451 | 0.590 | 0 | 1 | 0.52 |
| 6 | 0.629 | 0.054 | 11.58 | 0.000 | 100.5 | 0.522 | 0.735 | 0 | 2 | 0.52 |
| 7 | 0.726 | 0.067 | 10.87 | 0.000 | 89.0 | 0.595 | 0.856 | 0 | 3 | 0.52 |
Above, we can see that generally, as sensation seeking increases, so does the probability of AUD. For instance, for someone with a mean level of sensation seeking (ss=0), the estimated probability of AUD is 0.41. When sensation seeking is 1 SD above its mean, our model predicts that one is more likely than not to have AUD (0.52 probability).
Tying it together
Let’s take another look at our predicted probability values:
glmpreds%>%dplyr::select(ss,estimate)%>%kable()| ss | estimate |
|---|---|
| -3 | 0.155 |
| -2 | 0.222 |
| -1 | 0.308 |
| 0 | 0.410 |
| 1 | 0.521 |
| 2 | 0.629 |
| 3 | 0.726 |
Note that each row represents a value at different levels of sensation seeking.
How could I compute odds from these values?
Since odds are the ratio of the probability of an occurance compared to not occuring, we could compute the odds using the formula we encountered earlier:
\[\frac{\hat{\pi}}{1-\hat{\pi}}\] In R code, we can do this as follows:
glmpreds$odds<-glmpreds$estimate/(1-glmpreds$estimate)
glmpreds%>%dplyr::select(ss,estimate, odds)%>%kable()| ss | estimate | odds |
|---|---|---|
| -3 | 0.155 | 0.183 |
| -2 | 0.222 | 0.286 |
| -1 | 0.308 | 0.446 |
| 0 | 0.410 | 0.696 |
| 1 | 0.521 | 1.086 |
| 2 | 0.629 | 1.695 |
| 3 | 0.726 | 2.645 |
So, how do we now understand an odds ratio? We can simply take the ratios of these odds. E.g., comparing odds of AUD at levels of sensation seeking 1 unit apart:
#Odds at ss=-3 divided by ss=-2:
glmpreds$odds[2]/glmpreds$odds[1]## [1] 1.56#Odds at ss=-2 dividied by ss=-1:
glmpreds$odds[3]/glmpreds$odds[2]## [1] 1.56Note that comparing any two of these adjacent values yields the value 1.56 - the odds ratio value we observed earlier when we exponentiated our log-odds coefficient.
odds.ratio(logistic)## Waiting for profiling to be done...## OR 2.5 % 97.5 % p
## (Intercept) 0.760 0.567 1.01 0.063 .
## female 0.844 0.562 1.27 0.412
## ss 1.561 1.262 1.94 0.000052 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Implications: How do odds and probabilities complicate standard linear modeling approaches?
There is a reason we’re having a GLM workshop rather than just a regression workshop. This reason is non-linearity inherent within GLMs that we do not see in standard Gaussian regression models. Because we are now interpreting results in nonlinear space, two additional aspects need to be addressed:
1. Because each predictor has a non-linear relation with the outcome, this relation cannot be described by coefficients produced by standard (linear) model output. One must instead make specific comparisons along a given predictor to more accurately describe effects.
2. Similarly, the non-linear design of GLMs means that effects for a given predictor on the natural scale are contingent on all other predictors in the model. One must therefore make careful decisions on how these effects are generated with respect to covariates. This also has significant implications on understanding interaction effects within these models.
Both of these considerations can be addressed using a marginal effects approach.
Marginal Effects: Definition
Marginal effects are rates-of-change in a given outcome as a function of change in a predictor, holding all else constant in some fashion. Depending on the nature of the question or the nature of the predictor variable, these are typically defined as using either partial derivatives for continuous variables or finite (or first) differences for continuous or categorical variables.
For example, let’s generate marginal effects for our logistic model
describing how sensation seeking predicts the probability of AUD using
the plot_slopes function in
marginaleffects:
dydxplot<-plot_slopes(logistic,variables="ss",condition=list("ss"=seq(-3,3,.1),"female"=0.52))+
xlab("sensation seeking") +
ylim(0,0.2) +
scale_fill_manual(values="black")+
scale_color_manual(values="black")+
theme_bw() +
theme(legend.position="none")
dydxplot
We see that our “slope” increases through about the mean of sensation seeking, then levels off thereafter. Note that we can see this pattern in our predicted probabilities plot from earlier
predplot.probs
These two plots are directly related: you cna think of the “slopes” in
the first plot as tangent lines along the S-curve in the second.
So, what do these slopes represent mathematically? Because sensation seeking is a continuous variable, these are partial derivative values, or:
\(\frac{\partial AUD}{\partial SS}\).
This represents the instantaneous rate-of-change of sensation seeking
on AUD probability across the range of sensation seeking values.
marginaleffects provides this using the slope
function.
Plots of marginal effect slopes provide a nice summary of where (and how much) along a continuous predictor do we see an association with the outcome variable. Importantly, for nonlinear design models, it is actually possible for a slope to be positive, non-significant (i.e., zero), or of opposite signs (positive and negative) depending on the level of the predictor (Ai & Norton, 2003). As such, evaluating derivatives across the range of a continuous variable provides a comprehensive assessment of a predictor’s effect on the natural scale of an outcome.
So what can we say pragmatically? In tandem with the
predictions plot we made earlier, we can make this even more concrete by
using the comparisons function to compute first (or finite)
differences. These are basically just descriptors of how the predicted
probability changes for a specified-unit change in the predictors. When
working with categorical predictors, we would want to go directly to
this approach.
To start, we can see how the probability changes when sensation
seeking is at the mean and female is 0. First, we create a data frame
for these hypotheticals. Then, we supply it to the
comparisons function.
hyp1<-data.frame(ss=0,
female=0)
comps1<-comparisons(logistic,newdata=hyp1)
kable(comps1)| rowid | term | contrast | estimate | std.error | statistic | p.value | s.value | conf.low | conf.high | predicted_lo | predicted_hi | predicted | ss | female | aud |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | female | 1 - 0 | -0.041 | 0.050 | -0.821 | 0.412 | 1.28 | -0.139 | 0.057 | 0.432 | 0.391 | 0.432 | 0 | 0 | 0 |
| 1 | ss | +1 | 0.111 | 0.027 | 4.069 | 0.000 | 14.37 | 0.057 | 0.164 | 0.432 | 0.543 | 0.432 | 0 | 0 | 0 |
Each row of this table is providing us with a summary of what happens
to AUD probability as each variable changes by one unit at the
hypothetical range we supplied to comparisons(). This
describes, in straightforward terms, what the relation is between
sensation seeking and AUD probability for as you move up by one unit
from the hypothetical sensation seeking mean.
Although this is very helpful information, we need to be very thoughtful about how these quantities are represented. For instance, our comparisons may not be identical if we took different hypothetical values for generating these comparisons. say, for example, we chose a unit increase in sensation seeking at +2 SDs from the mean:
hyp2<-data.frame(ss=2,
female=0)
comps2<-comparisons(logistic,newdata=hyp2)
kable(comps2)| rowid | term | contrast | estimate | std.error | statistic | p.value | s.value | conf.low | conf.high | predicted_lo | predicted_hi | predicted | ss | female | aud |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | female | 1 - 0 | -0.040 | 0.049 | -0.813 | 0.416 | 1.26 | -0.135 | 0.056 | 0.649 | 0.610 | 0.649 | 2 | 0 | 0 |
| 1 | ss | +1 | 0.094 | 0.015 | 6.190 | 0.000 | 30.63 | 0.064 | 0.123 | 0.649 | 0.743 | 0.649 | 2 | 0 | 0 |
Our two implications become evident here. First, as we showed above,
our finite difference estimate for sensation seeking is different at 2
SDs compared to the mean. Second, and perhaps more surprisingly, we
actually see that our estimate for female as different as
well (albeit only slightly in this example).
What this means is that unlike the linear case, describing the effect of a predictor on an outcome requires that covariates be held constant in some principled fashion.
In the early examples, I used the “hypothetical values” approach. For
instance, if I’m interested in the effect of sensation seeking
specifically, I may consider holding the female variable at
its mean of 0.52:
hyp3<-data.frame(ss=0,
female=0.52)
comps3<-comparisons(logistic,newdata=hyp3)## Warning: The `female` variable is treated as a categorical (factor) variable, but
## the original data is of class numeric. It is safer and faster to convert
## such variables to factor before fitting the model and calling a
## `marginaleffects` function.
##
## This warning appears once per session.kable(comps3)| rowid | term | contrast | estimate | std.error | statistic | p.value | s.value | conf.low | conf.high | predicted_lo | predicted_hi | predicted | ss | female | aud |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | female | 1 - 0 | -0.041 | 0.050 | -0.821 | 0.412 | 1.28 | -0.139 | 0.057 | 0.432 | 0.391 | 0.41 | 0 | 0.52 | 0 |
| 1 | ss | +1 | 0.110 | 0.027 | 4.026 | 0.000 | 14.11 | 0.057 | 0.164 | 0.410 | 0.521 | 0.41 | 0 | 0.52 | 0 |
Perhaps a better option is to compute the average marginal effect
(AME) using the avg_comparisons() function. In contrast to
hypothetical values supplied to generate comparisons, average
comparisons will compute a given comparison for each row of the dataset,
then average these values:
ame_comps<-avg_comparisons(logistic)
ame_comps##
## Term Contrast Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## female 1 - 0 -0.0395 0.0483 -0.818 0.413 1.3 -0.1341 0.0551
## ss +1 0.1049 0.0240 4.376 <0.001 16.3 0.0579 0.1519
##
## Columns: term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high
## Type: responseIn my opinion, AMEs are the best way of summarizing effects because they are most beholden to the original data while also providing a simpler, singular summary of a given effect.
Sample write-up
Note that derivatives, prediction plots, and the marginal effects are all candidates for what we would want to include in a results section. E.g., if we needed to describe the effect of sensation seeking on AUD probability succinctly using AMEs, we could state the following:
“We used average marginal effects to estimate the association between sensation seeking and AUD probability. For instance, an average one-unit increase in sensation seeking corresponded with a 0.105 change in the probability of AUD (95% CI = [0.058, 0.152].”
dydxplot+
ggtitle("Figure 1: marginal effects of sensation seeking on AUD probability")
predplot.probs
Extra Content (if time)
Interaction Analysis
As I mentioned above, an implication of odds, probability, and count outcomes is that any given effect must be computed with respect to all covariates in a model. This has significant implications for how interaction effects are computed and interpreted in these models.
In short, unlike the linear case, interaction effects for probabilities and counts are not equal to the product between variables on the transformed scale (McCabe et al., 2022). In fact, interaction between variables can be present even without a product term included in the model, and it is even possible that the actual interaction effect is in the opposite direction of the product term.
Let’s revisit our logistic model, but this time, we’ll include a
product between ss and female and take a look
at the predicted probabilities for males and females separately.
logistic.int<-glm(aud~ss+female+ss*female, data=df,family="binomial")
summary(logistic.int)##
## Call:
## glm(formula = aud ~ ss + female + ss * female, family = "binomial",
## data = df)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.3695 0.1608 -2.30 0.022 *
## ss 0.7414 0.1858 3.99 0.000066 ***
## female -0.0976 0.2141 -0.46 0.649
## ss:female -0.4813 0.2312 -2.08 0.037 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 574.43 on 422 degrees of freedom
## Residual deviance: 549.11 on 419 degrees of freedom
## (4 observations deleted due to missingness)
## AIC: 557.1
##
## Number of Fisher Scoring iterations: 4plot_predictions(logistic.int,condition=c("ss","female"),
type="response") +
theme_bw()
Above, we can see evidence of a crossover interaction, such that sensation seeking appears to be more strongly associated with AUD among males. Moreover, females may have higher probability of AUD at low levels of sensation seeking but lower probability at higher levels. Because we are evaluating interaction the probability scale, the product term does not quantify this interaction effect. We therefore need to use marginal effects to redefine it.
Because we are working with a continuous-by-categorical variable
interaction, the simplest way to do this is to quantify the AME for
sensation seeking among males, the AME for females, and then take the
difference between the two. We can employ the modglm
package to achieve this:
require(modglm)## Loading required package: modglmmodglm.logistic<-modglm(model=logistic.int, vars=c("ss","female"), data=df, type="fd", hyps="means")
modglm.logistic$aie## aie.est aie.se.delta aie.ll aie.ul
## 1 -0.101 0.0452 -0.19 -0.0125Count Example
Let’s assume the same predictors as before, but instead we’ll predict
the count of alcohol frequency (qf) rather than probability
of AUD.
Poisson models are used often to model count outcomes. They involve a log-link function, which we are now familiar with having dealt with odds from a logistic model. The model we will estimate is as follows:
\[log(\hat{qf}) = b_0 + female + ss \] Above, the natural scale of our outcome would be predicting the count of quantity/frequency. we can describe this model by simply exponentiating both sides of our equation above:
\[\hat{qf} = exp(b_0 + female + ss) \]
Below is an example of how to estimate this:
pois<-glm(qf ~ ss + female, data=df,family="poisson")
summary(pois)##
## Call:
## glm(formula = qf ~ ss + female, family = "poisson", data = df)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.0188 0.0165 182.6 < 2e-16 ***
## ss 0.2421 0.0124 19.6 < 2e-16 ***
## female -0.1085 0.0231 -4.7 0.0000026 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 8122.7 on 397 degrees of freedom
## Residual deviance: 7634.4 on 395 degrees of freedom
## (29 observations deleted due to missingness)
## AIC: 9089
##
## Number of Fisher Scoring iterations: 5 plot_predictions(pois,condition=list("ss"=seq(-3,3,.1),
"female"=0.52),
type="response") +
xlab("Sensation Seeking")+
ylab("Count of Alcohol Quantity/Frequency") +
scale_fill_manual(values="black")+
scale_color_manual(values="black")+
theme(legend.position="none") +
theme_bw() +
theme(legend.position="none")
hyp1<-data.frame(ss=0,
female=0)
comparisons(pois,
newdata=hyp1)##
## Term Contrast Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 % ss
## female 1 - 0 -2.10 0.449 -4.69 <0.001 18.5 -2.99 -1.22 0
## ss +1 5.61 0.308 18.24 <0.001 244.4 5.00 6.21 0
## female
## 0
## 0
##
## Columns: rowid, term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted, ss, female, qf
## Type: responseavg_comparisons(pois)##
## Term Contrast Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## female 1 - 0 -2.15 0.457 -4.7 <0.001 18.5 -3.04 -1.25
## ss +1 5.42 0.317 17.1 <0.001 215.0 4.80 6.04
##
## Columns: term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high
## Type: responseHelpful References
McCabe, C. J., Halvorson, M. A., King, K. M., Cao, X., & Kim, D. S. (2022). Interpreting Interaction Effects in Generalized Linear Models of Nonlinear Probabilities and Counts. Multivariate Behavioral Research, 57(2–3), 243–263. https://doi.org/10.1080/00273171.2020.1868966
Halvorson, M. A., McCabe, C. J., Kim, D. S., Cao, X., & King, K. M. (2022). Making sense of some odd ratios: A tutorial and improvements to present practices in reporting and visualizing quantities of interest for binary and count outcome models. Psychology of Addictive Behaviors, 36(3), 284–295. https://doi.org/10.1037/adb0000669
Coxe, S., West, S. G., & Aiken, L. S. (2013). Generalized linear models. The Oxford handbook of quantitative methods, 2, 26-51.
Howell‐Moroney, M. (2023). Inconvenient truths about logistic regression and the remedy of marginal effects. Public Administration Review, puar.13786. https://doi.org/10.1111/puar.13786