Aplicação de Network Radial Base Function para modelagem de frequência cardíaca

Descrição do desafio

Construir um algoritmo em linguagem R, utilizando NRBF (Network Radial Base Funcion) para modelar o comportamento de um fenômeno descrito na base de dados no arquivo Anexo. Especificamente para a aplicação do NRBF, use:
1. Função radial Gaussiana como função de ativação
2. Técnica da matriz pseudo-inversa, para resolução do sistema de coeficientes das funções
3. Método de K-means para seleção dos polos
4. Número de polos igual ou maior à quantidade de variáveis (6)

Como resultado, esperamos receber um código em R contendo:
1. O algoritmo de modelagem
2. Rotinas para a visualização da aplicação do algoritmo
3. Análise de assertividade e erro
4. Aplicação do código para os dados do arquivo Anexo

Introdução

Regressão paramétrica é um tipo de modelagem cuja característica distintiva é a ausência (completa ou quase) de conhecimento a priori sobre a forma da função que está sendo estimada. Assim, o conjunto de formas que a função pode assumir é abrangente. Como consequência, haverá um número elevado de parâmetros estimados que não possuem uma interpretação física isolada.

As redes neurais para aprendizado supervisionado do tipo RBF (Radial Basis Function Neural Network) são um exemplo de modelo não-paramétrico baseado na teoria de aproximação de funções, que trata o problema genérico de aproximar uma função \(y(x)\) por uma função de aproximação \(f(w, x)\) dado um número fixo de parâmetros \(w\).

Uma função de ativação de base radial é caracterizada por apresentar uma resposta que decresce (ou cresce) monotomicamente com a distância em relação a um ponto central. O centro e a taxa de decrescimento em cada diração são alguns dos parâmetros a serem definidos. Uma função de base radial típica é a função gaussiana, dada na forma:

\(h_j(x) = \exp(\dfrac{-(x-c_j)^2}{{r_j^2}})\)

As funções de base radial podem ser utilizadas como funções-base em qualquer tipo de modelo de regressão não-linear e como função de ativação de qualquer tipo de rede multicamada. As redes RBF sempre apresentam uma única camada intermediária, os neurônios de saída são sempre lineares e os neurônios da camada intermediária têm uma função de base radial como função de ativação.

Um exemplo visual da rede RBF é mostrada na figura abaixo.

Uma das abordagens para o treinamento das redes neurais RBF é a adaptação supervisionada dos pesos da camada de saída, empregando técnicas como pseudo-inversão.

No caso de algoritmos que se ocupam apenas com o ajuste dos pesos da camada de saída de uma rede RBF, é necessário estabelecer algum critério para a fixação dos centros. Uma alternativa é auto-organizar os centros, de acordo com a distribuição ds dados de entrada, por meio do algoritmo k-means.

Em suma, existem três tipos de parâmetros que precisam ser aprendidos por uma rede neural RBF:
* Os centros da função de ativação;
* As dispersões (taxa de decrescimento) da função de ativação;
* Os pesos da camada intermediária para a camada de output.

Código

Para este exercício, são usados os seguintes pacotes:

library(readxl)
library(GGally)
library(tidyverse)

A leitura do arquivo e a mudança da ordem das colunas é feita abaixo:

arquivo <- "Base Projeto Daemon - Quant 2018.xlsx"
df <- readxl::read_excel(arquivo, sheet = "Sheet1",
                         col_names = TRUE, range = "B5:H128")
df <- as.data.frame(df)
# mudar ordem das colunas
df <- df[, rev(1:ncol(df))]
knitr::kable(head(df))

nuVHF	nuHF	nuLF	a1	RMSSD	SDNN	fc
0.0120016	0.2751218	0.2343089	0.2645995	0.0335986	0.0515247	1.393146
0.0347053	0.2907528	0.2066793	0.2863589	0.0263905	0.0358156	1.444002
0.0360863	0.2920406	0.1471003	0.2780337	0.0234379	0.0301764	1.457103
0.0184709	0.1712113	0.1472847	0.1994207	0.0245995	0.0459377	1.462573
0.0212941	0.1456406	0.1891565	0.1937349	0.0232357	0.0437853	1.467897
0.0168037	0.2222374	0.1856247	0.2385001	0.0278113	0.0452344	1.434446

