Getting started

This session, we will focus on modelling of count outcomes in epidemiology.

We will learn how to:

• Plot the relationships between count variables and other variables
• Run and interpret regression models
• Interpret the model fit and summary statistics of models
• Check models for over-dispersion
• Consider the effects of multi-collinearity on exposure-disease relationships
• Use regression nomograms to better understand and communicate the results of regression models
  
  

We will be leaning on the following libraries to do these analyses:

• regplot
• PerformanceAnalytics
• jtools
• corrplot
• finalfit and
• car

Background

We will be diving into a real-life research battle ground today.

Research battle by Sander Meyer from unsplash.com

There is some controversy about the cause of acute rheumatic fever, which is an overall rare disease, but it occurs much more commonly in Māori and Pacific children than children of any other ethnic group in NZ.

This disease often presents with fever and inflamed joints and can lead to life-limiting sequelae which include thickened and dysfunctional heart valves. If these aren't bad enough, it may lead to the need for cardiac surgery and in some cases, premature death.

Thickened heart valves associated with rheumatic heart disease

The most prevalent explanation in the health sector is that it is caused by bacterial streptococcal throat infection.

An alternative explanation is that it is caused by the common skin condition, scabies, which has the following skin appearance, which is easy not to notice. The little crops of raised lumps are scabies and they are caused by the microscopic mite pictured above.

Typical skin appearance of scabies lesions: image courtesy of Dr Daniel Engelman

The treatment for scabies in New Zealand is pictured below, and we will use records of pharmacy dispensing of this drug as an indirect measure of scabies prevalence in the community.

The proposed mechanism by which scabies causes streptococcal infection and rheumatic fever is given below:

Consequences of scabies infection

We will be looking for the evidence that prescribing for a certain treatment for scabies (the lotion permethrin) is linked with acute rheumatic fever at a suburb level (1 row of data for each suburb). Since the counts of acute rheumatic fever over a certain period constitute a 'count' variable, Poisson regression will be used. The population of the suburb will be used as an offset in the model. An offset is often used as the denominator when rate values are modelled. For example, here, our outcome is the count or number of permethrin scripts dispensed in a suburb divided by the population of that suburb. The count of prescriptions is naturally influenced by the population, so this information is incorporated in the model as an offset. The offset allows us to model rates rather than counts. Technically, it is a covariate or independent variable whose value of beta coefficient is forced to be one.

Covariates, or potential confounders in the analysis include:

• The proportion of Māori and Pacific people living in a suburb
• Socioeconomic deprivation

We will use the counts of acute rheumatic fever by suburb (or census area unit). This information is derived from routine datasets used in the health sector. The formal write-up of the analysis is here.

Set-up

To set-up our session for analysis and call the necessary libraries, we’ll use pacman::p_load() instead of install.packages() and library(). The pacman library shortens this procedure, and p_load() will look for a package and install it if it can't be found. It will also put it on the search path, like the library() command does. We won't run this code, as it will mess up our tutorial, but it is essential for getting this running on a local machine or posit cloud.

if(!require(pacman)) install.packages("pacman")
pacman::p_load(tidyverse, 
               epiDisplay, 
               visdat, # missing data visualisation
               skimr, # summarise data frame
               lattice, # basic histograms
               jtools, # exponentiate beta coefficients
               ggstance,
               broom.mixed,
               finalfit, # final tables
               dplyr,
               kableExtra,
               car,
               MixAll,
               regplot, # nomograms
               PerformanceAnalytics,
               corrplot, 
               performance) # overdispersion check

View the dataframe; called for_shiny. Notice that there’s 1 row for each suburb:

for_shiny

for_shiny |> names()

Rename our dataframe as df. Then rename the variables so they're more understandable. Although this seems tedious, it really makes life easier and reduces confusion.

df <- for_shiny # Rename dataframe

# Rename variables
names(df) <- c("arf_count", "population", "Pacific_proportion",
               "Maori_proportion", "Permethrin_rate", "Suburb",
               "nzdep_decile", "arf_rate")
df |> names() # Check that the variables have been renamed

stopifnot(df[df$Suburb |> duplicated(),] == 0)

df |> skim()

var(df$arf_count)

df |> vis_miss(sort_miss = TRUE)

Univariable distributions: histograms

Use the histogram() function for these.

df$arf_count |> histogram()

df$arf_rate |> histogram()

df$Permethrin_rate |> histogram()

df$Pacific_proportion 
df$Maori_proportion
df$nzdep_decile 

df$Pacific_proportion |> histogram()
df$Maori_proportion |> histogram()
df$nzdep_decile |> histogram()

# 5 groups instead of 10
df$nzDepQuint <- cut(df$nzdep_decile, breaks = c(-Inf, 2, 4, 6, 8, Inf),
                     labels = c("1 to 2 (least deprived)",
                                "3 to 4",
                                "5 to 6",
                                "7 to 8",
                                "9 to 10 (most deprived)"))

# Check the frequency of our 5 groups
epiDisplay::tabpct(df$nzDepQuint, df$nzdep_decile)

