Fuzzy c-means, um exemplo!

Pois bem, por sugestão do Professor Julio da Motta Singer apresentarei um script para realizar uma análise de agrupamento no escopo da lógica fuzzy. Seja um experimento fatorial (a mesma natureza dos dados aqui disponibilizados) a seguinte implementação realiza a clusterização de uma dado fator marginalmente ao outro fator.

A heurística da lógica fuzzy tem um incremento interpretativo de, ao contrário do agrupamento k-means, que é dicotômico, associar a cada cluster uma “probabilidade” de associação. Todavia, no método arbitramos uma constante () de fuzzyficação. Assim o seguinte script realiza o agrupamento fuzzy (por intermédio da função cmeans() da biblioteca e1071) para vários valores assim é possível observar o comportamento (decaimento) do erro associado a clusterização e também a disposição dos fatores em relação ao cluster.

Infelizmente (ou para alegria geral) não vou me atrever a discorrer da teoria envolvida nessa análise e um estudo específico deixo a cargo dos maiores interessados.

Todo esse mise en cene, me soa um pouco genérico demais, eu sei (–Me desculpem!). Mas acredito que o script possa ser útil de alguma forma.

Aproveito o ensejo para agradecer as sugestões do Prof Julio da Motta Singer ao trabalho e o incentivo a esse post. Também deixo um alô ao colegas que participaram, de uma forma ou de outra!

Vamos ao trabalho, primeiro instalando (para quem ainda não o fez) e carregando a biblioteca e o conjunto de dados!

## install.packages('e1071')
library(e1071) # carrega pacote com implementação 'c-means'

dados <- transform(read.table("http://dl.dropbox.com/u/38195533/fac_exemplo.txt", # arquivo
                              header = TRUE,
                              sep = '\t',
                              na.string = '.',
                              dec = ','),
                   fatorA = factor(fatorA),
                   fatorB = factor(fatorB))

str(dados)

ajuste <- lm(resp ~ fatorA + fatorB, data = dados)
anova(ajuste) # ANOVA

matriz <- tapply(dados$resp, list(dados$fatorB, dados$fatorA), mean)

A primeira etapa foi realizada. Não posso passar ou próximo trecho sem deixa-los a par do quem está por vir! O agrupamento que “tentamos” é com relação ao um padrão de comportamento do fator sob análise, em relação ao outro fator, assim constituímos ideótipos e assim procedemos. O modo como constituimos esses ideótipos é aplicado unicamente a natureza desse conjunto de dados, mas outros ideótipos, ao seu critério podem ser usados, fique a vontade!

Antes de passarmos agora a análise em si manipulamos os dados e constituímos os centróides:

FinWil <- apply(matriz, 2, mean) - mean(matriz) # Finlay & Wilkinson

altos <- which(FinWil >= 0) # ambientes acima da media
baixo <- which(FinWil < 0) # ambientes abaixo da media

medias <- cbind(apply(matriz[, altos], 1, mean),
                apply(matriz[, baixo], 1, mean))

##--------------------------------------------------------------------
## IDEOTIPOS
## I   - maximo em todos 
## II  - maximo nos altos minimo nos baixos
## III - minimo nos altos maximo nos baixos
## IV  - minimo em todos 
## ADICIONAIS:
## V   - media em todos
## VI  - maximo nos altos e media nos baixos
## VII  - media nos altos maximo nos baixos

I <- c(max(matriz[, altos]), max(matriz[, baixo]))
II <- c(max(matriz[, altos]), min(matriz[, baixo]))
III <- c(min(matriz[, altos]), max(matriz[, baixo]))
IV <- c(min(matriz[, altos]), min(matriz[, baixo]))
V <- c(mean(matriz[, altos]), mean(matriz[, baixo]))
VI <- c(max(matriz[, altos]), mean(matriz[, baixo]))
VII <- c(mean(matriz[, altos]), max(matriz[, baixo]))

centroides <- rbind(I, II, III, IV, V, VI, VII) # centróides

Resta agora realizar a análise em si, veja que o laço for() repete a análise para vários valores da constante e de cada uma, armazena as “probabilidades” de associação a cada centróide e o erro do ajustamento!

ms <- seq(1.2, 10, by = .1)
lista <- array(NA, dim = c(nlevels(dados$fatorB), 
                           nrow(centroides), 
                           length(ms)))
erros <- matrix(NA, nrow = length(ms), ncol = 1)
dimnames(lista) <- list(levels(dados$TRAT), 
                        c('I', 'II', 'III', 'IV', 'V', 'VI', 'VII'), 
                        ms)

for(i in 1:length(ms))
  {
    j <- ms[i]
    ## função de cluster 'fuzzy c-means'
    grupamentos <- cmeans(medias, # matriz dos dados
                          centroides, # centróides dos clusters
                          verbose = FALSE, # não imprime histórico
                          dist = 'euclidean', # distância empregada
                          m = j, # constante de desfuzzyficação
                          method = 'cmeans') # método
    lista[,,i] <- grupamentos$membership
    erros[i,] <- grupamentos$withinerror
  }

