# Need help understanding odds ratio over time example

I'm trying to recreate a paper that compares the frequency and characteristics of emergency department visits that are repeats (ie; same patient had another emergency department visit for the same reason in the prior 12 months) vs single. In it, "Logistic regressions were used to calculate the probability of repeat ED visits compared to single ED visits occurring across time. Unadjusted odds ratios were presented; and no additional variables were controlled for." No other details are provided.

In the section where they present the odds ratio, they show a clustered bar graph showing the percentage of ed visits in each of 2018-2022 that were either single or repeat. There is a slight increase over the 5 years in the percentage each year that are repeats. They give the unadjusted odds ratio as 1.108 "which indicates that with each year, the odds of an individual having a repeat ED visit rather than a single ED visit increase by 11% (p<.001)."

Does that mean that to obtain unadjusted odds ratio I should run a logistic regression with my independent variable being the year of the ed visit and the dependent variable 0/1 for single/repeat? I'm having trouble understanding what the independent variable is in this situation.

I did try this (using sklearn.LogisticRegression) and got an odds ratio of 0.99, which doesn't make sense. Am I correct in thinking I should transform the years vector to 1-5 instead of 2018-2022? Am I way off base?

• A citation would be helpful. Your problem description is not clear enough. Commented Jul 16 at 18:44
• Thanks. I've added links to the paper. Commented Jul 16 at 20:34
• Welcome to Cross Validated. You would likely want to treat Year as a categorical variable. But if it were kept numeric at 2018-2022 then there is only one coefficient to describe that variable, and I don't see how you obtained an odds ratio that is below 1 if there is an increase every year. Maybe you could share your code and output. Commented Jul 16 at 21:01
• @rolando2 This about all I can share. I realize you would need the full data to properly see what's going on. But perhaps I'm making an obvious mistake in my code. Commented Jul 16 at 22:19

First, LogisticRegression is a penalized variant of logistic regression. It applies a penalty to the coefficients, so there is a good chance you won't get the same answer.

Let's use R instead. I'll set up the data as the authors have

year <- 2018:2022
repeat_visits <- c(2626, 3090, 3364, 4693, 4244)
single_visits <- c(7642, 7675, 8073, 9475, 8197)
n <- repeat_visits + single_visits


Now, I'll run a logistic regression to determine how the odds of a repeat visit changes over years. I'll represent years as the variable yrs_since_2018 <- year - 2018.

library(modelsummary)
yrs_since_2018 <- year - 2018
fit <- glm(
cbind(repeat_visits, n-repeat_visits) ~ yrs_since_2018,
family = binomial()
)

modelsummary(fit)



The log odds ratio is 0.103, which when exponentiated yields 1.108, meaning with each passing year the odds increases by a factor of 10%.

If you are intent on using python, you may instead want to use statsmodels. A similar logistic regression can be fit in the following way

import statsmodels.api as sm
import statsmodels.formula.api as smf
import pandas as pd
import numpy as np

df = pd.DataFrame({
'year': np.arange(2018, 2023) - 2018,
'repeat_visits': [2626, 3090, 3364, 4693, 4244],
'single_visits': [7642, 7675, 8073, 9475, 8197]
})

mod = smf.glm('repeat_visits + single_visits ~ year', family=sm.families.Binomial(), data=df).fit()

print(mod.summary())

==============================================================================================
Dep. Variable:     ['repeat_visits', 'single_visits']   No. Observations:                    5
Model:                                            GLM   Df Residuals:                        3
Model Family:                                Binomial   Df Model:                            1
Method:                                          IRLS   Log-Likelihood:                -29.962
Date:                                Tue, 16 Jul 2024   Deviance:                       11.699
Time:                                        19:27:31   Pearson chi2:                     11.7
No. Iterations:                                     4   Pseudo R-squ. (CS):              1.000
Covariance Type:                            nonrobust
==============================================================================
coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept     -1.0463      0.017    -62.030      0.000      -1.079      -1.013
year           0.1025      0.006     15.850      0.000       0.090       0.115
==============================================================================
$$$$

• Elegant solution: no weighting and a more interpretable baseline! Commented Jul 17 at 1:44
• This was very helpful. Thank you! Commented Jul 19 at 18:19

You should transform your data into long format and run a logit of visit type (coded 0 if one, 1 if more than one) on year (entered continuously) with the number of visits as a frequency weight. The effect on the OR is the exponentiated logit coefficient.

Here's what that might look like in Stata and Python:

. clear

.
. // Input the data for the years 2018 to 2022
. input year visits0 visits1

year     visits0     visits1
1. 2018 7642 2626
2. 2019 7675 3090
3. 2020 8073 3364
4. 2021 9475 4693
5. 2022 8197 4244
6. end

.
. // Reshape the data from wide to long format
. reshape long visits, i(year) j(vis_type)
(j = 0 1)

Data                               Wide   ->   Long
-----------------------------------------------------------------------------
Number of observations                5   ->   10
Number of variables                   3   ->   3
j variable (2 values)                     ->   vis_type
xij variables:
visits0 visits1   ->   visits
-----------------------------------------------------------------------------

.
. // Sort the data by visit type and year
. sort vis_type year

.
. // List the data without observation numbers, separated by visit type
. list, noobs sepby(vis_type)