Bivariate associations: scatterplots

Before undertaking any regression analysis, it is useful to visualise the relationships between key exposure, confounding and outcome variables.

Use the car::scatterplot() function.

car::scatterplot(arf_rate ~ Permethrin_rate, data = df)

arf_rate ~ nzdep_decile
arf_rate ~ Pacific_proportion 
arf_rate ~ Maori_proportion 

car::scatterplot(arf_rate ~ nzdep_decile, data = df)
car::scatterplot(arf_rate ~ Pacific_proportion, data = df)
car::scatterplot(arf_rate ~ Maori_proportion, data = df)

## Select variables of interest into new data.frame
my_data <- df[, c("Pacific_proportion", 
                  "Maori_proportion", 
                  "Permethrin_rate", 
                  "nzdep_decile", 
                  "arf_rate")]

## Chart the correlations and regressions 
## and distributions of each variable.
PerformanceAnalytics::chart.Correlation(my_data, 
                                        histogram = TRUE, 
                                        pch = 19)

Regression

We will build a Poisson regression model with our outcome as arf_count and exposure as Permethrin_rate. To adjust for population size, we will include log_population as an offset. We will also adjust for NZDep quintile (nzDepQuint), Pacific_proportion and Maori_proportion.

Run the following code using the glm() (generalised linear model) function to make a Poisson regression model of arf_count with all other risk factors as independent variables or covariates.

## Make offset variable
df$log_population <- log(df$population)

## Build model
fit <- glm(arf_count ~ Permethrin_rate +
           nzDepQuint +
           Pacific_proportion +
           Maori_proportion, 
           offset = log_population,
           family = poisson(link = "log"),
           data = df)

## Print summary of model
summary(fit)

## check for overdispersion
performance::check_overdispersion(fit)

Look at the term, offset = log_population, in the model. With counts (arf_count) as the outcome, having log_population as the offset makes the model represent a rate (count per unit population). The overdispersion test is a test of the residual deviance differing substantially from the degrees of freedom in the model. Here, it greater than 0.05, so not significant. Hence, we can stick with our Poisson model.

Check your understanding.

jtools::summ(fit, exp = TRUE)

# fit model
fit1 <- glm(arf_count ~ I(Permethrin_rate/1000) +
             nzDepQuint +
             Pacific_proportion +
             Maori_proportion + offset(log(population)),
           family = poisson(link = "log"),
           data = df)

# view characteristics of model
summary(fit1)

# fit model
fit_qp <- glm(arf_count ~ I(Permethrin_rate/1000) +
             nzDepQuint +
             Pacific_proportion +
             Maori_proportion + offset(log(population)),
           family = quasipoisson(link = "log"),
           data = df)

# view characteristics of model
summary(fit_qp)

# jtools::summ(model_object, exp = logical_value)

jtools::summ(fit_qp, exp = TRUE)

jtools::plot_summs(fit_qp, exp = TRUE)

Multi-collinearity

In a model, when 2 or more predictors are strongly correlated with one other, this is called multi-collinearity. In a model, it may be problematic, since:

• increases the variance of estimates of regression coefficients
• makes these estimates unstable and unreliable
  

We can detect it by first looking for strong correlations between variables.

Use the cor() function to estimate the correlation between Pacific_proportion and Permethrin_rate:

cor(df$Pacific_proportion, df$Permethrin_rate)

You should get a correlation of 0.88, which is strong. Thus, multi-collinearity in the predictor variables may be a problem for our analysis.

To deal with this, we can remove one of these predictors from the model. So that we retain our exposure (Permethrin_rate) in our model, let’s remove Pacific_proportion by adding # before it.

Another change is that we'll divide Permethrin_rate by 1000. This will increase the coefficient for this variable by 1000 so that it’s easier to interpret. (This is like converting grams into kilograms; the latter is more meaningful for larger weights). Below, we use the I() function to do this calculation within the 'model' argument.
Now run our revised code below:

df$log_population <- log(df$population)

fit3 <- glm(arf_count ~ I(Permethrin_rate/1000) +
            nzDepQuint +
            # Pacific_proportion +
            Maori_proportion, offset = log(population),
            family = poisson(link = "log"),
            data = df)

summary(fit3)

What is the beta-coefficient or slope for the Permethrin_rate/1000 term?

Is it statistically significant?

What is the incidence rate ratio for I(Permethrin_rate/1000)? We need to exponentiate the beta-coefficient.

This next model has Pacific_proportion included again. so we can test for the effect of multi-collinearity. Run the code below.

fit5 <- glm(arf_count ~ I(Permethrin_rate/1000) +
              nzDepQuint +
              Pacific_proportion +
              Maori_proportion, offset = log_population,
              family = poisson(link = "log"),
              data = df)

summary(fit5)

Because of multicollinearity, we expect the standard error for I(Permethrin_rate/1000) to increase after Pacific_proportion has been added. Has it?

Regression nomogram

df$Suburb |> print()

The code below will run in posit, but unfortunately not on this platform, so we will just show you the code and the output. We have plotted the values for Otara East, highlighted on the figure in red.

