# Unstandardizing coefficients from a GLMM

Winter term grades are in and I can once again scrape together some time to write blog posts! đźŽ‰

The last post I did about making added variable plots led me to think about other â€śget model resultsâ€ť topics, such as the one Iâ€™m talking about today: unstandardizing coefficients.

I find this comes up particularly for generalized linear mixed models (GLMM), where models donâ€™t always converge if explanatory variables are left unstandardized. The lack of convergence can be caused by explanatory variables with very different magnitudes, and standardizing the variables prior to model fitting can be useful. In such cases, coefficients and confidence interval limits will often need to be converted to their unstandardized values for interpretation. I donâ€™t find thinking about the change in mean response for a 1 standard deviation increase in a variable to be super intuitive, which is the interpretation of a standardized coefficient.

The math for converting the standardized slope estimates to unstandardized ones turns out to be fairly straightforward. Coefficients for each variable need to be divided by the standard deviation of that variable (this is only true for slopes, not intercepts). The math is shown here.

The first time I went though this process was quite clunky. Since then Iâ€™ve managed to tidy things up quite a bit through work with students, and things are now much more organized.

# R packages

Model fitting will be done via **lme4**, which is where Iâ€™ve most often needed to do this. Data manipulation tools from **dplyr** will be useful for getting results tidied up. Iâ€™ll also use helper functions from **purrr** to loop through variables and **broom** for the tidy extraction of fixed-effects coefficients from the model.

`library(lme4) # v. 1.1-15`

`## Loading required package: Matrix`

```
suppressPackageStartupMessages( library(dplyr) ) # v. 0.7.4
library(purrr) # 0.2.4
library(broom) # 0.4.3
```

# The dataset

The dataset Iâ€™ll use is named `cbpp`

, and comes with **lme4**. It is a dataset that has a response variable that is counted proportions, so the data will be analyzed via a binomial generalized linear mixed model.

`glimpse(cbpp)`

```
## Observations: 56
## Variables: 4
## $ herd <fct> 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5...
## $ incidence <dbl> 2, 3, 4, 0, 3, 1, 1, 8, 2, 0, 2, 2, 0, 2, 0, 5, 0, 0...
## $ size <dbl> 14, 12, 9, 5, 22, 18, 21, 22, 16, 16, 20, 10, 10, 9,...
## $ period <fct> 1, 2, 3, 4, 1, 2, 3, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3...
```

This dataset has no continuous explanatory variables in it, so Iâ€™ll add some to demonstrate standardizing/unstandardizing. I create three new variables with very different ranges.

```
set.seed(16)
cbpp = mutate(cbpp, y1 = rnorm(56, 500, 100),
y2 = runif(56, 0, 1),
y3 = runif(56, 10000, 20000) )
glimpse(cbpp)
```

```
## Observations: 56
## Variables: 7
## $ herd <fct> 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5...
## $ incidence <dbl> 2, 3, 4, 0, 3, 1, 1, 8, 2, 0, 2, 2, 0, 2, 0, 5, 0, 0...
## $ size <dbl> 14, 12, 9, 5, 22, 18, 21, 22, 16, 16, 20, 10, 10, 9,...
## $ period <fct> 1, 2, 3, 4, 1, 2, 3, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3...
## $ y1 <dbl> 547.6413, 487.4620, 609.6216, 355.5771, 614.7829, 45...
## $ y2 <dbl> 0.87758754, 0.33596115, 0.11495285, 0.26466003, 0.99...
## $ y3 <dbl> 14481.07, 11367.88, 14405.16, 18497.73, 17955.66, 10...
```

# Analysis

## Unstandardized model

Here is my initial generalized linear mixed model, using the three continuous explanatory variables as fixed effects and â€śherdâ€ť as the random effect. A warning message indicates that standardizing might be necessary.

```
fit1 = glmer( cbind(incidence, size - incidence) ~ y1 + y2 + y3 + (1|herd),
data = cbpp, family = binomial)
```

```
## Warning: Some predictor variables are on very different scales: consider
## rescaling
```

```
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.947357 (tol =
## 0.001, component 1)
```

```
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?
```

## Standardizing the variables

Iâ€™ll now standardize the three explanatory variables, which involves subtracting the mean and then dividing by the standard deviation. The `scale()`

function is one way to do this in R.

I do the work inside `mutate_at()`

, which allows me to choose the three variables I want to standardize by name and add â€śsâ€ť as the suffix by using a name in `funs()`

. Adding the suffix allows me to keep the original variables, as I will need them later. I use `as.numeric()`

to convert the matrix that the `scale()`

function returns into a vector.

```
cbpp = mutate_at(cbpp, vars( y1:y3 ), funs(s = as.numeric( scale(.) ) ) )
glimpse(cbpp)
```

