Report 2 - Beta regression analysis

if(!require(fitdistrplus)) install.packages("fitdistrplus")
if(!require(EnvStats)) install.packages("EnvStats")
if(!require(betareg)) install.packages("betareg")
if(!require(ggplot2)) install.packages("ggplot2")
if(!require(car)) install.packages("car")

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

O primeiro passo que deve ser realizado ao modelar um problema por meio de uma regressão linear é analisar o comportamento das variáveis resposta e explicativa, verificando quais distribuições tais dados assumem. Logo, para isso analisaremos separadamente o histograma das variáveis consideradas no problema (entropia de Shannon e complexidade estatística).

Na geração dos modelos de regressão, vamos considerar duas hipóteses de análise: na primeira, apenas \(D\) e \(\tau\) influenciam na construção de uma região de intervalo de confiança, e na segunda além de tais parâmetros, o tamanho \(N\) das séries consideradas também possui um importante fator de impacto nesta construção.

Hypothesis 1: D and N matter

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.data.1 = list(data.frame("H" = numeric(400), "C" = numeric(400), stringsAsFactors=FALSE), 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.data.1$H[[cc]] = HC.BP$H[elements]
    HC.data.1$C[[cc]] = HC.BP$C[elements]
  }
  b = b + 4000
}

Facilitando a visualização, mostraremos neste relatório apenas o comportamento dos dados de 1 dos 16 subconjuntos utilizados.

hist(HC.data.1$H[[1]], breaks = 10)

hist(HC.data.1$C[[1]], breaks = 10)

Observando o comportamento das distribuições, podemos inferir que \(H_1 \sim \mathbb{\beta}(\mu_{11}, \phi_{11})\) e \(C_1 \sim \mathbb{\beta}(\mu_{12}, \phi_{12})\). Para testar tais suposições, aplicamos o fit sobre os dados e obtemos os resuldados abaixo.

fit.H.1 = fitdist(HC.data.1$H[[1]], "beta")

NaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs produced

plot(fit.H.1, las = 1) 

fit.C.1 = fitdist(HC.data.1$C[[1]], "beta")

NaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs produced

plot(fit.C.1, las = 1) 

Hypothesis 2: D, t and N matter

Seguindo os mesmos passos com a hipótese 2:

HC.data.2 = array(data.frame("H" = numeric(100), "C" = numeric(100), stringsAsFactors=FALSE), 160)
for(i in 1:160){
  elements = c((((i-1)*100)+1):(i*100))
  HC.data.2$H[[i]] = HC.BP$H[elements]
  HC.data.2$C[[i]] = HC.BP$C[elements]
}
hist(HC.data.2$H[[1]], breaks = 10)

hist(HC.data.2$C[[1]], breaks = 10)

fit.H.2 = fitdist(HC.data.2$C[[1]], "beta")

NaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs produced

plot(fit.H.2, las = 1) 

fit.C.2 = fitdist(HC.data.2$C[[1]], "beta")

NaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs produced

plot(fit.C.2, las = 1) 

Regressão beta

Os parâmetros desconhecidos de uma regressão linear são geralmente determinados pelo método da máxima verossimilhança, assumindo que dada as variáveis explicativas, a variável resposta correspondente segue uma distribuição gaussiana. No entando, verificamos acima que as variáveis do problema em questão assumem valores no intervalo unitário padrão \((0, 1)\) e seguem a distribuição beta. Tornando assim, as aproximações baseadas em distribuições gaussianas para estimativa de intervalos e teste de hipóteses bastante imprecisas.

Levando em consideração os fatores acima, analisamos a eficiência do modelo de regressão beta. Tal modelo inclui uma estrutura de variação mais flexível determinada por variáveis independentes distribuídas beta que estão relacionadas a um conjunto de variáveis independentes por meio de uma estrutura de regressão. É naturalmente heterocedástico, acomoda facilmente assimetrias e os seus parâmetros de regressão são interpretáveis em termos da média da variável de interesse.

As densidades beta podem ser expressas em uma parametrização em termos da média \(\mu\) e do parâmetro de precisão \(\phi\) do seguinte modo:

\[ f(y|\mu,\phi) = \frac{\Gamma(\phi)}{\Gamma(\mu\phi)\Gamma((1-\mu)\phi)} y^{\mu\phi-1} (1-y)^{(1-\mu)\phi-1}, \] sendo \(0 < y < 1\), \(0 < \mu < 1\) e \(\phi > 0\). Logo, temos que:

\[ \mathbb{E}(y) = \mu, \]

\[ \mathbb{Var}(y) = \frac{\mu(1-\mu)}{(1 + \phi)} , \]

O parâmetro \(\phi\) é conhecido como parâmetro de precisão, pois, para \(\mu\) fixo, quanto maior \(\phi\), menor a variação de y. \(\phi^{-1}\) é um parâmetro de dispersão.

Seja \(y_1, y_2, \dots, y_n\) um conjunto de variáveis aleatórias, tais que \(y_i \sim \mathbb{\beta}(\mu_i, \phi)\), \(i = 1, 2, \dots, n\). O modelo de regressão beta é definido como:

\[ g(\mu_i) = x_i^T \alpha = \eta_i, \]

onde \(\alpha = (\alpha_1, \dots, \alpha_n)^T\) é o vetor dos parâmetros desconhecidos da regressão e possui dimensão \(k \times 1\) \((k < n)\), \(x_i = (x_{i1}, \dots, x_{ik})^T\) é o vetor de \(k\) variáveis independentes e \(\eta_i\) é a variável linear preditora (isto é, \(\eta_i = \alpha_1 x_{i1} + \dots + \alpha_n x_{in}\), usualmente \(x_{i1} = 1\) para todo \(i\)).

A função \(g(.) : (0,1) \mapsto \mathbb{R}\) é uma função link, que é estritamente crescente e duplamente diferenciável. Existem duas principais motivações para o uso de uma função link na estrutura de regressão. Primeiramente, ambos os lados da equação de regressão assumem valores na linha real quando uma função de link é aplicada a \(\mu_i\). Em segundo lugar, há uma flexibilidade adicional, pois podemos escolher a função que produz o melhor ajuste.

As principais funções de link utilizadas são:

1. Logit: \(g(\mu) = log(\frac{\mu}{(1-\mu)})\);

2. Probit: \(g(\mu) = \phi^{-1}(\mu)\), onde \(\phi\) é a função de distribuição normal padrão;

3. Complementary log-log: \(g(\mu) = log(-log(1-\mu))\);

4. Log-log: \(g(\mu) = log(-log(\mu))\);

5. Cauchy: \(g(\mu) = tan(\pi(\mu - 0.5))\).

Seleção entre diferentes funções de link

A seleção de uma função de link apropriada pode melhorar muito o ajuste do modelo, especialmente se proporções extremas (próximas a 0 ou 1) foram observadas nos dados, como ocorre o caso em questão. Desse modo, modelamos nossas hipóteses com as principais variações da função link e calculamos seus respectivos valores de R-square. R-square consiste de uma medida estatística de quão próximos os dados estão da linha de regressão ajustada.

Hypothesis 1: D and t matter

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.data.1 = list(data.frame("H" = numeric(400), "C" = numeric(400), stringsAsFactors=FALSE), 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.data.1$H[[cc]] = HC.BP$H[elements]
    HC.data.1$C[[cc]] = HC.BP$C[elements]
  }
  b = b + 4000
}

HC = data.frame("H" = numeric(400), "C" = numeric(400), stringsAsFactors=FALSE)
lm.fit.logit.1 = lm.fit.probit.1 = lm.fit.cloglog.1 = lm.fit.loglog.1 = array(list(), 40)
r.squared.1 = data.frame("logit" = numeric(40), "probit" = numeric(40), "cloglog" = numeric(40),"loglog" = numeric(40), stringsAsFactors=FALSE)
for(i in 1:40){
  HC$H = HC.data.1$H[[i]]
  HC$C = HC.data.1$C[[i]]
  
  lm.fit.logit.1[[i]] = betareg(data = HC, formula = C ~ H, link = "logit")
  r.squared.1$logit[i] = lm.fit.logit.1[[i]]$pseudo.r.squared
  
  lm.fit.probit.1[[i]] = betareg(data = HC, formula = C ~ H, link = "probit")
  r.squared.1$probit[i] = lm.fit.probit.1[[i]]$pseudo.r.squared
  
  lm.fit.cloglog.1[[i]] = betareg(data = HC, formula = C ~ H, link = "cloglog")
  r.squared.1$cloglog[i] = lm.fit.cloglog.1[[i]]$pseudo.r.squared
  
  lm.fit.loglog.1[[i]] = betareg(data = HC, formula = C ~ H, link = "loglog")
  r.squared.1$loglog[i] = lm.fit.loglog.1[[i]]$pseudo.r.squared
} 

