BICadj: Ways to compare SITAR models for fit

Description

BICadj and AICadj calculate the BIC and AIC for SITAR models, adjusting the likelihood for Box-Cox transformed y variables. varexp calculates the variance explained by SITAR models, compared to the corresponding fixed effect models. getL is used by [AB]ICadj to find what power the y variable is raised to.

Usage

BICadj(..., pattern = NULL)
AICadj(..., k = 2, pattern = NULL)
varexp(..., pattern = NULL)
getL(expr)

Arguments

…

one or more SITAR models.

pattern

regular expression defining names of models.

numeric, the penalty per parameter to be used; the default k = 2 is the classical AIC.

expr

quoted or unquoted expression containing a single variable name.

Value

For BICadj and AICadj a named vector of deviances in increasing order. For varexp a named vector of percentages in decreasing order. For getL the power the variable in expr is raised to, or NA if expr is not a power of (a multiple of) the variable.

Details

The deviance is adjusted if the y variable is power-transformed, using the formula $$adjusted deviance = deviance - 2n ( (\lambda-1) * log(gm) + % log(abs(\lambda)) )$$ where $\lambda$ is the power transform, and $n$ and $gm$ are the length and geometric mean of y.

The variance explained is given by $$\% explained = 100 * (1 -% (\sigma_2/\sigma_1)^2)$$ where $\sigma_1$ is the fixed effects RSD and $\sigma_2$ the SITAR random effects RSD.

BICadj and AICadj accept non-sitar models with a logLik class. varexp ignores objects not of class sitar.

getL does not detect if the variable in expr, or its log, contains a multiplying constant, so that the expressions log(x) and 1 + 2 * log(3 * x) both return 0.

Examples

Run this code

# NOT RUN {
data(heights)
## fit sitar model for height
m1 <- sitar(x=age, y=height, id=id, data=heights, df=5)

## update it for log(height)
m2 <- update(m1, y=sqrt(height))

## compare variance explained in the two models
varexp(m1, m2)

## compare BIC adjusting for sqrt transform
## the pattern matches names starting with "m" followed by a digit
BICadj(pattern="^m[0-9]")

## find what power height is raised to
getL(quote(sqrt(sqrt(height))))
# }