```
## Observations: 56
## Variables: 10
## $ herd <fct> 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5...
## $ incidence <dbl> 2, 3, 4, 0, 3, 1, 1, 8, 2, 0, 2, 2, 0, 2, 0, 5, 0, 0...
## $ size <dbl> 14, 12, 9, 5, 22, 18, 21, 22, 16, 16, 20, 10, 10, 9,...
## $ period <fct> 1, 2, 3, 4, 1, 2, 3, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3...
## $ y1 <dbl> 547.6413, 487.4620, 609.6216, 355.5771, 614.7829, 45...
## $ y2 <dbl> 0.87758754, 0.33596115, 0.11495285, 0.26466003, 0.99...
## $ y3 <dbl> 14481.07, 11367.88, 14405.16, 18497.73, 17955.66, 10...
## $ y1_s <dbl> 0.34007250, -0.32531152, 1.02536895, -1.78352139, 1....
## $ y2_s <dbl> 1.4243017, -0.4105502, -1.1592536, -0.6520950, 1.837...
## $ y3_s <dbl> -0.2865770, -1.3443309, -0.3123665, 1.0781427, 0.893...
```

## Standardized model

The model with standardized variables converges without any problems.

```
fit2 = glmer( cbind(incidence, size - incidence) ~ y1_s + y2_s + y3_s + (1|herd),
data = cbpp, family = binomial)
```

# Unstandardizing slope coefficients

## Get coefficients and profile confidence intervals from model

If I want to use this model for inference I need to unstandardize the coefficients before reporting them to make them more easily interpretable.

The first step in the process is to get the standardized estimates and confidence intervals from the model. I use `tidy()`

from package **broom** for this, which returns a data.frame of coefficients, statistical tests, and confidence intervals. The help page is at `?tidy.merMod`

if you want to explore some of the options.

I use `tidy()`

to extract the fixed effects along with profile likelihood confidence intervals.

```
coef_st = tidy(fit2, effects = "fixed",
conf.int = TRUE,
conf.method = "profile")
```

`## Computing profile confidence intervals ...`

`coef_st`

```
## term estimate std.error statistic p.value conf.low
## 1 (Intercept) -2.0963512 0.2161353 -9.6992521 3.037166e-22 -2.565692724
## 2 y1_s 0.2640116 0.1395533 1.8918334 5.851318e-02 -0.009309633
## 3 y2_s 0.1031655 0.1294722 0.7968153 4.255583e-01 -0.151459144
## 4 y3_s 0.1569910 0.1229460 1.2769107 2.016338e-01 -0.086291865
## conf.high
## 1 -1.6536063
## 2 0.5472589
## 3 0.3635675
## 4 0.4033177
```

## Calculate standard deviations for each variable

I need the standard deviations for each variable in order to unstandardize the coefficients. If I do this right, I can get the standard deviations into a data.frame that I can then join to `coef_st`

. Once that is done, dividing the estimated slopes and confidence interval limits by the standard deviation will be straightforward.

I will calculate the standard deviations per variable with `map()`

from **purrr**, as it is a convenient way to loop through columns. I pull out the variables I want to calculate standard deviations for via `select()`

. An alternative approach would have been to take the variables from columns and put them in rows (i.e., put the data in *long* format), and then summarize by groups.

The output from `map()`

returns a list, which can be stacked into a long format data.frame via `utils::stack()`

. This results in a two column data.frame, with a column for the standard deviation (called `values`

) and a column with the variable names (called `ind`

).

```
map( select(cbpp, y1:y3), sd) %>%
stack()
```

```
## values ind
## 1 90.4430192 y1
## 2 0.2951881 y2
## 3 2943.2098667 y3
```

The variables in my model and in my output end with `_s`

, so Iâ€™ll need to add that suffix to the variable names in the â€śstandard deviationsâ€ť dataset prior to joining the two data.frames together.

```
sd_all = map( select(cbpp, y1:y3), sd) %>%
stack() %>%
mutate(ind = paste(ind, "s", sep = "_") )
sd_all
```

```
## values ind
## 1 90.4430192 y1_s
## 2 0.2951881 y2_s
## 3 2943.2098667 y3_s
```

## Joining the standard deviations to the coefficients table

Once the names of the variables match between the datasets I can join the â€śstandard deviationsâ€ť data.frame to the â€ścoefficientsâ€ť data.frame. Iâ€™m not unstandardizing the intercept at this point, so Iâ€™ll use `inner_join()`

to keep only rows that have a match in both data.frames. Notice that the columns Iâ€™m joining by have different names in the two data.frames.

```
coef_st %>%
inner_join(., sd_all, by = c("term" = "ind") )
```

