# Univariate Example, from Hawkins and Weisberg (2016)
m1 <- lm(I1L1 ~ pool, LoBD)
p1 <- powerTransform(m1, family="skewPower")
summary(p1)
# summary prints estimate, se and conf.ints for both parameters
# helper functions
c(p1$lambda, p1$gamma, LogLik=p1$llik)
vcov(p1) # Estimated covarinace from inverse Hessian
# tests are for lambda, maximizing over gamma (profile log-likelihoods
testTransform(p1, lambda=0.5)
# Contour plot of the log-likelihood
contour(p1, main="", levels=c(.5, .95, .99))
# the boxCox function can provide profile log-likelihoods for each of the two parameters:
boxCox(m1, family="skewPower", param="lambda", lambda=seq(0.25, 1.1, length=100))
boxCox(m1, family="skewPower", param="gamma", gamma=seq(3, 80, length=100))
# Fit with fixed gamma, in this case fixed at the estimate from p1
p2 <- powerTransform(m1, family="skewPower", gamma=p1$gamma)
# summary gives different tests because gamma is fixed rather than maximized over
summary(p2)
# Multivariate Response
p3 <- powerTransform(update(m1, as.matrix(cbind(LoBD$I1L2, LoBD$I1L1)) ~ .),
family="skewPower")
summary(p3)
# gamma fixed
p4 <- powerTransform(update(m1, as.matrix(cbind(LoBD$I1L2, LoBD$I1L1)) ~ .),
family="skewPower", gamma=p3$gamma)
summary(p4)
# mixed models fit with lmer - requires lmer and nloptr packages
data <- reshape(LoBD[1:20,], varying=names(LoBD)[-1], direction="long", v.names="y")
names(data) <- c("pool", "assay", "y", "id")
data$assay <- factor(data$assay)
require(lme4)
m2 <- lmer(y ~ pool + (1|assay), data)
f1 <- powerTransform(m2, family="skewPower")
summary(f1)
f2 <- powerTransform(m2, family="skewPower", gamma= 10)
summary(f2)
Run the code above in your browser using DataLab