Início > dados, delineamento, matemática, testes > PCA de três dimensões

PCA de três dimensões

Olá a todos. O objetivo desse post é mostrar como fazer a análise de componentes principais em três dimensões. Mas esse post tem o destaque não ter sido escrito sozinho… Sem a participação da Msc Thalita (ou “Hermínia”, ninguém vai entender!) não teria post. Agradeçam a ela, mas reclamem comigo! Pois bem, Vamos usar como exemplo dados de um estudo com fatorial duplo (4 \times 4), em que foram tomadas quatro características (dados disponíveis para teste!).

Vou aproveitar e deixar aqui também um alô para o Dsc Guilherme pela ajuda nos detalhes do script, Valeu!

Um pouco de teoria:

Experimentos fatoriais (n \times n) caracterizam-se pela apresentação dos dados em tabelas de duas entradas (matriz), sendo que cada casela da tabela contém a resposta média de cada combinação dos fatores sob análise. Todas possibilidades desses experimentos estão fora do escopo desse post, por isso vamos um pouco além…

Acrescentando-se um outro fator, ou ainda outra característica (n \times n \times n, por exemplo), com 3 interações duplas e uma tripla, isso fica um pouco mais complexo. Neste caso, os dados são organizados em “CUBO” (arranjo de três entradas).

Nesta situação não é possível aplicar a análise de componentes principais usual. Uma alternativa é a utilização da decomposição PCAn, proposta por Tuckey em 1964.

Portanto nesse post pretendemos (afinal esse foi um trabalho em equipe) mostrar como realizar a decomposição de um arranjo de três entradas (CUBO). E mais que isso, mostrar uma função que realiza a escolha entre os vários modelos possíveis que variam conforme a dimensão dos dados.

Vamos ao código… mas calma, primeiro preparamos os dados!

Veja que usamos um pacote (PTAk, disponível no CRAN). Caso você não o tenha, proceda com install.packages(‘PTAk’)!

require(PTAk) # pacote analises pscicrometricas

## le e transforma os dados em fatores
dados <- transform(read.table("http://dl.dropbox.com/u/38195533/dados.txt")
                              na.string = '.',
                              header = TRUE,
                              dec = ',',
                              sep = '\t'),
                   elemento = factor(elemento),
                   residuo = factor(residuo),
                   metodo = factor(metodo))

CUBO.inicial <- with(dados,
                     tapply(conc, list(residuo, elemento, metodo),
                            mean, na.rm = TRUE))
## padronização dos dados
CUBO.standard <- apply(CUBO.inicial, 2, scale, center = FALSE)

## reconstrucao da planilha
estrutura <- expand.grid(elemento = factor(1:4),
                         residuo = factor(1:4),
                         metodo = factor(1:4))

plan <- data.frame(estrutura, conc = c(CUBO.standard)) ## dados

ajuste <- aov(conc ~ elemento * residuo * metodo - elemento:residuo:metodo,
              data = plan)
erros <- ajuste$residuals

## CUBO
Z <- tapply(erros, list(plan$residuo, plan$elemento, plan$metodo), mean) 

Nesse trecho são lidos os dados, na forma usual de uma planilha de dados (uma coluna para cada variável). Pela natureza dos dados, estes são padronizados (N(\mu = 0, \sigma^2 = 1)). Depois a panilha de dados é “reconstruída” e feito o ajuste do modelo linear respectivo (sem a interação tripla — objeto do nosso estudo!). Do resultado da análise extraímos os resíduos e construímos o CUBO.

Quase tudo pronto! O próximo trecho é a função que efetivamente faz os vários ajustes possíveis (isso é função da dimensão dos dados).

otimim.PCAn <- function(cubo) {
  dimensoes <- dim(Z) - 1
  diagnostico <- vector('list', length = prod(dimensoes - 1))
  ajustes <- vector('list', length = prod(dimensoes - 1))
  contador <- 0
  for(i in 2:dimensoes[ 1 ])
    {
      for(j in 2:dimensoes[ 2 ])
        {
          for(k in 2:dimensoes[ 3 ])
            {
              modelo <- c(i, j, k)
              tamanho <- sum(modelo)
              ## decomposicao Tucker -- PTAk
              tucker <- PCAn(Z, # array
                             dim = modelo, # dimensoes
                             test = 1E-12,
                             Maxiter = 1000,
                             smoothing = FALSE,
                             verbose = FALSE,
                             file = NULL)
              ## diagnostico
              porcentagem <- tucker[[ 3 ]]$pct # % explicada
              contador <- contador + 1
              ajustes[[ contador ]] <- tucker
              diagnostico[[ contador ]] <- c(modelo, tamanho, porcentagem)
            }
        }
    }
  ## diagnosticos dos ajustes
  tab.ajustes <- data.frame(do.call(rbind, diagnostico))
  tab.ajustes <- tab.ajustes[order(-tab.ajustes[, 4], -tab.ajustes[, 5]), ]
  por.total <- split(tab.ajustes, tab.ajustes[, 4])
  melhores <- vector('list', length = length(por.total))
  for(i in 1:length(melhores))
    {
      melhores[[i]] <- por.total[[i]][1, ]
    }
  ## melhores ajustes por dimensao
  tab.melhores <- do.call(rbind, melhores)
  tab.melhores <- tab.melhores[order(-tab.melhores[, 4]), ]
  ## melhor modelo
  melhor <- tab.melhores[(which((tab.melhores[,5]-
                          c(tab.melhores[,5][-1],0))>=5)[1]+1),]
  ## resposta -- melhor ajuste
  return(list(escolhido = melhor,
              resposta = ajustes[[ as.numeric(rownames(melhor)) ]]))
}

UFA!

Feita a função agora é só usa-lá (simples não?). Um uso possível dessa decomposição é representar e interpretar os dados graficamente em um biplot. Nessa condição de arranjos de três entradas são possíveis diferentes formas de gráficos biplots. Mas o objetivo do nosso post não é apresentá-los.

O trecho final aplica a função recém-construída (e retorna um controle de tempo de cada ajuste). Mostramos também como extrair as matrizes (A, B e C) e o arranjo G do melhor ajuste.

As matrizes de resposta e o arranjo G é que são úteis na construção dos gráficos. Aproveitem!

ajuste.PCA <- otimim.PCAn(Z)

A <- ajuste.PCA$resposta[[1]]$v
B <- ajuste.PCA$resposta[[2]]$v
C <- ajuste.PCA$resposta[[3]]$v

G <- ajuste.PCA$resposta[[3]]$coremat

summary(ajuste) # ANOVA do modelo
summary.PTAk(ajuste.PCA$resposta) # quanto cada componente explica

Voi lá…, Até mais pessoal!

About these ads
  1. Nenhum comentário ainda.
  1. No trackbacks yet.

Deixe um comentário

Preencha os seus dados abaixo ou clique em um ícone para log in:

Logotipo do WordPress.com

Você está comentando utilizando sua conta WordPress.com. Sair / Alterar )

Imagem do Twitter

Você está comentando utilizando sua conta Twitter. Sair / Alterar )

Foto do Facebook

Você está comentando utilizando sua conta Facebook. Sair / Alterar )

Foto do Google+

Você está comentando utilizando sua conta Google+. Sair / Alterar )

Conectando a %s

Seguir

Obtenha todo post novo entregue na sua caixa de entrada.

Junte-se a 52 outros seguidores

%d blogueiros gostam disto: