Report 3 - Box-Cox Transformation in Regression Analysis

if(!require(geoR)) install.packages("geoR")
if(!require(stats)) install.packages("stats")
if(!require(ggplot2)) install.packages("ggplot2")
if(!require(nortest)) install.packages("nortest")
if(!require(graphics)) install.packages("graphics")
if(!require(EnvStats)) install.packages("EnvStats")
if(!require(rcompanion)) install.packages("rcompanion")

Para realizar este estudo, utilizamos um dataset com 16.000 valores de entropia de permutação de Shannon e complexidade estatística referentes a ruídos brancos gerados artificialmente com amostras de tamanhos \(N \in \{10.000, 20.000, 30.000, 40.000, 50.000, 60.000, 70.000, 80.000, 90.000, 100.000\}\), aplicando todas as possíveis configurações de \(D \in \{3, 4, 5, 6\}\) e \(\tau \in \{1, 2, 3, 4\}\) aos descritores.

A disposição dos dados no plano \(HC\) pode ser observada abaixo:

HC.BP = data.frame("H" = numeric(16000), 
                   "C" = numeric(16000),
                   "Dist" = numeric(16000),
                   "D" = numeric(16000),
                   "t" = numeric(16000), 
                   "N" = numeric(16000), 
                   stringsAsFactors=FALSE)
HC.BP$N = as.factor(rep(c(rep(1e+04, 100), rep(2e+04, 100), rep(3e+04, 100), rep(4e+04, 100), rep(5e+04, 100), rep(6e+04, 100), rep(7e+04, 100), rep(8e+04, 100), rep(9e+04, 100), rep(1e+05, 100)), 16))
HC.BP$t = as.factor(rep(c(rep(1,1000), rep(2,1000), rep(3,1000), rep(4,1000)), 4))
file.csv = data.frame(read.csv("../Data/HC_series_fk0_16000.csv"))
HC.BP$H = file.csv[,1]
HC.BP$C = file.csv[,2]
HC.BP$Dist = HC.BP$C / HC.BP$H
HC.BP$D= as.factor(file.csv[,3])

Na geração dos modelos de regressão, consideramos três hipóteses de análise:

1. Os fatores \(D\) e \(N\) influenciam na construção de uma região de intervalo de confiança, logo são considerados na equação de regressão utilizando o operador (\(*\)).

2. Os fatores \(D\) e \(N\) influenciam na construção de uma região de intervalo de confiança e são considerados na equação de regressão usando o operador (\(+\)).

3. Consideraremos na equação de regressão apenas a variável \(H\), gerando diferentes modelos para as diferentes configurações de \(D\) e \(N\).

Hipótese 1 - Calculando a regressão como uma função de H, D e N (usando o operador *)

A seguir vamos fazer o ajuste do modelo e inspecionar os resíduos.

lm.alternative.1 = lm(data = HC.BP, formula = C ~ H * D * N)
summary(lm.alternative.1)

Call:
lm(formula = C ~ H * D * N, data = HC.BP)

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0045843 -0.0000738 -0.0000140  0.0000634  0.0035562 

Coefficients:
              Estimate Std. Error  t value Pr(>|t|)    
