[R-sig-ME] Unstandardizing GLMM Regression Coefficients

[Timothy Lau]

Hello Thierry,

My apologies for the HTML formatting, I am guessing gmail did that
because I italicized the the code sections for readability.

Below is the unitalicized version:

**************************************
Hello,

lme4 mentioned that I should consider rescaling some of my predictor
variables because of the scale differences:

fit warnings:
Some predictor variables are on very different scales: consider rescaling



I have since been trying to expand Ben Bolker's function:

# http://stackoverflow.com/questions/23642111/how-to-unscale-the-coefficients-from-an-lmer-model-fitted-with-a-scaled-respon
rescale.coefs <- function(beta,mu,sigma) {
  beta2 <- beta ## inherit names etc.
  beta2[-1] <- sigma[1]*beta[-1]/sigma[-1]
  beta2[1]  <- sigma[1]*beta[1]+mu[1]-sum(beta2[-1]*mu[-1])
  beta2
}
# regular model
m1 <- lm(Illiteracy~.,as.data.frame(state.x77))
b1 <- coef(m1)
# Make a scaled version of the data
ss <- scale(state.x77)
# Scaled coefficients:
m1S <- update(m1,data=as.data.frame(ss))
b1S <- coef(m1S)
# rescaling
icol <- which(colnames(state.x77)=="Illiteracy")
p.order <- c(icol,(1:ncol(state.x77))[-icol])
m <- colMeans(state.x77)[p.order]
s <- apply(state.x77,2,sd)[p.order]
all.equal(b1,rescale.coefs(b1S,m,s))  ## TRUE



to work with GLMM models (e.g., count data) that can't have the
outcome standardized, only the predictors:
# function assumes scaled data is within model; no NA's and only used
variables are in orig.data
unscale.coef.mer <- function(model, orig.data) {
  require(lme4)
  ran.ef <- as.character(attr(attr(model at frame,
"terms"),"predvars.random"))[-1:-2] # random effects
  offset <- attr(attr(model at frame, "terms"), "offset")
  scaled.data <- model at frame[,!names(model at frame) %in% c(ran.ef,
offset)] # the scaled covariate(s) and outcome used in the model
  orig.data <- orig.data[,!names(orig.data) %in% c(ran.ef, offset)]
  beta <- fixef(model) # the fixed effects
  outvar <- names(model at frame)[attr(attr(m3 at frame, "terms"),
"response")] # the outcoem var. name
  icol <- which(colnames(orig.data) == outvar) # outcome variable column number
  p.order <- c(icol, (1:ncol(orig.data))[-icol]) # full list of
variable names with outcome 1st
  mu <- colMeans(x = orig.data)[p.order] # variable means with outcome 1st
  sigma <- apply(X = orig.data, MARGIN = 2, FUN = sd)[p.order] #
variable SDs with outcome 1st
  beta2 <- beta # inherit names etc.
  beta2[-1] <- sigma[1] * beta[-1] / sigma[-1] # the fixed effects
except intercept
  beta2[1]  <- sigma[1] * beta[1] + mu[1] - sum(beta2[-1] * mu[-1]) #
the FE intercept
  return(beta2)
}


Could someone help me with backing out or unscaling the coefficeints
so that m2 coef matches m1 even though m2 used the standardized data
(example data):
# warnings due to scaling
m1 <- glmer(formula = awards ~ 1 + write + (1 | cid), family =
poisson, data = foreign::read.dta(file =
"http://www.ats.ucla.edu/stat/data/hsbdemo.dta"))
# fixed scaling but need to back out coefficients to original scale
for interpretative reasons
m2 <- glmer(formula = awards ~ 1 + scale(write) + (1 | cid), family =
poisson, data = foreign::read.dta(file =
"http://www.ats.ucla.edu/stat/data/hsbdemo.dta"))



Best,
Timothy



"The will to win means nothing without the will to prepare."
http://en.wikipedia.org/wiki/Juma_Ikangaa