Hypothesis 2: D, t and N matter

HC = data.frame("H" = numeric(100), "C" = numeric(100), stringsAsFactors=FALSE)
lm.fit.logit.2 = lm.fit.probit.2 = lm.fit.cloglog.2 = lm.fit.loglog.2 = array(list(), 160)
r.squared.2 = data.frame("logit" = numeric(160), "probit" = numeric(160), "cloglog" = numeric(160),"loglog" = numeric(160), stringsAsFactors=FALSE)
for(i in 1:160){
  HC$H = HC.data.2$H[[i]]
  HC$C = HC.data.2$C[[i]]
  
  lm.fit.logit.2[[i]] = betareg(data = HC, formula = C ~ H, link = "logit")
  r.squared.2$logit[i] = lm.fit.logit.2[[i]]$pseudo.r.squared
  
  lm.fit.probit.2[[i]] = betareg(data = HC, formula = C ~ H, link = "probit")
  r.squared.2$probit[i] = lm.fit.probit.2[[i]]$pseudo.r.squared
  
  lm.fit.cloglog.2[[i]] = betareg(data = HC, formula = C ~ H, link = "cloglog")
  r.squared.2$cloglog[i] = lm.fit.cloglog.2[[i]]$pseudo.r.squared
  
  lm.fit.loglog.2[[i]] = betareg(data = HC, formula = C ~ H, link = "loglog")
  r.squared.2$loglog[i] = lm.fit.loglog.2[[i]]$pseudo.r.squared
} 
summary(r.squared.2)

     logit            probit          cloglog           loglog      
 Min.   :0.6430   Min.   :0.7119   Min.   :0.6399   Min.   :0.7644  
 1st Qu.:0.7304   1st Qu.:0.7858   1st Qu.:0.7281   1st Qu.:0.8270  
 Median :0.7619   Median :0.8124   Median :0.7598   Median :0.8501  
 Mean   :0.7600   Mean   :0.8107   Mean   :0.7581   Mean   :0.8469  
 3rd Qu.:0.7926   3rd Qu.:0.8359   3rd Qu.:0.7909   3rd Qu.:0.8673  
 Max.   :0.8527   Max.   :0.8859   Max.   :0.8511   Max.   :0.9091  

Como visto acima, os melhores ajustes são obtidos ao aplicar a função de link loglog sobre o modelo de regressão. No entanto, em geral um modelo ajusta bem os dados se as diferenças entre os valores observados e os valores previstos do modelo são pequenas e imparciais. Ao observar os gráficos de diagnótico para as regressões vemos que tal modelo ainda não é o adequado para os nossos dados.

Hypothesis 1: D and N matter

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

Hypothesis 2: D, t and N matter

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

Concluímos que o modelo ainda não obtem um bom ajuste após verificarmos os plots de suas respectivas regressões.

Hypothesis 1: D and N matter

HC = data.frame("H" = HC.data.1$H[[1]], "C" = HC.data.1$C[[1]], stringsAsFactors=FALSE)
ggplot(HC, aes(x = H, y = C)) +
  geom_point() +
  scale_fill_grey() +
  geom_line(aes(y = predict(lm.fit.loglog.1[[1]], HC))) 

Hypothesis 2: D, t and N matter

HC = data.frame("H" = HC.data.2$H[[1]], "C" = HC.data.2$C[[1]], stringsAsFactors=FALSE)
ggplot(HC, aes(x = H, y = C)) +
  geom_point() +
  scale_fill_grey() +
  geom_line(aes(y = predict(lm.fit.loglog.2[[1]], HC)))