# Univariate Example from Hawkins and Weisberg (2015)
m1 <- lm(I1L1 ~ pool, LoBD)
p1 <- powerTransform(m1, family="skewPower")
summary(p1)
# summary prints estimate, se and conf.ints for both parameters
# helper functions
coef(p1)
vcov(p1)
p1$value # maximum value of log-likelihood surface
# 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))
# Create new variate corresponding to the transformation, then update the regression
LoBD$z <- skewPower(LoBD$I1L1, lambda=p1$lambda, gamma=p1$gamma, jacobian.adjusted=FALSE)
m1.t <- update(m1, z ~ . -1)
# The example in Hawkins and Weisberg(2015) requires computing the inverse of the skew power
# transformation and certain summaries of the updated regression
mu0 <- mean(coef(m1.t)[1:4])
sd <- sigmaHat(m1.t)
LoB <- mu0 + 1.645*sd
LoD <- LoB + 1.645*sd
skewPowerInverse <- function(z, lambda, gamma){
q <- if(abs(lambda) <= 1.e-9) {2 * log(z)} else {2*(z*lambda + 1) ^(1/lambda)}
(q^2 - gamma^2)/(2 * q)
}
y.untransformed <- round(skewPowerInverse(c(LoB, LoD), p1$lambda, p1$gamma), 2)
# Multivariate Response
p3 <- powerTransform(update(m1, as.matrix(cbind(LoBD$I1L2, LoBD$I1L1)) ~ .), family="skewPower")
summary(p3)
# Construct test for the same transformations for each response. Get new data file by stacking
d2 <- data.frame(assay = c(LoBD$I1L2, LoBD$I1L1), pool=factor(c(LoBD$pool, LoBD$pool)),
set=factor(rep(c("A", "B"), c(84, 84))))
# the following allows for different within set means, but common transformation parameters
p5 <- powerTransform(lm(assay ~ pool*set, d2), family="skewPower") # allow grps to differ by set.
c(df=(df <- length(coef(p3)) - length(coef(p5))),
test = (test <- -.5 * (p5$value - p3$value)),
pvalue = 1 - pchisq(test, df))
Run the code above in your browser using DataCamp Workspace