On Mon, Jan 18, 2016 at 3:10 AM, Thierry Onkelinx
<thierry.onkelinx at inbo.be> wrote:
> Dear Timothy,
>
> Please don't post in HTML. It makes code unreadable.
> Please post code that can be easily reproduced: a simple copy/paste should
> work.
>
> You might want to consider a very simple scaling by just using expressing
> the variables in another unit (= multiplying by some power of 10). E.g.
> kilometers to meters or gram to decagram.
>
> Best regards,
>
> Thierry
>
>
> ir. Thierry Onkelinx
> Instituut voor natuur- en bosonderzoek / Research Institute for Nature and
> Forest
> team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
> Kliniekstraat 25
> 1070 Anderlecht
> Belgium
>
> To call in the statistician after the experiment is done may be no more than
> asking him to perform a post-mortem examination: he may be able to say what
> the experiment died of. ~ Sir Ronald Aylmer Fisher
> The plural of anecdote is not data. ~ Roger Brinner
> The combination of some data and an aching desire for an answer does not
> ensure that a reasonable answer can be extracted from a given body of data.
> ~ John Tukey
>
> 2016-01-18 7:52 GMT+01:00 Timothy Lau <timothy.s.lau at gmail.com>:
>>
>> Hello,
>>
>> lme4 mentioned that I should consider rescaling some of my predictor
>> variables because of the scale differences:
>>
>> *fit warnings:*
>> *Some predictor variables are on very different scales: consider
>> rescaling*
>>
>>
>>
>> I have since been trying to expand Ben Bolker's function:
>>
>> *#
>>
>> http://stackoverflow.com/questions/23642111/how-to-unscale-the-coefficients-from-an-lmer-model-fitted-with-a-scaled-respon
>>
>> <http://stackoverflow.com/questions/23642111/how-to-unscale-the-coefficients-from-an-lmer-model-fitted-with-a-scaled-respon>*
>> *rescale.coefs <- function(beta,mu,sigma) {*
>> *  beta2 <- beta ## inherit names etc.*
>> *  beta2[-1] <- sigma[1]*beta[-1]/sigma[-1]*
>> *  beta2[1]  <- sigma[1]*beta[1]+mu[1]-sum(beta2[-1]*mu[-1])*
>> *  beta2*
>> *}*
>> *# regular model*
>> *m1 <- lm(Illiteracy~.,as.data.frame(state.x77))*
>> *b1 <- coef(m1)*
>> *# Make a scaled version of the data*
>> *ss <- scale(state.x77)*
>> *# Scaled coefficients:*
>> *m1S <- update(m1,data=as.data.frame(ss))*
>> *b1S <- coef(m1S)*
>> *# rescaling*
>> *icol <- which(colnames(state.x77)=="Illiteracy")*
>> *p.order <- c(icol,(1:ncol(state.x77))[-icol])*
>> *m <- colMeans(state.x77)[p.order]*
>> *s <- apply(state.x77,2,sd)[p.order]*
>> *all.equal(b1,rescale.coefs(b1S,m,s))  ## TRUE*
>>
>>
>>
>> to work with GLMM models (e.g., count data) that can't have the outcome
>> standardized, only the predictors:
>> *# function assumes scaled data is within model; no NA's and only used
>> variables are in orig.data*
>> *unscale.coef.mer <- function(model, orig.data) {*
>> *  require(lme4)*
>> *  ran.ef <- as.character(attr(attr(model at frame,
>> "terms"),"predvars.random"))[-1:-2] # random effects*
>> *  offset <- attr(attr(model at frame, "terms"), "offset")*
>> *  scaled.data <- model at frame[,!names(model at frame) %in% c(ran.ef, offset)]
>> # the scaled covariate(s) and outcome used in the model*
>> *  orig.data <- orig.data[,!names(orig.data) %in% c(ran.ef, offset)]*
>> *  beta <- fixef(model) # the fixed effects*
>> *  outvar <- names(model at frame)[attr(attr(m3 at frame, "terms"), "response")]
>> # the outcoem var. name*
>> *  icol <- which(colnames(orig.data) == outvar) # outcome variable column
>> number*
>> *  p.order <- c(icol, (1:ncol(orig.data))[-icol]) # full list of variable
>> names with outcome 1st*
>> *  mu <- colMeans(x = orig.data)[p.order] # variable means with outcome
>> 1st*
>> *  sigma <- apply(X = orig.data, MARGIN = 2, FUN = sd)[p.order] # variable
>> SDs with outcome 1st*
>> *  beta2 <- beta # inherit names etc.*
>> *  beta2[-1] <- sigma[1] * beta[-1] / sigma[-1] # the fixed effects except
>> intercept*
>> *  beta2[1]  <- sigma[1] * beta[1] + mu[1] - sum(beta2[-1] * mu[-1]) # the
>> FE intercept*
>> *  return(beta2)*
>> *}*
>>
>>
>> Could someone help me with backing out or unscaling the coefficeints so
>> that m2 coef matches m1 even though m2 used the standardized data (example
>> data):
>> *# warnings due to scaling*
>> *m1 <- glmer(formula = awards ~ 1 + write + (1 | cid), family = poisson,
>> data = foreign::read.dta(file =
>> "http://www.ats.ucla.edu/stat/data/hsbdemo.dta
>> <http://www.ats.ucla.edu/stat/data/hsbdemo.dta>"))*
>>
>> *# fixed scaling but need to back out coefficients to original scale for
>> interpretative reasonsm2 <- glmer(formula = awards ~ 1 +
>> scale(write) + (1 | cid), family = poisson, data =
>> foreign::read.dta(file = "http://www.ats.ucla.edu/stat/data/hsbdemo.dta
>> <http://www.ats.ucla.edu/stat/data/hsbdemo.dta>"))*
>>
>>
>>
>> Best,
>> Timothy
>>
>>
>>
>> "The will to win means nothing without the will to prepare."
>> http://en.wikipedia.org/wiki/Juma_Ikangaa
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>