In GHBolstad/lmmHeterosced: Linear mixed model with heteroscedastisity

knitr::opts_chunk$set(echo = TRUE)
library(lmmHeterosced)

Simulated data

To illustrate the use of the model we simulate data with two explanatory variables and 10 categories with 50 data points within each category. The two explanatory variables influence both the expectation and variance of y.

n_cat <- 10
n <- 50
d <- data.frame(x1 = rnorm(n*n_cat, 2, 1), 
                x2=rnorm(n*n_cat, 1, 1), 
                cat = rep(1:n_cat, each = n),
                SE = runif(n*n_cat,0.1,0.9))
d$y <- rnorm(nrow(d), 1+0.5*d$x1+0.1*d$x2, exp(0.5*(1+0.5*d$x1+0.6*d$x2))) + # within category effect including heteroscedastity
        rep(rnorm(n_cat), each = n) + # random effect of category
        rnorm(nrow(d),0, d$SE) # measurment error

par(mfrow=c(2,1))
with(d, plot(x1,y))
with(d, plot(x2,y))
par(mfrow=c(1,1))

Centring

If there is more than one explanatory variable, it is recommendable to mean center the all explanatory variable before the analysis. Otherwise the output of the plotting functions will look strange because of the intercept. With only one explanatory variable it makes no difference.

d$x1c <- d$x1 - mean(d$x1)
d$x2c <- d$x2 - mean(d$x2)

The model

mod <- lmmHeterosced(formula = y~ x1c + x2c +(1|cat), heterosced_formula = ~ x1c + x2c, SE = d$SE, data = d)

Results of average change in y

mod$Fixef

Results of Heteroscedasticity analysis

The intercept gives the log residual variance at x1 = x2 = 0 and the linear parameters (x1 and x2) gives the change in log residual variance with x

mod$Heterosced

AICc

mod$AICc

Plot of the model

The plotting function only handles one explanatory variable at the time. The black solid line give the change in the average, and blue lines shows how the heteroscedasticity changes by showing +/- one standard deviation of the residual distribution.

par(mfrow=c(2,1))
plot(mod, "x1c", col_data = "grey", col_expectation = "black", col_sd = "blue")
plot(mod, "x2c", col_data = "grey", col_expectation = "black", col_sd = "blue")
par(mfrow=c(1,1))

Change in absolute y-values

Plot that show the expected change in the absolute values of the response variable. The average change expected change over all observed values of x is given by the inset with 95\% confidence interval in parenthesis.

par(mfrow=c(2,1))
plot_abs_y(mod, x = "x1c", col = "grey", line_col = "black", cex.legend = 0.9) 
plot_abs_y(mod, x = "x2c", col = "grey", line_col = "black", cex.legend = 0.9) 
par(mfrow=c(1,1))