Report 5 - Duda kernel function

if(!require(drc)) install.packages("drc")
if(!require(geoR)) install.packages("geoR")
if(!require(stats)) install.packages("stats")
if(!require(ggplot2)) install.packages("ggplot2")
if(!require(ggrepel)) install.packages("ggrepel")
if(!require(nortest)) install.packages("nortest")
if(!require(graphics)) install.packages("graphics")
if(!require(EnvStats)) install.packages("EnvStats")
if(!require(rcompanion)) install.packages("rcompanion")
if(!require(aomisc)) install_github("OnofriAndreaPG/aomisc")

Neste documento, utilizamos um modelo de regressão usando uma função base não-linear. Para isto, primeiramente aplicamos a transformação de Box-Cox na variável resposta do nosso problema, ou seja, nos valores de complexidade estatística:

HC.BP = data.frame("H" = numeric(20800), 
                   "C" = numeric(20800),
                   "D" = numeric(20800),
                   "t" = numeric(20800), 
                   "N" = numeric(20800),
                   "index" = numeric(20800),
                   stringsAsFactors=FALSE)
HC.BP$index = as.factor(c(1:20800))
HC.BP$N = as.factor(rep(c(rep(1e+04, 400), 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,1300), rep(2,1300), rep(3,1300), rep(4,1300)), 4))
file.csv = data.frame(read.csv("../Data/HC_series_fk0.csv"))
HC.BP$H = file.csv[,1]
HC.BP$C = file.csv[,2]
HC.BP$D= as.factor(file.csv[,3])
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
HC.BP$C = T.norm
plot(x = HC.BP$H, y = HC.BP$C)

Para selecionar o kernel que melhor ajuste o nosso modelo adotamos uma abordagem empírica. Ou seja, analisamos o problema, plotamos os dados e observamos que eles seguem um certo padrão. Um problema com a regressão não linear é que ele funciona iterativamente: precisamos fornecer suposições iniciais para os parâmetros do modelo e o algoritmo os ajusta passo a passo, até que (esperemos) converja na solução aproximada de mínimos quadrados. Portanto, usamos as funções do R, incluindo as rotinas de autoinicialização adequadas, o que simplificou bastante o processo de ajuste. Usamos como suporte o pacote ‘aomisc’.

Sabemos por testes anteriores que tais dados obtiveram resultados positivos utilizando a função de link loglog utilizada na regressão beta. Desse modo, baseando-se em tal usamos a seguinte parametrização:

\[ C' = a + b \times log(-log(1-H)) \]

Uma vez que \(\lim_{H \to 1}a + b \times log(-log(1-H)) = - \infty\), a seguinte modificação foi feita:

\[ C_i' = a + b \times (1_{(H_i < 1)}log(-log(1-H)) + 1_{(H_i = 1)}) \]

"LinkFunctionDuda" <- function(fixed = c(NA, NA), names = c("a", "b"))
{
  ## Checking arguments
  numParm <- 2
  if (!is.character(names) | !(length(names) == numParm)) {stop("Not correct 'names' argument")}
  if (!(length(fixed) == numParm)) {stop("Not correct 'fixed' argument")}
  
  ## Fixing parameters (using argument 'fixed')
  notFixed <- is.na(fixed)
  parmVec <- rep(0, numParm)
  parmVec[!notFixed] <- fixed[!notFixed]
  
  ## Defining the non-linear function
  fct <- function(x, parm)
  {
    
    
    parmMat <- matrix(parmVec, nrow(parm), numParm, byrow = TRUE)
    parmMat[, notFixed] <- parm
    
    a <- parmMat[, 1]; b <- parmMat[, 2]
    y = rep(0, length(x))
    for(i in 1:length(x)){
      if(x[i] == 1)
       y[i] = a[i] + b[i]
      else
       y[i] = a[i] + b[i]*log(-log(1-x[i]))
    }
    y
  }
  
  ## Defining self starter function
  ssfct <- function(dataf)
  {
    x <- dataf[, 1]
    y <- dataf[, 2]
    
    #regression on pseudo y values
    pseudoY <- (y)
    pseudoX <- (x)
    coefs <- coef( lm(pseudoY ~ pseudoX) )
    a <- (coefs[1])
    b <- (coefs[2])
    return(c(a, b)[notFixed])
  }
  
  ## Defining names
  pnames <- names[notFixed]
  
  ## Defining descriptive text
  text <- "Power curve (Freundlich equation)"
  
  ## Returning the function with self starter and names
  returnList <- list(fct = fct, ssfct = ssfct, names = pnames, text = text, noParm = sum(is.na(fixed)))
  
  class(returnList) <- "drcMean"
  invisible(returnList)
}

Calculando a regressão como uma função apenas de H

dimension = c(3,4,5,6)
lm.fit = array(list(), 40)
N = c(1e+04, 2e+04, 3e+04, 4e+04, 5e+04, 6e+04, 7e+04, 8e+04, 9e+04, 1e+05)
HC.test.1 = data.frame("H" = numeric(1600), 
                      "C" = numeric(1600),
                      "D" = numeric(1600),
                      "N" = numeric(1600),
                      "index" = numeric(1600),
                      stringsAsFactors=FALSE)
HC.test = data.frame("H" = numeric(400), 
                     "C" = numeric(400),
                     "D" = numeric(400),
                     "N" = numeric(400),
                     "index" = numeric(400),
                     stringsAsFactors=FALSE)
b = cc = 0
for(i in 1:length(dimension)){
  for(j in 1:length(N)){
      
      if(i*j == 1){
        cc = cc + 1
        a = c((((j - 1) * 400) + 1):(j * 400))
        elements = c(a + b, a + b + 1300, a + b + 2600, a + b + 3900)
        
        HC.test.1$H = HC.BP$H[elements]
        HC.test.1$C = HC.BP$C[elements]
        HC.test.1$D = HC.BP$D[elements]
        HC.test.1$N = HC.BP$N[elements]
        HC.test.1$index = HC.BP$index[elements]
        
        lm.fit[[cc]] = drm(C ~ H, data = HC.test.1, fct = LinkFunctionDuda())
        
      }
      else{
        cc = cc + 1
        a = c((((j - 1) * 100) + 1):(j * 100))
        elements = c(a + b, a + b + 1300, a + b + 2600, a + b + 3900)
      
        HC.test$H = HC.BP$H[elements]
        HC.test$C = HC.BP$C[elements]
        HC.test$D = HC.BP$D[elements]
        HC.test$N = HC.BP$N[elements]
        HC.test$index = HC.BP$index[elements]
        
        lm.fit[[cc]] = drm(C ~ H, data = HC.test, fct = LinkFunctionDuda())
      }
  } 
  b = b + 5200
}
a = c(1:400)
elements = c(a, a + 1300, a + 2600, a + 3900)
HC.test.1$H = HC.BP$H[elements]
HC.test.1$C = HC.BP$C[elements]
HC.test.1$D = HC.BP$D[elements]
HC.test.1$N = HC.BP$N[elements]
HC.test.1$index = HC.BP$index[elements]
model = drm(C ~ H, data = HC.test.1, fct = LinkFunctionDuda())
      
ggplot(HC.test.1, aes(x = H, y = C)) +
  geom_point() +
  scale_fill_grey() +
  geom_line(aes(y = predict(model, HC.test.1, type = 'response'))) 

No entanto, como podemos verificar abaixo, tal problemática ainda não é solucionada pois para valores muito próximos de \(1\) os resíduos ainda continuam aumentando:

plot(fitted(model), residuals(model, typeRes = "standard"))

pr = predict(model, HC.test.1, type = 'response')
ggplot(HC.test.1, aes(x = H, y = C - pr)) +
  geom_point() +
  scale_fill_grey() +
  geom_line(aes(y =0)) 

\(Q-Qplot\) dos resíduos

qqPlot(residuals(model, typeRes = "standard")) 

[1] 1591 1126

Histograma dos resíduos

plotNormalHistogram(residuals(model, typeRes = "standard")) 

hist(residuals(model, typeRes = "standard"), breaks = 200)

Logo, concluímos que com o kernel não-linear proposto conseguimos obter resíduos menos variantes e um bom ajuste em relação aos dados transformados. Entretanto, tal função ainda apresenta problemas quando \(H \to 1\).