(Intercept)  0.9792309  0.0017126  571.771  < 2e-16 ***
H           -0.9792041  0.0017305 -565.859  < 2e-16 ***
D4          -0.0050475  0.0023441   -2.153 0.031311 *  
D5          -0.0134205  0.0023628   -5.680 1.37e-08 ***
D6          -0.0056669  0.0023251   -2.437 0.014808 *  
N20000      -0.0008542  0.0025504   -0.335 0.737686    
N30000      -0.0012288  0.0026234   -0.468 0.639523    
N40000       0.0012452  0.0025938    0.480 0.631171    
N50000      -0.0018326  0.0025542   -0.717 0.473083    
N60000      -0.0025412  0.0025032   -1.015 0.310032    
N70000      -0.0087140  0.0025459   -3.423 0.000621 ***
N80000      -0.0010086  0.0024629   -0.410 0.682158    
N90000      -0.0023638  0.0024987   -0.946 0.344145    
N1e+05       0.0027722  0.0024988    1.109 0.267275    
H:D4         0.0051042  0.0023680    2.155 0.031140 *  
H:D5         0.0135594  0.0023872    5.680 1.37e-08 ***
H:D6         0.0057087  0.0023491    2.430 0.015101 *  
H:N20000     0.0008580  0.0025759    0.333 0.739065    
H:N30000     0.0012488  0.0026495    0.471 0.637408    
H:N40000    -0.0012586  0.0026204   -0.480 0.631004    
H:N50000     0.0018453  0.0025802    0.715 0.474496    
H:N60000     0.0025688  0.0025288    1.016 0.309722    
H:N70000     0.0087846  0.0025716    3.416 0.000637 ***
H:N80000     0.0010270  0.0024880    0.413 0.679783    
H:N90000     0.0023836  0.0025239    0.944 0.344976    
H:N1e+05    -0.0027910  0.0025236   -1.106 0.268770    
D4:N20000    0.0009338  0.0035390    0.264 0.791888    
D5:N20000    0.0131977  0.0036219    3.644 0.000269 ***
D6:N20000    0.0027692  0.0034098    0.812 0.416718    
D4:N30000    0.0076083  0.0034998    2.174 0.029726 *  
D5:N30000    0.0143683  0.0036483    3.938 8.24e-05 ***
D6:N30000   -0.0078594  0.0035804   -2.195 0.028169 *  
D4:N40000    0.0015956  0.0035188    0.453 0.650233    
D5:N40000    0.0022151  0.0034557    0.641 0.521530    
D6:N40000    0.0035702  0.0035790    0.998 0.318509    
D4:N50000    0.0017496  0.0034519    0.507 0.612274    
D5:N50000    0.0143471  0.0036178    3.966 7.35e-05 ***
D6:N50000    0.0097167  0.0035058    2.772 0.005584 ** 
D4:N60000    0.0096611  0.0035758    2.702 0.006904 ** 
D5:N60000    0.0260962  0.0035571    7.336 2.30e-13 ***
D6:N60000    0.0027184  0.0034530    0.787 0.431145    
D4:N70000    0.0022586  0.0034393    0.657 0.511375    
D5:N70000    0.0214660  0.0034798    6.169 7.05e-10 ***
D6:N70000   -0.0019506  0.0033035   -0.590 0.554901    
D4:N80000    0.0111703  0.0034921    3.199 0.001383 ** 
D5:N80000    0.0128356  0.0034877    3.680 0.000234 ***
D6:N80000    0.0007265  0.0034909    0.208 0.835150    
D4:N90000   -0.0015940  0.0035375   -0.451 0.652278    
D5:N90000    0.0112703  0.0036562    3.082 0.002056 ** 
D6:N90000    0.0109826  0.0034801    3.156 0.001604 ** 
D4:N1e+05   -0.0059590  0.0035017   -1.702 0.088828 .  
D5:N1e+05    0.0075672  0.0034866    2.170 0.029993 *  
D6:N1e+05    0.0050481  0.0034998    1.442 0.149212    
H:D4:N20000 -0.0009663  0.0035739   -0.270 0.786866    
H:D5:N20000 -0.0133241  0.0036572   -3.643 0.000270 ***
H:D6:N20000 -0.0027803  0.0034443   -0.807 0.419546    
H:D4:N30000 -0.0076897  0.0035347   -2.175 0.029609 *  
H:D5:N30000 -0.0145001  0.0036845   -3.935 8.34e-05 ***
H:D6:N30000  0.0079203  0.0036157    2.191 0.028500 *  
H:D4:N40000 -0.0016186  0.0035546   -0.455 0.648859    
H:D5:N40000 -0.0022610  0.0034915   -0.648 0.517259    
H:D6:N40000 -0.0035788  0.0036150   -0.990 0.322197    
H:D4:N50000 -0.0017700  0.0034868   -0.508 0.611728    
H:D5:N50000 -0.0144810  0.0036539   -3.963 7.43e-05 ***
H:D6:N50000 -0.0097876  0.0035413   -2.764 0.005719 ** 
H:D4:N60000 -0.0097524  0.0036112   -2.701 0.006929 ** 
H:D5:N60000 -0.0263366  0.0035931   -7.330 2.42e-13 ***
H:D6:N60000 -0.0027489  0.0034884   -0.788 0.430692    
H:D4:N70000 -0.0022986  0.0034737   -0.662 0.508157    
H:D5:N70000 -0.0216719  0.0035152   -6.165 7.21e-10 ***
H:D6:N70000  0.0019819  0.0033369    0.594 0.552558    
H:D4:N80000 -0.0112826  0.0035278   -3.198 0.001386 ** 
H:D5:N80000 -0.0129847  0.0035230   -3.686 0.000229 ***
H:D6:N80000 -0.0007253  0.0035259   -0.206 0.837032    
H:D4:N90000  0.0015893  0.0035727    0.445 0.656444    
H:D5:N90000 -0.0113930  0.0036925   -3.085 0.002036 ** 
H:D6:N90000 -0.0110636  0.0035153   -3.147 0.001651 ** 
H:D4:N1e+05  0.0059638  0.0035362    1.687 0.091713 .  
H:D5:N1e+05 -0.0076538  0.0035213   -2.174 0.029754 *  
H:D6:N1e+05 -0.0050977  0.0035346   -1.442 0.149258    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.0002934 on 15920 degrees of freedom
Multiple R-squared:  0.9986,    Adjusted R-squared:  0.9986 
F-statistic: 1.465e+05 on 79 and 15920 DF,  p-value: < 2.2e-16

