R Notebook

# ref1 : https://www.isixsigma.com/tools-templates/normality/making-data-normal-using-box-cox-power-transformation/
library(MASS)
par(mfrow=c(3,2))
plot(lm(dist~speed,data=cars))
boxcox(lm(dist~speed,data=cars),lambda=seq(0,1,by=.1))

library(Matrix)
library(car)
library(stats)
library(MASS)
library(graphics)

# successful case
par(mfrow=c(3,2))
x <- rexp(100)
qqnorm(x)                 # Normal probability plot for original variable
boxcox(x~1)               # Illustration of Log-Likelihood profile

p <- box.cox.powers(x)    # Estimaton of Box-Cox lambda

'box.cox.powers' is deprecated.
Use 'powerTransform' instead.
See help("Deprecated") and help("car-deprecated").

y <- box.cox(x, p$lambda) # Box-Cox transformation

'box.cox' is deprecated.
Use 'bcPower' instead.
See help("Deprecated") and help("car-deprecated").

car::qqPlot(x); car::qqPlot(y) 
plot(density(x)); plot(density(y))

# unsuccessful case
par(mfrow=c(3,2))
x<-rnorm(1000)^3        
x<-x-min(x)+sd(x)
qqnorm(x)                     # Normal probability plot for original variable
boxcox(x~1)             # Illustration of Log-Likelihood profile

p <- box.cox.powers(x)    # Estimaton of Box-Cox lambda

'box.cox.powers' is deprecated.
Use 'powerTransform' instead.
See help("Deprecated") and help("car-deprecated").

y <- box.cox(x, p$lambda)   # Box-Cox transformation

'box.cox' is deprecated.
Use 'bcPower' instead.
See help("Deprecated") and help("car-deprecated").

car::qqPlot(x); car::qqPlot(y) 
plot(density(x)); plot(density(y))

###################################################################################################################
# transforming skewed data                                                                                        #
# ref1 : http://rcompanion.org/handbook/I_12.html                                                                 #
# ref2 : http://fmwww.bc.edu/repec/bocode/t/transint.html                                                         #
###################################################################################################################
#install.packages('rcompanion')
par(mfrow=c(4,4))
library(rcompanion)
Turbidity = c(1.0, 1.2, 1.1, 1.1, 2.4, 2.2, 2.6, 4.1, 5.0, 10.0, 4.0, 4.1, 4.2, 4.1, 5.1, 4.5, 5.0, 15.2, 10.0, 20.0, 1.1, 1.1, 1.2, 1.6, 2.2, 3.0, 4.0, 10.5)
plotNormalHistogram(Turbidity, main = "original data")
qqnorm(Turbidity, ylab="Sample Quantiles for Turbidity")

qqline(Turbidity, col="red")
# some transformations
T_sqrt = sqrt(Turbidity)
plotNormalHistogram(T_sqrt, main ="square root transformation")

T_cub = sign(Turbidity) * abs(Turbidity)^(1/3)
plotNormalHistogram(T_cub, main ="cube root transformation")
T_log = log(Turbidity)
plotNormalHistogram(T_log, main ="log transformation")

T_tuk = transformTukey(Turbidity, plotit=FALSE)

    lambda     W Shapiro.p.value
397   -0.1 0.935         0.08248

if (lambda >  0){TRANS = x ^ lambda} 
if (lambda == 0){TRANS = log(x)} 
if (lambda <  0){TRANS = -1 * x ^ lambda} 

plotNormalHistogram(T_tuk, main = "tukey's ladder of power transformation")
Box = boxcox(Turbidity ~ 1, lambda = seq(-6,6,0.1))

Cox = data.frame(Box$x, Box$y)            
Cox2 = Cox[with(Cox, order(-Cox$Box.y)),] 
Cox2[1,]                                  

lambda = Cox2[1, "Box.x"]                 # Extract that lambda
T_box = (Turbidity ^ lambda - 1)/lambda   # Transform the original data
plotNormalHistogram(T_box, main = "box-cox transformation")
# box-cox T for anova model 
Input =("
        Location Turbidity
        a        1.0
        a        1.2
        a        1.1
        a        1.1
        a        2.4
        a        2.2
        a        2.6
        a        4.1
        a        5.0
        a       10.0
        b        4.0
        b        4.1
        b        4.2
        b        4.1
        b        5.1
        b        4.5
        b        5.0
        b       15.2
        b       10.0
        b       20.0
        c        1.1
        c        1.1
        c        1.2
        c        1.6
        c        2.2
        c        3.0
        c        4.0
        c       10.5
        ")
Data = read.table(textConnection(Input),header=TRUE)
model = lm(Turbidity ~ Location, data=Data)
Anova(model, type = "II")

Anova Table (Type II tests)

Response: Turbidity
          Sum Sq Df F value  Pr(>F)  
Location  132.63  2  3.8651 0.03447 *
Residuals 428.95 25                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

x <- residuals(model)
plotNormalHistogram(x, main = "residuals of un-transformed data")

qqnorm(residuals(model), ylab="Sample Quantiles for residuals")
qqline(residuals(model), col="red")

plot(fitted(model),residuals(model))
Box = boxcox(Turbidity ~ Location,
             data = Data,
             lambda = seq(-6,6,0.1)
)

Cox = data.frame(Box$x, Box$y)
Cox2 = Cox[with(Cox, order(-Cox$Box.y)),]
Cox2[1,]

lambda = Cox2[1, "Box.x"]
Data$Turbidity_box = (Data$Turbidity ^ lambda - 1)/lambda   
boxplot(Turbidity_box ~ Location,
        data = Data,
        ylab="Box–Cox-transformed Turbidity",
        xlab="Location")
model = lm(Turbidity_box ~ Location, data=Data)
Anova(model, type="II")

Anova Table (Type II tests)

Response: Turbidity_box
          Sum Sq Df F value Pr(>F)   
Location  4.1643  2  6.6929 0.0047 **
Residuals 7.7774 25                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

x = residuals(model)
plotNormalHistogram(x, main = "box-cox T model")

qqnorm(residuals(model), ylab="Sample Quantiles for residuals")
qqline(residuals(model), col="red")

plot(fitted(model),residuals(model))