O objetivo deste trabalho é modelar a variável de frequência cardíaca em função de todas as outras variáveis do conjunto de dados.

Modelagem

No código abaixo, é criada uma função que realiza as seguintes tarefas:
* Extrai uma quantidade k de centros da matriz de input X. Por padrão, k corresponde à quantidade de colunas na matriz de input X;
* Aplica uma função de ativação gaussiana para obter a matriz \(\phi\), onde cada elemento da matriz N x k corresponde a \(\gamma \times h(\|x_i - t_i\|)\), onde \(\gamma\) corresponde à dispersão da função de ativação;
* Obtem pela técnica da matriz pseudo inversa o vetor de pesos, a partir da matriz \(\phi\) e do vetor de output (ou target) y;
* Retorna valores ajustados para o vetor de output y a partir de uma matriz de input X_test.

nrbf <- function(X, Y, X_test = X, k = ncol(X), gamma = 1.0, seed = 123, plot = TRUE){
  library(corpcor)
  library(neuralnet)
  set.seed(seed)
  #### Definição dos argumentos:
  # X: Matriz de input 
  # Y: Matriz de output
  # k: número de centros (polos)
  # gamma: parametro de aprendizado
  
  N <- nrow(X) # numero de observacoes
  # criar funcao de ativacao gaussiana
  ativ_gaussiana <- function(gamma, x, y){
    # x e y sao vetores numericos
    
    v <- as.matrix(x - y)
    modulo <- norm(v, "F")
    exp(-gamma * modulo^2)
  }
  
    
  # criar clusteres em matriz de input normalizada
  cl <- kmeans(scale(X), centers = k)
  cl_centros <- cl$centers # (Matriz de dimensões k x ncol(X))
  # criar matriz vazia Phi para preencher  posteriormente
  matriz_phi <- matrix(NA_real_, nrow = N, ncol = k + 1)
  
  #browser()
  # iterar em cada linha
  for (l in 1:N){
    # Preencher a primeira coluna com 1 (coluna de bias)
    matriz_phi[l, 1] <- 1
    # iterar em cada coluna 
    for (c in 1:k){
      # calcular modulo do vetor das diferencas
      # preencher celula (com excecao da primeira coluna, que eh a de bias)
      matriz_phi[l, c + 1] <- ativ_gaussiana(gamma = gamma, x = X[l, ], y = cl_centros[c, ])
    }
  }
  
  w <- corpcor::pseudoinverse(t(matriz_phi) %*% matriz_phi) %*% t(matriz_phi) %*% Y 
  
  # obter previsoes para Y
  # inicializar vetor de previsoes a partir do bias do vetor de pesos
  prev <- rep(w[1], N)
  for (n in 1:N){
    for (j in 1:nrow(cl_centros)){
      prev[n] <- prev[n] + w[j + 1] * ativ_gaussiana(gamma = gamma, x = X_test[n, ], y = cl_centros[j, ])
    }
  }
  
  # plotar rede neural
  if (plot){
    # combinar input e target em um dataframe
    df_nnet <- data.frame(input, target = target)
    model_formula <- paste0("target ~ ", paste0(colnames(input), collapse = " + "))
    nn <- neuralnet(model_formula, data = df_nnet, hidden = k)
    plot(nn, information = FALSE, show.weights = FALSE)
  }
  
  list(weights=w, centers=cl_centros, gamma=gamma, prev = prev)  
}

Após criar a função, o código abaixo mostra como é simples ajustar o modelo à matriz de input deste exercício:

input <- as.matrix(df[, -7]) # a ultima coluna corresponde ao vetor y que se deseja prever
target <- as.matrix(df[, 7])
modelo <- nrbf(input, target)

# visualizar resultados do modelo
modelo$weights

                 [,1]
[1,]     0.4299523943
[2,] -7643.2332381372
[3,]     4.9617618704
[4,]    -1.8867322271
[5,]    69.6572399773
[6,]    12.6253134136
[7,]     0.1175875511

modelo$centers

          nuVHF          nuHF           nuLF            a1         RMSSD          SDNN