qqPlot(lm.alternative.1$residuals)

[1] 13535 11088

No gráfico do \(QQ-plot\) podemos ver que o comportamento dos dados se afasta muito de uma distribuição normal. Assim, realizamos testes para verificar o desvio dos pressupostos.

cat("Total: ", 1, " Teste de normalidade: ", length(which(ad.test(lm.alternative.1$residuals)$p.value > 0.05)), '\n')

Total:  1  Teste de normalidade:  0 

cat("P-value: ", ad.test(lm.alternative.1$residuals)$p.value, '\n')

P-value:  3.7e-24 

Hipótese 2 - Calculando a regressão como uma função de H, D e N (usando o operador +)

A seguir vamos fazer o ajuste do modelo e inspecionar os resíduos.

lm.alternative.2 = lm(data = HC.BP, formula = C ~ H + D + N)
summary(lm.alternative.2)

Call:
lm(formula = C ~ H + D + N, data = HC.BP)

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0046729 -0.0000665 -0.0000209  0.0000609  0.0037258 

Coefficients:
              Estimate Std. Error   t value Pr(>|t|)    
(Intercept)  9.757e-01  2.871e-04  3398.001   <2e-16 ***
H           -9.756e-01  2.899e-04 -3364.867   <2e-16 ***
D4          -5.609e-06  6.631e-06    -0.846    0.398    
D5           2.783e-06  6.631e-06     0.420    0.675    
D6          -4.425e-06  6.631e-06    -0.667    0.505    
N20000      -5.218e-06  1.049e-05    -0.498    0.619    
N30000       6.341e-06  1.049e-05     0.605    0.545    
N40000      -2.668e-06  1.049e-05    -0.254    0.799    
N50000       2.738e-06  1.049e-05     0.261    0.794    
N60000       6.757e-06  1.049e-05     0.644    0.519    
N70000      -1.466e-05  1.049e-05    -1.398    0.162    
N80000       6.549e-06  1.049e-05     0.625    0.532    
N90000      -3.320e-06  1.049e-05    -0.317    0.752    
N1e+05      -7.244e-06  1.049e-05    -0.691    0.490    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.0002965 on 15986 degrees of freedom
Multiple R-squared:  0.9986,    Adjusted R-squared:  0.9986 
F-statistic: 8.714e+05 on 13 and 15986 DF,  p-value: < 2.2e-16

qqPlot(lm.alternative.2$residuals)

[1] 13535 11088

No gráfico do \(QQ-plot\) podemos ver que o comportamento dos dados se afasta muito de uma distribuição normal. Assim, realizamos testes para verificar o desvio dos pressupostos.

cat("Total: ", 1, " Teste de normalidade: ", length(which(ad.test(lm.alternative.2$residuals)$p.value > 0.05)), '\n')

Total:  1  Teste de normalidade:  0 

cat("P-value: ", ad.test(lm.alternative.2$residuals)$p.value, '\n')

P-value:  3.7e-24 

Hipótese 3 - Calculando a regressão como uma função apenas de H

Levando em consideração na equação apenas a entropia de shannon também não conseguimos garantir a condição de normalidade dos resíduos, gerando assim intervalos de confiança sem precisão.

dimension = c(3,4,5,6)
N = c(1e+04, 2e+04, 3e+04, 4e+04, 5e+04, 6e+04, 7e+04, 8e+04, 9e+04, 1e+05)
HC.BP.regression.data = data.frame("H" = numeric(400), 
                                   "C" = numeric(400),
                                   stringsAsFactors=FALSE)
shapiro.alternative.3 = rep(0, 40)
lm.alternative.3 = array(list(), 40)
b = cc = 0
for(i in 1:length(dimension)){
  for(j in 1:length(N)){
    cc = cc + 1
    a = c((((j - 1) * 100) + 1):(j * 100))
    elements = c(a + b, a + b + 1000, a + b + 2000, a + b + 3000)
    HC.BP.regression.data$H = HC.BP$H[elements]
    HC.BP.regression.data$C = HC.BP$C[elements]
    lm.alternative.3[[cc]] = lm(data = HC.BP.regression.data, formula = C ~ H)
    shapiro.alternative.3[cc] = shapiro.test(lm.alternative.3[[cc]]$residuals)$p.value
  }
  b = b + 4000
}
cat("Total: ", 40, " Teste de normalidade: ", length(which(shapiro.alternative.3 > 0.05)), '\n')