Pronto… feito, quer dizer, quase! Realizamos portanto a análise para os vários valores de resta agora “escolher” um deles e partir para o resto! Até a próxima…

Equação junto ao gráfico na lattice

Reposta à salinidade para quatro espécies de craca. A regulação osmótica está ligada à capacidade adaptativa do organismo ao ambiente. O gráfico mostra os valores observados e a curva de regressão estimada acompanhada de bandas de confiança (95%). As equações em cada gráfico representam o modelo ajustado que descreve a curva.

Os polinômios são os modelos de regressão mais empregados em algumas áreas do conhecimento. Isso não quer dizer que eles sejam modelos bons. Na verdade, isso depende do cada um define como um bom modelo. Os polinômios são uma aproximação local. Isso significa que eles não devem ser usados para predições fora da amplitude observada das covariáveis (região experimental). Em outras palavras, não se deve fazer extrapolações, mas pode-se, com certos cuidados, fazer interpolações. Os polinômios são desprovidos de interpretação a partir do segundo grau. A partir do segundo grau, apenas o intercepto tem interpretação pois é o valor da função quando a covariável (x) é zero. Por isso, interpretações em geral são perdas de tempo. O que pode-se fazer são julgamentos quando ao sinal do coeficiente de maior grau que por exemplo, num modelo quadrático, indica se a curvatura é para cima (ou não). Outro ponto é que os polinômios apresentam problema de colinearidade. Para superar essa deficiência pode-se centrar a covariável ou usar polinômios ortogonais. Apesar do polinômio já ser sem interpretação, o uso de polinômios ortogonais é menos interpretável ainda pois os coeficientes não são taxas em relação à x, x², x³…, mas para uma transformação linear de x, essa que induz a ortogonalidade. Resumindo, os polinômios ortogonais são úteis pois minimizam o problema de colineridade. Então, do ponto de vista de reportar uma equação de predição, aquela que descreve o valor esperado de y (resposta) em função de x, é mais útil usar o polinômio cru. Isso porque, mesmo sem ter interpretação, toda pessoa sabe traçar a curva de um polinômio, é algo familiar, visto desde o ensino médio, e novamente, não que isso seja um ponto positivo para um modelo. De qualquer forma, o post de hoje ensina a ajustar polinômios e representar no gráfico os valores observados, os valores ajustados com intervalo de confiança e a equação dentro do gráfico. Os procedimentos abaixo eu desenvolvi para atender demandas de consultoria. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# definições da sessão

require(lattice)
require(latticeExtra)
require(grid)

#-----------------------------------------------------------------------------
# leitura

osmo <- read.table("http://www.leg.ufpr.br/~walmes/ridiculas/osmo.txt",
                   header=TRUE, sep="\t")
str(osmo)

xyplot(Osm~Salinidade|Espécie, data=osmo, type=c("p","smooth"))

#-----------------------------------------------------------------------------
# vou desconsiderar o fato dos dados serem heterocedasticos e ajustar um lm
# em outra matéria trataremos os desvios de pressupostos
# ajuste usando polinômios ortogonais

m0 <- aov(Osm~Espécie*poly(Salinidade,2), data=osmo)
par(mfrow=c(2,2)); plot(m0); layout(1)

qqmath(~residuals(m0)|osmo$Espécie)
xyplot(sqrt(abs(residuals(m0)))~fitted(m0)|osmo$Espécie, type=c("p","smooth"))
# desvio de pressupostos, trataremos em outra matéria do blog

summary(m0) # inferências sem validade pois existe desvio de pressupostos

#-----------------------------------------------------------------------------
# partição das somas de quadrados, declara Salinidade dentro de Espécie

m1 <- aov(Osm~Espécie/poly(Salinidade,2), data=osmo)

names(coef(m1))
attr(terms(m1), "term.labels")
ncoef <- 2 # número de efeitos por expécie=2 (linear e quadrático)
ngrup <- 4 # número de níveis de espécie=4
split.list <- split(matrix(1:(ncoef*ngrup), ngrup, ncoef), f=1:ngrup)
names(split.list) <- levels(osmo$Espécie)
split.list <- lapply(split.list,
                     function(x){ names(x) <- c("lin","qua"); x})
split.list <- list(a=split.list)
names(split.list) <- attr(terms(m1), "term.labels")[2]

summary(m1, split=split.list)

split.list <- as.list(1:(ncoef*ngrup))
names(split.list) <- c(outer(levels(osmo$Espécie), c("lin","qua"),
                             paste, sep=": "))
split.list <- list(a=split.list)
names(split.list) <- attr(terms(m1), "term.labels")[2]

summary(m1, split=split.list)

#-----------------------------------------------------------------------------
# prepara as equações para colocar no gráfico, função para polinômiais

coefs <- coef(lm(Osm~-1+Espécie/(Salinidade+I(Salinidade^2)), data=osmo))
coefs <- as.numeric(formatC(coefs, format="f", digits=3, flag="-"))