1  0.5837561752  0.6350442465 -0.66467743516  1.4381165768  2.2253072683  2.1569218102
2  0.4543580595  1.4390765402 -0.72753512269  0.4754167243  0.4489823963  0.2195427664
3 -0.6797732920 -0.8771208817 -0.36038859061 -1.0627150483 -0.6210457527 -0.2637736570
4  1.9941167775 -0.2407254829  0.87993133176 -0.1132062724 -0.3152785222 -0.5166736536
5 -0.5319770927 -0.8128502961  1.74751201173 -0.9378453344 -0.5700864149 -0.3475751648
6 -0.5577582063 -0.2105481031 -0.05084140639  0.2103505167 -0.4882146668 -0.5692230101

Análise dos resultados do modelo

O código abaixo mostra a qualidade do ajuste do modelo:

prev <- modelo$prev
plot(target)
lines(prev, col = "blue", lty = 2)

O gráfico acima mostra que o ajuste obtido é relativamente bom, mas falhou em ajustar os pontos mais à esquerda no gráfico. A função abaixo calcula o MAPE (Mean Absolute Percentage Error) do modelo de função de base radial.

# obter acuracia
mape <- function(real, previsto){
  erro <- real - previsto
  mean(abs(erro)/real)
}
mape(target, prev)

[1] 0.07508241996

Portanto, o modelo apresenta um erro de 7,5%.

Comparação com o modelo de regressão linear múltipla (MRLM)

O modelo de regressão linear múltipla abaixo é criado para comparar sua acurácia com o modelo de função de base radial:

# regressao nultipla linear
modelo_lin <- lm(fc ~ ., data = df)
prev_lin <- fitted(modelo_lin)
mape_lin <- mape(target, prev_lin)
mape_lin

[1] 0.05113152375

O erro do modelo de regressão múltipla é de 5%, o que significa que sua qualidade de ajuste é superior ao modelo de função de base radial. O gráfico abaixo mostra essa comparação visualmente:

plot(target)
lines(prev, col = "blue", lty = 2)

lines(fitted(modelo_lin), col = "red", lty = 2)
legend("topright", lty = 1, col = c("blue", "red"),
       legend = c("NRBF", "MRLM"))

Otimização do número de polos (centros)

É possível com certa facilidade obter a qualidade de ajuste do modelo NRBF para diferentes valores do número de polos (k) usados no modelo para descobrir qual número maximiza a acurácia do modelo.

O código abaixo obtem o MAPE do modelo NRBF para valores de K de 1 a 60. Os resultados são mostrados no gráfico abaixo, no qual a linha vertical corresponde ao MAPE do modelo de regressão linear múltipla.

# simular valores de k para observar se o rmse cai
vetor_k <- 1:60
# criar dataframe para armazenar resultados do rmse
df_mape <- data.frame(k = vetor_k, mape = rep(NA, length(vetor_k)))
for (i in 1:length(vetor_k)){
  prev_loop <- nrbf(input, target, k = vetor_k[i], plot = FALSE)$prev
  df_mape[i, 2] <- mape(target, prev_loop)
}
# analisar resultados visualmente
ggplot(df_mape, aes(x = k, y = 100 * mape)) +
  geom_point() +
  scale_x_continuous(breaks = seq(0, 60, 5)) +
  scale_y_continuous(breaks = scales::pretty_breaks(n = 10)) +
  geom_hline(yintercept = 100 * mape_lin, linetype = "dashed") + 
  labs(x = "Quantidade de polos", y = "MAPE (%)",
       title = "Relação entre quantidade de polos e o erro do modelo \n de função de base radial")

Uma informação muito útil extraída do gráfico é que a partir de \(k = 20\) o modelo NRBF se torna melhor que o de regressão linear múltipla. O erro do modelo é reduzido sensivelmente de \(k = 30\) para \(k = 31\) e se estabiliza desde então, o que nos leva a acreditar que este seja o número ideal de polos para o modelo.

O gráfico abaixo mostra a qualidade do ajuste do modelo NRBF com 31 polos no conjunto de dados do exercício.

# obter modelo com 31 clusteres
melhor_modelo <- nrbf(input, target, k = 31, plot = TRUE)

plot(target)
lines(melhor_modelo$prev, col = "blue", lty = 2)

Referências

Notas de aula dos professores Fernando J. Von Zuben & Levy Boccato

Explicação teórica de NRBF com aplicação de código R