Total:  40  Teste de normalidade:  0 

qqPlot(lm.alternative.3[[1]]$residuals)

[1] 162 325

Logo, em casos como o aqui apresentado (resíduos não normais) a literatura propõe duas soluções:

1. A aplicação de transformações sobre os dados ou,

2. A escolha de um modelo de regressão mais adequado.

Como a tranformação de dados é uma das possíveis formas de contarnar o problema de séries que não obedecem os pressupostos da análise de variância e normalidade, neste relatório optamos por testar o impacto da transformação de Box Cox e avaliar os resíduos gerados pelo método de regressão linear.

Transformação de Box-Cox

Hipótese 1 - Calculando a regressão como uma função de H, D e N (usando o operador *)

Quando analisamos os valores de complexidade estatística do dataset isoladamente conseguimos verificar que estes não seguem uma distribuição normal:

plotNormalHistogram(HC.BP$C)

qqPlot(HC.BP$C)

[1] 12601   740

Para tentar contornar o problema vamos usar a transformação Box-Cox, que consiste em transformar os dados de acordo com a expressão:

\[ y(\lambda) = \left\{ \begin{eqnarray*} \frac{y^{\lambda} - 1}{\lambda}, \lambda \neq 0\\ \\ log(y), \lambda = 0.\\ \end{eqnarray*} \right. \]

Logo, aplicamos a transformação sobre a variável resposta do nosso problema, ou seja, nos valores de complexidade estatística do dataset de ruídos branco.

parameters = boxcoxfit(HC.BP$C, lambda2 = F)
lambda = parameters$lambda[1]
if(lambda==0){T.norm = log(HC.BP$C)}
if(lambda!=0){T.norm = ((HC.BP$C)^lambda - 1)/lambda}
plotNormalHistogram(T.norm)

Após aplicar a transformação nossos dados passam a se apresentar dispostos do seguinte modo:

plot(x = HC.BP$H, y = T.norm)

E obtemos o seguinte \(QQ-plot\):

qqPlot(T.norm)

[1] 12601   740

Model fitting (After transformation)

HC.BP.regression.data = data.frame("H" = numeric(16000), 
                                   "C" = numeric(16000),
                                   "D" = numeric(16000),
                                   "N" = numeric(16000), 
                                   stringsAsFactors=FALSE)
HC.BP.regression.data$H = HC.BP$H
HC.BP.regression.data$C = T.norm
HC.BP.regression.data$D = HC.BP$D
HC.BP.regression.data$N = HC.BP$N
  
lm.1.box.cox = lm(data = HC.BP.regression.data, formula = C ~ H * D * N)
ggplot(HC.BP.regression.data, aes(x = H, y = C)) +
  geom_point() +
  scale_fill_grey() +
  geom_line(aes(y = predict(lm.1.box.cox, HC.BP.regression.data, type = 'response'))) 

cat("Total: ", 1, " Teste de normalidade: ", length(which(ad.test(lm.1.box.cox$residuals)$p.value > 0.05)), '\n')

Total:  1  Teste de normalidade:  0 

cat("P-value: ", ad.test(lm.1.box.cox$residuals)$p.value, '\n')

P-value:  3.7e-24 

Analisando o plot do modelo:

plot(lm.1.box.cox)

Warning messages:
1: In dir.create(dirname(dest), recursive = TRUE) :
  cannot create dir '/home/eduarda/../Data', reason 'Permission denied'
2: In file.create(to[okay]) :
  cannot create file '/home/eduarda/../Data/HC_series_fk0_16000.csv', reason 'No such file or directory'

Examinando o histograma dos resíduos:

hist(lm.1.box.cox$residuals, breaks = 100)

Examinando o \(QQ-plot\) dos resíduos:

qqPlot(rstandard(lm.1.box.cox)) 

[1]   740 11088

Como podemos observar abaixo, com a transformação Box Cox obtemos um menor ajuste R-quadrado:

adj.rsq1 = summary(lm.alternative.1)$adj.r.squared
adj.rsq2 = summary(lm.1.box.cox)$adj.r.squared
cat("Adjusted R-Square (before transformation):", adj.rsq1, "Adjusted R-Square (after transformation):", adj.rsq2, sep = "\n")

Adjusted R-Square (before transformation):
0.9986192
Adjusted R-Square (after transformation):
0.8394782

Hipótese 2 - Calculando a regressão como uma função de H, D e N (usando o operador +)

  
lm.2.box.cox = lm(data = HC.BP.regression.data, formula = C ~ H + D + N)
ggplot(HC.BP.regression.data, aes(x = H, y = C)) +
  geom_point() +
  scale_fill_grey() +
  geom_line(aes(y = predict(lm.2.box.cox, HC.BP.regression.data, type = 'response'))) 

cat("Total: ", 1, " Teste de normalidade: ", length(which(ad.test(lm.2.box.cox$residuals)$p.value > 0.05)), '\n')

Total:  1  Teste de normalidade:  0 

cat("P-value: ", ad.test(lm.2.box.cox$residuals)$p.value, '\n')

P-value:  3.7e-24 

Analisando o plot do modelo:

plot(lm.2.box.cox)

Examinando o histograma dos resíduos:

hist(lm.2.box.cox$residuals, breaks = 100)

Examinando o \(QQ-plot\) dos resíduos:

qqPlot(rstandard(lm.2.box.cox)) 

[1] 12601   740

Como podemos observar abaixo, com a transformação Box Cox obtemos um menor ajuste R-quadrado:

adj.rsq1 = summary(lm.alternative.2)$adj.r.squared
adj.rsq2 = summary(lm.2.box.cox)$adj.r.squared
cat("Adjusted R-Square (before transformation):", adj.rsq1, "Adjusted R-Square (after transformation):", adj.rsq2, sep = "\n")

Adjusted R-Square (before transformation):
0.9985897
Adjusted R-Square (after transformation):
0.8356933

Hipótese 3 - Calculando a regressão como uma função apenas de H

dimension = c(3,4,5,6)
N = c(1e+04, 2e+04, 3e+04, 4e+04, 5e+04, 6e+04, 7e+04, 8e+04, 9e+04, 1e+05)
HC.test.3 = data.frame("H" = numeric(400), 
                       "C" = numeric(400),
                       stringsAsFactors=FALSE)
shapiro.3.box.cox = rep(0, 40)
lm.3.box.cox = array(list(), 40)
b = cc = 0
for(i in 1:length(dimension)){
  for(j in 1:length(N)){
    cc = cc + 1
    a = c((((j - 1) * 100) + 1):(j * 100))
    elements = c(a + b, a + b + 1000, a + b + 2000, a + b + 3000)
    
    HC.test.3$H = HC.BP.regression.data$H[elements]
    HC.test.3$C = HC.BP.regression.data$C[elements]
    
    lm.3.box.cox[[cc]] = lm(data = HC.test.3, formula = C ~ H)
    shapiro.3.box.cox[cc] = shapiro.test(lm.3.box.cox[[cc]]$residuals)$p.value
  }
  b = b + 4000
}
cat("Total: ", 40, " Teste de normalidade: ", length(which(shapiro.3.box.cox > 0.05)), '\n')

Total:  40  Teste de normalidade:  0 

HC.test.3$H = HC.BP.regression.data$H[1:100]
HC.test.3$C = HC.BP.regression.data$C[1:100]
plot(x = HC.test.3$H, y = HC.test.3$C)

ggplot(HC.test.3, aes(x = H, y = C)) +
  geom_point() +
  scale_fill_grey() +
  geom_line(aes(y = predict(lm.3.box.cox[[1]], HC.test.3, type = 'response'))) 

Analisando o plot do modelo:

plot(lm.3.box.cox[[1]])

Warning messages:
1: In dir.create(dirname(dest), recursive = TRUE) :
  cannot create dir '/home/eduarda/../Data', reason 'Permission denied'
2: In file.create(to[okay]) :
  cannot create file '/home/eduarda/../Data/HC_series_fk0_16000.csv', reason 'No such file or directory'

hist(lm.3.box.cox[[1]]$residuals, breaks = 200)

Examinando o \(QQ-plot\) dos resíduos:

qqPlot(rstandard(lm.3.box.cox[[1]])) 

[1] 100 336

Como podemos observar abaixo, além de não alcançar a normalidade dos resíduos a regressão feita com a transformação Box Cox obteve um menor ajuste R-quadrado:

adj.rsq1 = summary(lm.alternative.3[[1]])$adj.r.squared
adj.rsq2 = summary(lm.3.box.cox[[1]])$adj.r.squared
cat("Adjusted R-Square (before transformation):", adj.rsq1, "Adjusted R-Square (after transformation):", adj.rsq2, sep = "\n")

Adjusted R-Square (before transformation):
0.9985105
Adjusted R-Square (after transformation):
0.837164