get.expression <- function(coef, digits){
  coef <- as.numeric(formatC(coef, format="f", digits=digits, flag="-"))
  postext <- c("","*x", paste("*x^", 2:(length(coef)-1), sep=""))
  sig <- c(ifelse(coef[1]<0, "-", ""), ifelse(coef[-1]<0, "-", "+"))
  txt <- paste(paste(sig, abs(coef), postext, sep=""), collapse=" ")
  eval(parse(text=paste("expression(y==", txt, ")", sep="")))
}

coefs <- split(matrix(coefs, ncol=3), f=1:4)
eqn <- lapply(coefs, get.expression, digits=3)
eqn

#-----------------------------------------------------------------------------
# nomes científicos em itálico nas tarjas

fl <- as.expression(lapply(levels(osmo$Espécie),
    function(x){ bquote(italic(.(x))) }))
fl

#-----------------------------------------------------------------------------
# predição

pred <- expand.grid(Espécie=levels(osmo$Espécie), Salinidade=seq(0,40,by=2))
pred <- cbind(pred, predict(m1, newdata=pred, interval="confidence"))

#-----------------------------------------------------------------------------
# carrega funções que serão necessárias para fazer as bandas de confiança

source("http://www.leg.ufpr.br/~walmes/ridiculas/ridiculas_functions.R")
grep("ciH", ls(), value=TRUE) # nome das funções (painéis) que vamos usar

#-----------------------------------------------------------------------------
# gráfico final

#png("f035.png", width=600, height=500)
xyplot(Osm~Salinidade|Espécie, data=osmo, col=1,
       xlab=expression(Salinidade~(ppm)),
       ylab=expression(Osmolalidade~(mOsm~kg^{-1}~de~H[2]*O)),
       strip=strip.custom(bg="gray90", factor.levels=fl))+
  as.layer(xyplot(fit~Salinidade|Espécie, data=pred, type="l", col=1,
                  ly=pred$lwr, uy=pred$upr,
                  prepanel=prepanel.ciH, panel=panel.ciH))
dimLay <- trellis.currentLayout()
for(col in 1:ncol(dimLay)){
  for(row in 1:nrow(dimLay)){
    trellis.focus("panel", col, row, highlight=FALSE)
    grid.text(label=eqn[[dimLay[row,col]]],
              x=0.03, y=1-0.03, just=c("left","top"))
  }
}
trellis.unfocus()
#dev.off()

#-----------------------------------------------------------------------------

Gráfico tridimensional com contornos de nível

Peso de 100 grãos em função do nível de saturação de água no solo e conteúdo de potássio fornecido pela adubação. Linhas no plano inferior representam a projeção das linhas de mesmo valor para peso de 100 grãos.

Nesse post vou apresentar uma função para adicionar contornos de nível em gráficos tridimensionais gerados pela lattice::wireframe(). Para ilustrar o uso da função eu fiz um ajuste de modelo de regressão múltipla aos dados do experimento com soja. Esses dados eu já usei em várias matérias aqui no Ridículas. Uma vez ajustado o modelo, os valores preditos são então representados pela wireframe() usando o panel.3d.contour que foi uma função que aprimorei das leituras que fiz das mensagens da r-help.

Um gráfico tridimensional como este impressiona quem o vê, porém esteja ciente das más impressões que um gráfico tridimensional pode dar pois ele não é de fato tridimensional e sim uma projeção tridimensional em um plano, o que acarreta algumas deformações. Todavia, achei por bem documentar e expor a função mas isso não significa que esta é a melhor representação. A função permite colocar os contornos no plano abaixo, acima e sobre a superfície (bottom, top, on). Tenha cuidado ao escolher o ângulo de observação do gráfico (screen=) e controle a amplitude do eixo z (zlim=) para que a superfície não esconda a projeção abaixo. Confira os resultados no CMR abaixo. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# leitura dos dados

soja <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/cnpaf/soja.txt",
                   header=TRUE, sep="\t", dec=",")
str(soja)

#-----------------------------------------------------------------------------
# ver o peso de 100 grãos

require(lattice)

xyplot(pesograo~potassio, groups=agua, data=soja, type=c("p","a"))
xyplot(pesograo~potassio|agua, data=soja, type=c("p","a"))

#-----------------------------------------------------------------------------
# ajuste de um modelo polinomial nos dois fatores

m0 <- lm(pesograo~bloco+poly(agua,2)*poly(potassio,2), data=soja)
par(mfrow=c(2,2)); plot(m0); layout(1)

anova(m0)
summary(m0)

m1 <- lm(pesograo~bloco+(agua+I(agua^2))*potassio+I(potassio^2), data=soja)
par(mfrow=c(2,2)); plot(m1); layout(1)

anova(m0, m1)
summary(m1)

#-----------------------------------------------------------------------------
# fazendo a predição

pred <- expand.grid(bloco="I", agua=seq(37.5,62.5,l=30),
                    potassio=seq(0,180,l=30))
pred$y <- predict(m1, newdata=pred)

#-----------------------------------------------------------------------------
# vendo a superfície ajustada

source("http://www.leg.ufpr.br/~walmes/ridiculas/ridiculas_functions.R")
args(panel.3d.contour)
body(panel.3d.contour)

require(RColorBrewer)