```
## term estimate std.error statistic p.value conf.low conf.high
## 1 y1_s 0.2640116 0.1395533 1.8918334 0.05851318 -0.009309633 0.5472589
## 2 y2_s 0.1031655 0.1294722 0.7968153 0.42555829 -0.151459144 0.3635675
## 3 y3_s 0.1569910 0.1229460 1.2769107 0.20163378 -0.086291865 0.4033177
## values
## 1 90.4430192
## 2 0.2951881
## 3 2943.2098667
```

With everything in one data.frame I can easily divide `estimate`

, `conf.low`

and `conf.high`

by the standard deviation in `values`

via `mutate_at()`

. I will round the results, as well, although Iâ€™m ignoring the vast differences in the variable range when I do this rounding.

```
coef_st %>%
inner_join(., sd_all, by = c("term" = "ind") ) %>%
mutate_at( vars(estimate, conf.low, conf.high), funs(round( ./values, 4) ) )
```

```
## term estimate std.error statistic p.value conf.low conf.high
## 1 y1_s 0.0029 0.1395533 1.8918334 0.05851318 -0.0001 0.0061
## 2 y2_s 0.3495 0.1294722 0.7968153 0.42555829 -0.5131 1.2316
## 3 y3_s 0.0001 0.1229460 1.2769107 0.20163378 0.0000 0.0001
## values
## 1 90.4430192
## 2 0.2951881
## 3 2943.2098667
```

Iâ€™ll get rid of the extra variables via `select()`

, so I end up with the unstandardized coefficients and confidence interval limits along with the variable name. I could also get the variable names cleaned up, possibly removing the suffix and/or capitalizing and adding units, etc., although I donâ€™t do that with these fake variables today.

```
coef_unst = coef_st %>%
inner_join(., sd_all, by = c("term" = "ind") ) %>%
mutate_at( vars(estimate, conf.low, conf.high), funs(round( ./values, 4) ) ) %>%
select(-(std.error:p.value), -values)
coef_unst
```

```
## term estimate conf.low conf.high
## 1 y1_s 0.0029 -0.0001 0.0061
## 2 y2_s 0.3495 -0.5131 1.2316
## 3 y3_s 0.0001 0.0000 0.0001
```

## Estimates from the unstandardized model

Note that the estimated coefficients are the same from the model where I manually unstandardize the coefficients (above) and the model fit using unstandardized variables.

`round( fixef(fit1)[2:4], 4)`

```
## y1 y2 y3
## 0.0029 0.3495 0.0001
```

Given that the estimates are the same, couldnâ€™t we simply go back and fit the unstandardized model and ignored the warning message? Unfortunately, the convergence issues can cause problems when trying to calculate profile likelihood confidence intervals, so the simpler approach doesnâ€™t always work.

In this case there are a bunch of warnings (not shown), and the profile likelihood confidence interval limits arenâ€™t successfully calculated for some of the coefficients.

```
tidy(fit1, effects = "fixed",
conf.int = TRUE,
conf.method = "profile")
```

`## Computing profile confidence intervals ...`

```
## term estimate std.error statistic p.value
## 1 (Intercept) -4.582545e+00 1.094857e+00 -4.1855190 2.845153e-05
## 2 y1 2.919191e-03 1.520407e-03 1.9200071 5.485700e-02
## 3 y2 3.495238e-01 4.370778e-01 0.7996834 4.238943e-01
## 4 y3 5.334512e-05 4.049206e-05 1.3174218 1.876973e-01
## conf.low conf.high
## 1 NA NA
## 2 NA NA
## 3 NA NA
## 4 -2.931883e-05 0.0001370332
```

# Further (important!) work

These results are all on the scale of log-odds, and I would exponentiate the unstandardized coefficients to the odds scale for reporting and interpretation.

Along these same lines, one thing I didnâ€™t discuss in this post that is important to consider is the appropriate and interesting unit increase for each variable. Clearly the effect of a â€ś1 unitâ€ť increase in the variable is likely not of interest for at least `y2`

(range between 0 and 1) and `y3`

(range between 10000 and 20000). In the first case, 1 unit encompasses the entire range of the variable and in the second case 1 unit appears to be much smaller than the scale of the measurement.

The code to calculate the change in odds for a practically interesting increase in each explanatory variable would be similar to what Iâ€™ve done above. I would create a data.frame with the unit increase of interest for each variable in it, join this to the â€ścoefficientsâ€ť dataset, and multiply all estimates and CI by those values. The multiplication step can occur before or after unstandardizing but must happen before doing exponentiation/inverse-linking. Iâ€™d report the unit increase for each variable in any tables of results so the reader can see that the reported estimate is a change in estimated odds/mean for the given practically important increase in the variable.