Report 4 - Beta regression using Box-Cox Transformation

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(betareg)) install.packages("betareg")
if(!require(glmmTMB)) install.packages("glmmTMB")
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])

Transformação de Box-Cox

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. Após aplicar a transformação nossos dados passam a se apresentar dispostos do seguinte modo:

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}
T.norm = (T.norm - min(T.norm))/(max(T.norm) - min(T.norm))
n.obs = sum(!is.na(T.norm))
T.norm = (T.norm * (n.obs - 1) + 0.5)/n.obs
plot(x = HC.BP$H, y = T.norm)

Regressão beta

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 *)

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
r.squared.1 = data.frame("logit" = numeric(1), "probit" = numeric(1), "cauchy" = numeric(1), "cloglog" = numeric(1),"loglog" = numeric(1),    
                         stringsAsFactors=FALSE)
lm.fit.logit.1 = betareg(data = HC.BP.regression.data, formula = C ~ H * D * N, link = "logit")
r.squared.1$logit = lm.fit.logit.1$pseudo.r.squared
  
lm.fit.probit.1 = betareg(data = HC.BP.regression.data, formula = C ~ H * D * N, link = "probit")
r.squared.1$probit = lm.fit.probit.1$pseudo.r.squared
lm.fit.cauchit.1 = betareg(data = HC.BP.regression.data, formula = C ~ H * D * N, link = "cauchit")
r.squared.1$cauchy = lm.fit.cauchit.1$pseudo.r.squared
  
lm.fit.cloglog.1 = betareg(data = HC.BP.regression.data, formula = C ~ H * D * N, link = "cloglog")
r.squared.1$cloglog = lm.fit.cloglog.1$pseudo.r.squared
  
lm.fit.loglog.1 = betareg(data = HC.BP.regression.data, formula = C ~ H * D * N, link = "loglog")
r.squared.1$loglog = lm.fit.loglog.1$pseudo.r.squared
summary(r.squared.1)

     logit            probit           cauchy            cloglog           loglog      
 Min.   :0.5379   Min.   :0.6963   Min.   :0.007184   Min.   :0.3925   Min.   :0.8519  
 1st Qu.:0.5379   1st Qu.:0.6963   1st Qu.:0.007184   1st Qu.:0.3925   1st Qu.:0.8519  
 Median :0.5379   Median :0.6963   Median :0.007184   Median :0.3925   Median :0.8519  
 Mean   :0.5379   Mean   :0.6963   Mean   :0.007184   Mean   :0.3925   Mean   :0.8519  
 3rd Qu.:0.5379   3rd Qu.:0.6963   3rd Qu.:0.007184   3rd Qu.:0.3925   3rd Qu.:0.8519  
 Max.   :0.5379   Max.   :0.6963   Max.   :0.007184   Max.   :0.3925   Max.   :0.8519  

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

r.squared.2 = data.frame("logit" = numeric(1), "probit" = numeric(1), "cauchy" = numeric(1), "loglog" = numeric(1),    
                         stringsAsFactors=FALSE)
lm.fit.logit.2 = betareg(data = HC.BP.regression.data, formula = C ~ H + D + N, link = "logit")
r.squared.2$logit = lm.fit.logit.2$pseudo.r.squared
lm.fit.cauchit.2 = betareg(data = HC.BP.regression.data, formula = C ~ H + D + N, link = "cauchit")
r.squared.2$cauchy = lm.fit.cauchit.2$pseudo.r.squared
  
lm.fit.probit.2 = betareg(data = HC.BP.regression.data, formula = C ~ H + D + N, link = "probit")
r.squared.2$probit = lm.fit.probit.2$pseudo.r.squared
  