colr <- brewer.pal(11, "Spectral")
colr <- colorRampPalette(colr, space="rgb")

zlab <- "Peso de 100 grãos"
xlab <- "Potássio no solo"
ylab <- "Nível de saturação de água"

#png(file="f034.png", width=500, height=500)
wireframe(y~potassio*agua, data=pred,
          scales=list(arrows=FALSE), zlab=list(zlab, rot=90),
          xlab=list(xlab, rot=24), ylab=list(ylab, rot=-37),
          zlim=c(5,16), col="gray50", col.contour=1,
          panel.3d.wireframe="panel.3d.contour", type=c("on","bottom"),
          col.regions=colr(100),  drape=TRUE,
          screen=list(z=40, x=-70))
#dev.off()

#-----------------------------------------------------------------------------

Emacs personalizado

Printscreen da minha sessão Emacs com as atuais configurações do arquivo .emacs.

Trocando umas idéias com o Fernando Toledo (co-autor do Ridículas) aprendi diversas opções disponíveis para personalizar o Emacs, meu editor de scripts R. A personalização torna a experiência mais agradável, de forma que o usuário pode habilitar ferramentas do editor que proporcionem mais conforto e agilidade. Algumas opções de personalização são os padrões de cores. Eu por exemplo prefiro um fundo preto pois acredito que cansa menos. Numeração na margem é outro aspecto interessante, principalmente para situar pessoas que acompanham algum curso com o R usando o Emacs. O realçador de parenteses é outro recurso que evita perda de tempo na depuração de um código. E coisas como especificar o tamanho de letra e definir como padrão os comandos ctrl-zxcv para desfazer-recortar-copiar-colar evitam de você ter que fazer isso graficamente em toda sessão. Para usar o esquema de cores você vai ter que instalar o emacs-goodies-el

$ sudo apt-get install emacs-goodies-el

Uma vez instalado, o restante das opções você consegue copiando para o seu arquivo .emacs, que fica na sua home (se não existe, crie um), os trechos desejados do meu atual arquivo .emacs

;;;--------------------------------------------------------------------------
;;; arquivo de configuração do Emacs por Walmes Zeviani

;;---------------------------------------------------------------------------
;; faz com que o emacs23 fique no modo de texto visual
;; (quebra de linha visual) quando edita porção tex de um Rnw
(add-hook 'text-mode-hook 'turn-on-visual-line-mode)