+--------------------------+
| year   vis_type   visits |
|--------------------------|
| 2018          0     7642 |
| 2019          0     7675 |
| 2020          0     8073 |
| 2021          0     9475 |
| 2022          0     8197 |
|--------------------------|
| 2018          1     2626 |
| 2019          1     3090 |
| 2020          1     3364 |
| 2021          1     4693 |
| 2022          1     4244 |
+--------------------------+

.
. // Estimate a logistic regression model using GLM
. // The dependent variable is vis_type (assumed to be binary with values 0 or 1)
. // The independent variable is year
. // The 'fw' option specifies frequency weights, where 'visits' is the weight variable
. // The 'family(bernoulli)' option specifies the Bernoulli family for logistic regression
. // The 'nolog' option suppresses the iteration log
. // Uncomment eform to get exponentiated coefficients directly
. glm vis_type c.year [fw = visits], family(bernoulli) link(logit) nolog // eform

Generalized linear models                         Number of obs   =     59,079
Optimization     : ML                             Residual df     =     59,077
Scale parameter =          1
Deviance         =  72415.48807                   (1/df) Deviance =   1.225781
Pearson          =  59076.72434                   (1/df) Pearson  =   .9999953

Variance function: V(u) = u*(1-u)                 [Bernoulli]
Link function    : g(u) = ln(u/(1-u))             [Logit]

AIC             =   1.225808
Log likelihood   = -36207.74404                   BIC             =  -576641.7

------------------------------------------------------------------------------
|                 OIM
vis_type | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
year |   .1025484   .0064701    15.85   0.000     .0898673    .1152296
_cons |  -207.9891   13.07094   -15.91   0.000    -233.6076   -182.3705
------------------------------------------------------------------------------

.
. // Compute and display the exponentiated coefficient for the year variable
. nlcom (eff_on_OR:exp(_b[year]))

eff_on_OR: exp(_b[year])

------------------------------------------------------------------------------
vis_type | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
eff_on_OR |   1.107991   .0071688   154.56   0.000      1.09394    1.122042
------------------------------------------------------------------------------

.
. python
----------------------------------------------- python (type end to exit) -------------------------------------------------------------------------------------------------------------------------------------------
>>> import pandas as pd
>>> import numpy as np
>>> import statsmodels.api as sm
>>>
>>> # Input the data for the years 2018 to 2022
... data = {
...     'year': [2018, 2019, 2020, 2021, 2022],
...     'visits0': [7642, 7675, 8073, 9475, 8197],
...     'visits1': [2626, 3090, 3364, 4693, 4244]
... }
>>>
>>> # Create a DataFrame
... df = pd.DataFrame(data)
>>>
>>> # Reshape the data from wide to long format
... df_long = pd.melt(df, id_vars=['year'], value_vars=['visits0', 'visits1'], var_name='vis_type', value_name='visits')
>>>
>>> # Convert vis_type to binary (0 for visits0, 1 for visits1)
... df_long['vis_type'] = df_long['vis_type'].apply(lambda x: 0 if x == 'visits0' else 1)
>>>
>>> # Sort the data by visit type and year
... df_long = df_long.sort_values(by=['vis_type', 'year'])
>>>
>>> # Add a constant (intercept) to the independent variables
... df_long['constant'] = 1
>>>
>>> # Define the independent variables (constant and year)
... X = df_long[['constant', 'year']]
>>> y = df_long['vis_type']
>>> weights = df_long['visits']
>>>
>>> # Fit the model using statsmodels
>>> model_sm = sm.GLM(y, X_sm, family=sm.families.Binomial(), freq_weights=weights).fit()
>>>
>>> # Print the summary of the model
... print(model_sm.summary())
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable:               vis_type   No. Observations:                   10
Model:                            GLM   Df Residuals:                    59077
Model Family:                Binomial   Df Model:                            1
Method:                          IRLS   Log-Likelihood:                -36208.
Date:                Tue, 16 Jul 2024   Deviance:                       72415.
Time:                        17:56:43   Pearson chi2:                 5.91e+04
No. Iterations:                     5   Pseudo R-squ. (CS):              1.000
Covariance Type:            nonrobust
==============================================================================
coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const       -207.9891     13.071    -15.912      0.000    -233.608    -182.370
year           0.1025      0.006     15.850      0.000       0.090       0.115
==============================================================================
>>>
>>> # Exponentiate the coefficient for the year variable to get the odds ratio
... odds_ratio = np.exp(model_sm.params['year'])
>>> print("\nExponentiated Year Coefficient is the Multiplicative Effect of Year on Odds Ratio:", odds_ratio)

Exponentiated Year Coefficient is the Multiplicative Effect of Year on Odds Ratio: 1.1079909694224508
>>> end
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

.
end of do-file
`

The multiplicative effect on the OR is $$\exp(.1025484) = 1.1079909$$. Rounding, that's a 11% increase.

Testing that the coefficient on year is zero is equivalent to testing that the exponentiated coefficient is equal to 1 since adding zero and multiplying by 1 are equivalent. This means that you can use the p-value from the model output directly. It is $$\ll 0.01$$, which agrees with the paper.

• This was very helpful. Thank you! Commented Jul 19 at 18:19