lm.fit.loglog.2 = betareg(data = HC.BP.regression.data, formula = C ~ H + D + N, link = "loglog")
r.squared.2$loglog = lm.fit.loglog.2$pseudo.r.squared
summary(r.squared.2) 

     logit            probit           cauchy             loglog      
 Min.   :0.5351   Min.   :0.6941   Min.   :0.006429   Min.   :0.8531  
 1st Qu.:0.5351   1st Qu.:0.6941   1st Qu.:0.006429   1st Qu.:0.8531  
 Median :0.5351   Median :0.6941   Median :0.006429   Median :0.8531  
 Mean   :0.5351   Mean   :0.6941   Mean   :0.006429   Mean   :0.8531  
 3rd Qu.:0.5351   3rd Qu.:0.6941   3rd Qu.:0.006429   3rd Qu.:0.8531  
 Max.   :0.5351   Max.   :0.6941   Max.   :0.006429   Max.   :0.8531  

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)
lm.fit.logit.3 = lm.fit.probit.3 = lm.fit.cloglog.3 = lm.fit.loglog.3 = array(list(), 40)
r.squared.3 = data.frame("logit" = numeric(40), "probit" = numeric(40), "cloglog" = numeric(40),"loglog" = numeric(40), stringsAsFactors=FALSE)
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.fit.logit.3[[cc]] = betareg(data = HC.test.3, formula = C ~ H, link = "logit")
      r.squared.3$logit[cc] = lm.fit.logit.3[[cc]]$pseudo.r.squared
      
      lm.fit.probit.3[[cc]] = betareg(data = HC.test.3, formula = C ~ H, link = "probit")
      r.squared.3$probit[cc] = lm.fit.probit.3[[cc]]$pseudo.r.squared
      
      lm.fit.cloglog.3[[cc]] = betareg(data = HC.test.3, formula = C ~ H, link = "cloglog")
      r.squared.3$cloglog[cc] = lm.fit.cloglog.3[[cc]]$pseudo.r.squared
      
      lm.fit.loglog.3[[cc]] = betareg(data = HC.test.3, formula = C ~ H, link = "loglog")
      r.squared.3$loglog[cc] = lm.fit.loglog.3[[cc]]$pseudo.r.squared
  } 
  b = b + 4000
}
summary(r.squared.3)

     logit            probit          cloglog           loglog      
 Min.   :0.3506   Min.   :0.5413   Min.   :0.2428   Min.   :0.7712  
 1st Qu.:0.4743   1st Qu.:0.6551   1st Qu.:0.3450   1st Qu.:0.8328  
 Median :0.5841   Median :0.7184   Median :0.4362   Median :0.8669  
 Mean   :0.5935   Mean   :0.7144   Mean   :0.4688   Mean   :0.8596  
 3rd Qu.:0.6948   3rd Qu.:0.7777   3rd Qu.:0.5271   3rd Qu.:0.8841  
 Max.   :0.8642   Max.   :0.8684   Max.   :0.7993   Max.   :0.9279  

Análise dos plots dos modelos obtidos

Como visto acima, os melhores ajustes são obtidos ao aplicar a função de link loglog sobre o modelo de regressão. Logo, analisaremos os gráficos de diagnótico para tais regressões.

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

plot(lm.fit.loglog.1) 

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

plot(lm.fit.loglog.2) 

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

plot(lm.fit.loglog.3[[1]]) 

Visualizando os ajustes obtidos pela função de link loglog

ggplot(HC.BP.regression.data, aes(x = H, y = C)) +
  geom_point() +
  scale_fill_grey() +
  geom_line(aes(y = predict(lm.fit.loglog.1, HC.BP.regression.data, type = 'response'))) 

ggplot(HC.BP.regression.data, aes(x = H, y = C)) +
  geom_point() +
  scale_fill_grey() +
  geom_line(aes(y = predict(lm.fit.loglog.2, HC.BP.regression.data, type = 'response'))) 

a = c(1:100)
elements = c(a, a + 1000, a + 2000, a + 3000)
HC.test.3$H = HC.BP.regression.data$H[elements]
HC.test.3$C = HC.BP.regression.data$C[elements]
      
ggplot(HC.test.3, aes(x = H, y = C)) +
  geom_point() +
  scale_fill_grey() +
  geom_line(aes(y = predict(lm.fit.cloglog.3[[1]], HC.test.3, type = 'response'))) 

Ajuste R-quadrado de cada hipótese analisada

adj.rsq1 = lm.fit.loglog.1$pseudo.r.squared
adj.rsq2 = lm.fit.loglog.2$pseudo.r.squared
adj.rsq3 = lm.fit.loglog.3[[1]]$pseudo.r.squared
cat(" Adjusted R-Square (hypothesis 1):", adj.rsq1, "\n", 
    "Adjusted R-Square (hypothesis 2):", adj.rsq2, "\n", 
    "Adjusted R-Square (hypothesis 3):", adj.rsq3, "\n")

 Adjusted R-Square (hypothesis 1): 0.8519265 
 Adjusted R-Square (hypothesis 2): 0.8530805 
 Adjusted R-Square (hypothesis 3): 0.8407004