;;---------------------------------------------------------------------------
;; faz com que apareceça os argumentos das funções do R no minibuffer
(require 'ess-eldoc)

;;---------------------------------------------------------------------------
;; numeração das linhas na margem esquerda
(global-linum-mode 1)

;;---------------------------------------------------------------------------
;; inicia Emacs com ctrl-{zxcv} abilitado para desf/recor/cop/colar
(cua-mode t)

;;---------------------------------------------------------------------------
;; realçador de pareamento de parenteses, chaves, colchetes, aspas...
(show-paren-mode 1)

;;---------------------------------------------------------------------------
;; tamanho da fonte (120~12pt)
;; comandos C-x C-+ e C-x C-- para aumentar e diminiur o tamanho da fonte
(set-face-attribute 'default nil :height 120)

;;---------------------------------------------------------------------------
;; estilo de cores de fundo e fontes, mais disponíveis em
;; http://themes.sweyla.com/
;; http://jasonm23.github.com/emacs-theme-editor/

(defun sweyla838549 ()
  "Theme generated by Sweyla: http://themes.sweyla.com/seed/838549/"
  (interactive)
  (color-theme-install
   '(sweyla838549
     ((background-color . "#010000")
      (foreground-color . "#FFFFFF")
      (background-mode . dark)
      (border-color . "#323232")
      (cursor-color . "#FFFFFF")
      (mouse-color . "#323232"))
     (mode-line ((t (:foreground "#FFFFFF" :background "#323232"))))
     (region ((t (:background "#323232"))))

     (font-lock-comment-face ((t (:foreground "#00FF73"))))
     (font-lock-constant-face ((t (:foreground "#58D3A8"))))
     (font-lock-builtin-face ((t (:foreground "#16AF3D"))))
     (font-lock-function-name-face ((t (:foreground "#82FFD5"))))
     (font-lock-variable-name-face ((t (:foreground "#00FFD2"))))
     (font-lock-keyword-face ((t (:foreground "#3AFFC4"))))
     (font-lock-string-face ((t (:foreground "#04F069"))))
     (font-lock-doc-string-face ((t (:foreground "#04F069"))))
     (font-lock-type-face ((t (:foreground "#00FFFF"))))
     )))

(provide 'sweyla838549)

(require 'color-theme)
(color-theme-initialize)
(sweyla838549)

;;---------------------------------------------------------------------------

Visite os dois links disponíveis na parte de cores. Nesses sites você pode definir um esquema de cores próprios para sua sessão Emacs. Daí é só aproveitar a boa experiência que o Emacs dá à edição de código R e LaTex. A figura no início do post é um printscreen do meu Emacs com essas configurações. Atualmente eu uso o Emacs23 no Ubuntu 12.04. Não esqueça para usar o R com Emacs você precisa instalar o ess

$ sudo apt-get install emacs23 ess

Meus agradecimentos ao Fernando Toledo pelas preciosas dicas. Até a próxima ridícula.

Gerar provas diferentes com Sweave

Método japonês para evitar a cola.

Há muito tempo que eu procurava uma forma de gerar provas diferentes por meio de um processo de sorteio aleatório. Eu não queria fazer duas provas, tirar cópias delas e distribuir. Eu queria uma prova diferente da outra, sem cópias, entende? Bem, no começo pensei que isso fosse algo extremamente difícil mas não é! Nos vamos fazer isso com Sweave, que é o mecanismo de acoplar o R (ambiente de computação estatística) ao LaTex (editor de textos).

O procedimento é bem simples e eu vou descrevê-lo em etapas:

Etapa 1. Criamos um arquivo Sweave (extensão *.Rnw) e colocamos os enunciados da nossa prova. Esse arquivo não precisar ter preâmbulo, apenas o conteúdo da prova. Como exemplo, considere o conteúdo abaixo que deve ser salvo em um arquivo do Sweave (prova.Rnw). Note que a função sample() é a responsável pela aleatoriedade das questões. Com isso, em cada compilação teremos enunciados diferentes, cada um deles com mesma probabilidade de ocorrer;

<<echo=false, results=hide>>=
r <- sample(1:2, 1)
conf <- list(A="da mesma cor", B="de cor diferente")
x <- conf[[r]]
@ 

\noindent 1. Pega-se um baralho e coloca-se os quatro ases na mesa, virados 
para baixo. Dois deles (A$\clubsuit$, A$\spadesuit$) são pretos, os outros 
dois (A$\heartsuit$, A$\diamondsuit$) são vermelhos. Seja o experimento
retirar duas dessas quatro cartas aleatoriamente. Qual a probabilidade das 
cartas serem \Sexpr{x}?

<<echo=false, results=hide>>=
r <- sample(1:2, 1)
conf <- list(A=c(2,2), B=c(3,2.5))
symb <- c(0,0,1,conf[[r]][1])
symbb <- paste("|", paste(symb, collapse="|"), "|", sep="")
cust <- conf[[r]][2]
@ 

\vspace{2cm}
\noindent 2. Uma máquina caça níquel possui 2 cilíndros e em cada uma deles
está gravado a sequência de números \Sexpr{symbb}. A máquina premia o jogador
com o valor da soma dos números observados, ou seja, se sair | 1 | 0 | o
jogador ganha 1 dólar. O valor de cada jogada é \Sexpr{cust} dólares. Seja
$X$ o \textbf{lucro} do jogador em cada jogada. Determine:
\begin{compactenum}[a)]
  \item A distribuição de probabilidades de $X$;
  \item O lucro médio do jogador, ou seja, o valor esperado de $X$.
\end{compactenum}

Etapa 2. Criamos um arquivo R para fazermos a compilação repetida do arquivo prova.Rnw. Cada resultado de compilação será salvo em um diretório dentro ao atual (./cada_tex) com nome de arquivos diferentes (prova1.tex, prova2.tex, etc). Nessa mesma etapa nos geramos um arquivo todos_inputs.tex que nada mais é que a chamada de cada prova para compilação que será usada no arquivo *.tex da etapa a seguir;

# cria o diretório para os arquivos compilados
dir.create("cada_tex")
# número de provas e o nome das provas
n <- 10
names <- paste("cada_tex/prova", 1:n, ".tex", sep="")
# função para compilar com nomes de saída diferente
do.tex <- function(name){
  Sweave("prova.Rnw", encoding="utf8", output=name)
}
sapply(names, do.tex)
# para gerar um arquivo tex com os inputs
sink("todos_inputs.tex")
for(i in names){
  cat(paste("\\input{", i, "}\n", sep=""))
  cat("\\newpage\n")
}
sink()
# para compilar o arquivo tex
system("pdflatex provas.tex")

Etapa 3. Agora nos compilamos de *.Rnw para *.pdf. Esse arquivo (provas.tex) contém o preâmbulo que precisamos e o input de cada uma das provas, uma por página.

\documentclass[a4paper,10pt]{article}
\usepackage{Sweave}
\usepackage[utf8]{inputenc}
\usepackage[brazil]{babel}
\usepackage[T1]{fontenc}
\usepackage{enumerate}
\usepackage{paralist}

\begin{document}
\pagestyle{empty}

\input{todos_inputs.tex}

\end{document}

Ao executar essas etapa no seu computador, verifique a codificação dos arquivos. Como eu uso linux, minha codificação é a padrão do SO (utf-8). No mais é só aproveitar. Até a próxima ridícula.

Ajuste simultâneo de regressão/anova para várias respostas

Valores observados e ajustados acompanhados do intervalo de 95% em função da dose de potássio e nível de água para as variáveis resposta estudadas.

Na realização de experimento é comum a observação de várias variáveis respostas em uma mesma unidade amostral. Em termos estatísticos, o que observamos é um vetor aleatório. Por exemplo, em experimento de com soja para fins de ensaio de competição de cultivares o pesquisador pode observar: altura da planta, altura da primeira vagem, massa seca da planta, número de vagens por planta, teores foliares de N, P, K, Mg, Ca, peso de 100 grãos, proporção de sementes viáveis, dias para fechar o ciclo, todas elas observadas em cada unidade experimental. Nesse caso então, observamos um vetor aleatório de 12 variáveis aleatórias. O ideal seria empregarmos um modelo multivariado, porém como as respostas podem ser contínuas não limitadas, contagens, contínuas limitadas, contínuas com censura, proporções, etc, ainda não dispomos de um modelo de distribuição multivariado que contemple essa situação. O mais comum é fazer a análise univariada e como todas essas variáveis estão sob o efeito das mesmas fontes de variação, a análise estatística vai empregar o mesmo modelo (pelo menos no tocante a parte determinística). Sendo assim, tem como fazer o estudo de todas as variáveis reposta de uma vez só? Sim. Para saber como continue lendo.

Existem diversas formas de aplicar tarefas em grupos ou colunas. A família apply do stats oferece recurso para um número bem grande de situações. A família ply do plyr acrescenta mais opções. Além destes temos o doBy com mais algumas funcionalidades.

Bem, vamos ao que interessa. No CMR eu faço ajustes de regressão com polinômios para dados de um experimento com a cultura da soja (vasos em casa de vegetação). Foram estudados o nível de água (umidade do solo em relação à capacidade de campo) e potássio aplicado ao solo. O objetivo desse experimento foi verificar se o fornecimento de potássio pode minimizar impactos do deficit hídrico. Para acessar o artigo completo clique aqui: Umidade do solo e doses de potássio na cultura da soja. Não deixe de conferir a lista de links abaixo. Em próximos artigos do Rídiculas vamos aprofundar o estudo nesses dados comparando curvas. Até lá.

A brief introduction to “apply” in R;

plyr – The split-apply-combine strategy for R;

plyr: Tools for splitting, applying and combining data;

R Grouping functions: sapply vs. lapply vs. apply. vs. tapply vs. by vs. aggregate vs.;

Say it in R with “by”, “apply” and friends;

The doBy package;

doBy: doBy – Groupwise summary statistics, general linear contrasts, population means (least-squares-means), and other utilities;

#-----------------------------------------------------------------------------
# importa dados
soja <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/soja.txt",
                   header=TRUE, sep="\t", dec=",")
str(soja)
soja <- transform(soja, a=agua, A=factor(agua),
                  k=potassio, K=factor(potassio))
resp <- 4:7 # índice das variáveis resposta
str(soja)

#-----------------------------------------------------------------------------
# ver

soja2 <- cbind(stack(soja[,resp]), soja[,-resp])

require(lattice)

xyplot(values~k|ind, groups=A, data=soja2,
       scales="free", jitter.x=TRUE, type=c("p","a"))

#-----------------------------------------------------------------------------
# ajuste do modelo saturado de efeitos categóricos (modelo de médias)

m0 <- apply(soja[,resp], 2,
            function(r){ aov(r~bloco+A*K, soja) })
lapply(m0, anova)

#-----------------------------------------------------------------------------
# ajuste de polinômio para potassio por nível de água

m1 <- apply(soja[,resp], 2,
            function(r){ aov(r~bloco+A/poly(k, 4), soja) })
lapply(m1, anova)
names(coef(m1[[1]]))

inter <- m1[[1]]$assign==3        # indices das interações
eff <- names(coef(m1[[1]]))[inter]
togrep <- outer(c("A37.5","A50","A62.5"), c("1$","2$","(3|4)$"),
                paste, sep=":.*") # nomes de busca
or <- order(togrep)
tolist <- outer(c("A37.5","A50.0","A62.5"), c("lin","qua","lof"),
                paste, sep="|")   # nomes de rótulo
list.desd <- sapply(togrep, grep, eff)
names(list.desd) <- tolist
list.desd <- list.desd[or]
list.desd                         # lista com os índices dos desdobramentos

# anovas com desdobramento de efeito de potássio em cada nível de água
lapply(m1, summary, split=list("A:poly(k, 4)"=list.desd))

#-----------------------------------------------------------------------------
# com base nas anovas acima selecionamos os modelos

form <- list( rengrao~bloco+A/poly(k, 2, raw=TRUE),
             pesograo~bloco+A/poly(k, 2, raw=TRUE),
                kgrao~bloco+A/poly(k, 2, raw=TRUE),
                pgrao~bloco+A/poly(k, 3, raw=TRUE))
names(form) <- names(soja)[resp]

m2 <- lapply(form, aov, data=soja, contrasts=list(bloco=contr.sum))
lapply(m2, coef)

#-----------------------------------------------------------------------------
# predição a partir do modelo

str(soja)
pred <- expand.grid(bloco="I", A=levels(soja$A), k=seq(0,180,5))

#aux <- lapply(m2, predict, newdata=pred, interval="confidence")
aux <- lapply(m2, function(r){ # remove o efeito do bloco na predição
  predict(r, newdata=pred, interval="confidence")-coef(r)["bloco1"]})
aux <- lapply(aux, as.data.frame)
str(aux)

require(plyr)
pred <- cbind(ldply(aux), pred)
pred$.id <- factor(pred$.id, levels=levels(soja2$ind))

#-----------------------------------------------------------------------------
# gráfico final

fl <- c("K nos grãos (mg)", "Peso 100 grãos (g)",
        "P nos grãos (mg)", "Rendimento de grãos (g)")
cbind(levels(pred$.id), levels(soja2$ind), fl)

require(latticeExtra)
display.brewer.all()
mycol <- brewer.pal(6, "Set1")

panel.ciH <- function(x, y, ly, uy, subscripts, ...){
  y <- as.numeric(y)
  x <- as.numeric(x)
  or <- order(x)
  ly <- as.numeric(ly[subscripts])
  uy <- as.numeric(uy[subscripts])
  panel.polygon(c(x[or], x[rev(or)]), c(ly[or], uy[rev(or)]),
                col=1, alpha=0.05, border=NA)
  panel.lines(x[or], ly[or], lty=3, lwd=0.5, col=1)
  panel.lines(x[or], uy[or], lty=3, lwd=0.5, col=1)
  panel.xyplot(x, y, subscripts=subscripts, ...)
}

#png(file="f031.png", width=500, height=500)
trellis.par.set(superpose.symbol=list(col=mycol),
                superpose.line=list(col=mycol),
                strip.background=list(col="gray90"))
xyplot(values~k|ind, groups=A, data=soja2,
       scales="free", jitter.x=TRUE,
       xlab="Potássio no solo", ylab=NULL,
       strip=strip.custom(factor.levels=fl),
       auto.key=list(title="Níveis de água (%)", cex.title=1.2,
         columns=3, points=FALSE, lines=TRUE))+
  as.layer(xyplot(fit~k|.id, groups=A, pred, type="l", lwd=1.5,
                  ly=pred$lwr, uy=pred$upr,
                  panel.groups=panel.ciH,
                  panel=panel.superpose))
#dev.off()

#-----------------------------------------------------------------------------
# pode-se tentar comparar as curvas ajustando modelos apropriados sob H0
# exemplo, será que as curvas vermelha e verde para kgrao e pesograo não
# podem ser representados pela mesma curva?
# as curvas azul e ver para pgrao diferem de uma reta horizontal?

#-----------------------------------------------------------------------------

Reparametrização do modelo quadrático para ponto crítico

Perfil de log-verossimilhaça representada pelo módulo da raíz da função deviance para cada um dos parâmetros do modelo quadrático reparametrizado. Linhas tracejadas indicam os limites dos intervalos de confiança de 99, 95, 90, 80 e 50%. Padrões de perfil diferentes de um V indicam intervalos assimétricos.

O emprego de funções polinomiais é muito frequente em ciências agrárias. O polinômio de segundo grau ou modelo quadrático é a escolha mais natural e imediata de modelo para quando a variável dependente apresentam uma tendência/relação não linear com a variável independente. Não raramente, o emprego desse modelo motiva também a estimação do ponto crítico (mínimo ou máximo). Para isso, o que se faz é estimar os parâmetros do modelo e via regras do cálculo obter a estimativa pontual do ponto crítico. O problema é que normalmente se necessita e não se faz, por desconhecimento de métodos, associação de incerteza à essa estimativa. É sobre isso que vou tratar agora: como associar incerteza ao ponto crítico de um modelo quadrático.

O modelo quadrático tem essa equação

,

em que  é o intercepto, ou seja, ,  é valor da derivada do modelo na origem, ou seja, , e mede a curvatura. Derivando o modelo, igualando a zero e resolvendo, obtemos que o ponto crítico é

,

portanto,  é uma função das estimativas dos parâmetros, e como tal, deve apresentar incerteza associada.

Uma forma de associar incerteza é empregar o método delta para encontrar erros padrões. No R, o mesmo está disponível em msm::deltamethod(). Outra forma é usar um simulação boostratp e aplicar nas amostras bootstrap a relação acima. E o terceiro, e que eu mais gosto pois tem uma série de vantagens, é usar o perfil de log-verossimilhança.

Para usar o perfil de verossimilhança eu tive que escrever o modelo de forma que o ponto crítico fosse um parâmetro dele. Reparametrizar um modelo nem sempre é uma tarefa fácil. Nesse caso o procedimento escrito é curto e simples de entender. Pórem eu demorei cerca de uma hora com tentativa-erro até usar a seguinte abordagem: imagine uma função que tenha ponto de máximo. Podemos aproximar essa função ao redor de um valor por série de Taylor. Para uma aproximação de segundo grau temos

.

Agora imagine que é um ponto crítico da função. O termo linear da aproximação é 0 no ponto crítico por definição. Se a função que estamos aproximando for perfeitamente quadrática, esse procedimento nos leva a reparametrização que desejamos. Assim, o modelo quadrático reparametrizado é

.

em que é o nosso parâmetro de interesse, o ponto crítico, é o valor da função no ponto crítico, ou seja, , é a métade da segunda derivada da função no ponto crítico, ou seja, , é uma curvatura. Então nos reparametrizamos o model quadrático para ter o ponto crítico como parâmetro. Nessa reparametrização obtvemos 3 parâmetros com interpretação simples. Em termos de interpretação, se não temos curvatura, portanto, não temos ponto crítico e devemos ajustar um modelo linear da reta. Se temos curvatura podemos estimar o ponto crítico e o valor da função no ponto crítico. Então imagine que você estuda doses de nutrientes para uma cultura, não é fundamental saber qual a dose que confere o máximo e qual a produção esperada nessa dose? Melhor que isso é poder associar incerteza à essas quantidades.

O modelo reparametrizado, diferente do original, é um modelo não linear. Por isso, para estimar os parâmetros no R usaremos a função nls() e não a lm(). (Visite os posts passados para ver mais coisas relacionadas à modelos não lineares.) Para obter os intervalos de perfil de log-verossimilhança usaremos os métodos confint() e profile(). Consulte a documentação dessas funções para saber o que elas de fato fazem.

No CMR abaixo eu ilustro a aplicação desses métodos. É recomendando ao leitor que compare os intervalos de confiança, que repita os procedimento simulando os dados com outros parâmetros, tamanhos e espaços, e se for o caso, que use dados reais. Futuros artigos do Rídiculas vão incluir comparação de modelos quadráticos (teste de hipótese) e o modelo quadrático para duas regressoras (bidimensional). Até a próxima rídicula.

#-----------------------------------------------------------------------------
# dados simulados

set.seed(23456)
x <- 1:12
y <- 3+0.17*x-0.01*x^2+rnorm(x,0,0.05)
plot(x, y)                         # dados observados
curve(3+0.17*x-0.01*x^2, add=TRUE) # curva que gerou os dados

#-----------------------------------------------------------------------------
# ajuste do y = a+b*x+c*x², é um modelo linear

m0 <- lm(y~x+I(x^2))
summary(m0)

#-----------------------------------------------------------------------------
# método delta

require(msm)
coef <- coef(m0)[2:3]; vcov <- vcov(m0)[2:3,2:3]
xm <- -0.5*coef[1]/coef[2]; xm                         # estimativa
sd.xm <- deltamethod(~(-0.5*x1/x2), coef, vcov); sd.xm # erro padrão
xm+c(-1,1)*qt(0.975, df=df.residual(m0))*sd.xm         # IC 95%

#-----------------------------------------------------------------------------
# simulação (bootstrap paramétrico)

require(MASS)
aa <- mvrnorm(1000, mu=coef, Sigma=vcov)
xm.boot <- apply(aa, 1, function(i) -0.5*i[1]/i[2])
mean(xm.boot)                      # média bootstrap
sd(xm.boot)                        # desvio-padrão bootstrap
quantile(xm.boot, c(0.025, 0.975)) # IC boostrap

#-----------------------------------------------------------------------------
# ajuste do y = ym+c*(x-xm)², é um modelo não linear

plot(x, y)
abline(h=3.7, v=8.5, lty=3) # aproximação visual do chute
start <- list(ym=3.7, xm=8.5, c=coef(m0)[3])
max(fitted(m0))             # boa aproximação de chute
x[which.max(fitted(m0))]    # boa aproximação de chute
n0 <- nls(y~ym+c*(x-xm)^2, start=start)
summary(n0) # estimativa do xm que é parâmetro do modelo, estimado diretamente

#-----------------------------------------------------------------------------
# gráfico dos observados com a curva ajustada

plot(x, y)
with(as.list(coef(n0)),{
  curve(ym+c*(x-xm)^2, add=TRUE)
  abline(h=ym, v=xm, lty=3)
})

#-----------------------------------------------------------------------------
# veja a igualdade dos modelos

cbind(coef(m0),
      with(as.list(coef(n0)), c(ym+c*xm^2, -2*c*xm, c)))
deviance(m0)
deviance(n0)
plot(fitted(m0), fitted(n0)); abline(0,1)

#-----------------------------------------------------------------------------
# intervalos de confiança

confint.default(n0) # IC assintóticos, são simétricos por construção
confint(n0)         # IC baseado em perfil de log-verossimilhança

#-----------------------------------------------------------------------------
# gráficos de perfil de verossimilhança

#png(file="f030.png", width=500, height=175)
par(mfrow=c(1,3))
plot(profile(n0))
layout(1)
#dev.off()

#-----------------------------------------------------------------------------
# representação do ajuste

pred <- data.frame(x=seq(0,13,l=30))
pred <- cbind(as.data.frame(predict(m0, newdata=pred, interval="confidence")),
              pred)
str(pred)
ic <- confint(n0); ic

plot(x, y, type="n")              # janela vazia
with(pred,                        # bandas de confiança
     polygon(c(x,rev(x)), c(lwr,rev(upr)), col="gray90", border=NA))
points(x, y)                      # valores observados
with(pred, lines(x, fit))         # valores preditor
abline(v=ic[2,], h=ic[1,], lty=2) # IC para xm e ym

#-----------------------------------------------------------------------------