regplot::regplot(fit, 
                 observation = df[df$Suburb == "Otara East",], clickable = FALSE)

## Error in `[.data.frame`(observation, reord): undefined columns selected

Regplot

For Otara East, for example, we see that Pacific_proportion is about 0.8 or 80%. Also, the dispensing rate of permethrin Permethrin_rate is relatively high, as expected from the correlation coefficient between this variable and Pacific_proportion as we discussed earlier. Each variable contributes a score (β(X-m) terms) and these add up as a type of weighted average to produce a predicted rate of acute rheumatic fever (arf_count). Note that for Otara East, the predicted arf_count is very high.

Please attempt the next 2 questions.

Question 19

Hint: In each density plot or histogram, look at where the red dot (for Otara East) is.

regplot::regplot(fit, observation = df[df$Suburb == "Remuera South",])

## Error in `[.data.frame`(observation, reord): undefined columns selected

Regplot2

Model fit

How well a model fits the data can be assessed with a 'marginal model plot', which shows the relationship between model predicted (model) and observed data (data). In a good fitting model, we expect to see the two lines closely superimposed on each other. Poor fit is identified by seeing separation between the two lines, which may indicate over or under prediction of the model compared to the data.

## The assign commands are only needed in 
## the shiny app to get the marginalModelPlot function working
assign("fit", fit, envir = .GlobalEnv)
assign("df", df, envir = .GlobalEnv)

## model fit plot----
car::marginalModelPlots(fit)

What is the nature of the overall fit of the data to the model?
- Good fit is indicated by good alignment between Data (blue line) and the Model (red line)
- Poor fit is indicated by poor alignment between these lines
- Often, for more extreme or outlying values of an exposure, there is deviation of the Data from the Model. This may not invalidate your model, but would be important to be aware of, and may be worthwhile discussing in your results.

Publication-quality regression tables

In published articles, it's common to present univariable (unadjusted) and multivariable (adjusted) regression coefficients in a table. To create these, run the following code.

df$Permethrin_rate_per_1000 <- df$Permethrin_rate/1000  # Like converting grams into kilograms
df$log_population <- log(df$population)  # Offset

dependent <- "arf_count"
explanatory <- c("Pacific_proportion", "Maori_proportion", 
                 "Permethrin_rate_per_1000","nzdep_decile")

# `fit_uni` creates univariable (unadjusted) models
fit_uni = df |> 
  glmuni(dependent, explanatory, family = quasipoisson,
         offset = df$log_population) |> 
  fit2df(estimate_name = "Rate ratio (univariable)")

# `fit_multi` creates a multivariable (adjusted) model
fit_multi = df |> 
  glmmulti(dependent, explanatory, family = quasipoisson, 
           offset = df$log_population) |> 
  fit2df(estimate_name = "Rate ratio (multivariable)")

df |>
  ## Crosstable
  summary_factorlist(dependent, explanatory, cont = "median", fit_id=TRUE)  |> 
  
  #### Add univariable table
  ff_merge(fit_uni, estimate_name = "Rate ratio") |> 
  
  #### Add multivariable table
  ff_merge(fit_multi, estimate_name = "Rate ratio",
           last_merge = TRUE) |> 
dependent_label(df, dependent) -> regression_table

regress_tab <- knitr::kable(regression_table, row.names = FALSE, 
                            align = c("l", "l", "r", "r", "r"))

#### Puts output into viewer so that it can be cut and pasted into Excel, 
#### then Word.
regress_tab |> kableExtra::kable_styling()

Notice that the adjusted Rate ratio (multivariable) for Permethrin_rate_per_1000 is not statistically significant. Does that mean that permethrin or scabies is not related to acute rheumatic fever incidence?

Summary

We have learned a number of things related to count models.

• The outcome is a count (non-negative integer) and this is transformed in the model using a natural logarithm.
• We use offsets to convert count outcomes to rates by incorporating the denominator into our model.
• The beta coefficients or slopes for each exposure or confounder term need to be exponentiated to be readily interpretable as incidence rate ratios which are interpreted similarly to risk ratios or odds ratios (jtools::summ(model, exp = TRUE)).
• We usually start with a Poisson model and check for over or under dispersion (mean and variance unequal) and may use a quasi-Poisson or negative binomial to accommodate this feature of the data.
• We can check for model fit using marginal model plots.
• Collinearity or strong associations within exposure or independent variables can be a feature of any regression analysis and can be checked for using correlation coefficients and scatterplots (PerformanceAnalytics::chart.Correlation()). The inclusion or exclusion of colinear variables can make a big difference to the results of our regression!!
• We often present models in scientific papers by comparing crude and adjusted modelled incidence rate ratios (finalfit).

Homework

Select only suburbs that have a socioeconomic deprivation of 6 or more:
df_dep <- df[df$nzdep_decile >= 6, ]

Then, for this subgroup, repeat the analyses we’ve covered in this tutorial.

Please have a go at these tasks in posit (desktop) or posit cloud.