Crescimento micelial por análise de imagens

Tamanho da colônia micelial (cm²) determinada pela análise de imagem com funções do pacote EBImage.

Nessa matéria vou apresentar os procedimentos em R para determinar o tamanho de uma colônia de micélios em placa de petri. Os resultados são baseados em numa calibração que consistiu em escolher melhor fundo, iluminação e demais parâmetros fotográficos para uma boa recuperação da informação via análise de imagem pelo pacote EBImage. Em outras palavras, como obter uma boa foto para cálculo do crescimento micelial. Esta foi uma epata preliminar de um estudo que vai comparar a determinação do crescimento micelial por análise de imagens e pelo método padrão, via uso de paquímetros.

As fotos foram obtidas com uma câmera fotográfica Canon Ti preza em um suporte de madeira com distância fixa da placa de petri. A iluminção foi ambiente. Fotografamos a placa de petri aberta. Inicialmente nós usamos um fundo azul mas que não foi adequado por não apresentar diferenças em tons na escala cinza da região micelial. Então substituímos por um fundo preto para aumentar esse contraste de tons. Superado isso, tivemos problemas com o brilho da borda da placa de petri. Isso foi resolvido com seleção de pixels dentro de um círculo que não contivesse tal borda. Vai ficar mais claro no CMR abaixo cada etapa do processo. Por fim, determinamos a área, perímetro, e diâmetros da colônia com as funções do pacote EBImage.

Por fim, tenho que deixar os créditos ao acadêmico de Doutorado em Agronomia (UFPR) Paulo Lichtemberg que é o responsável pelo estudo de crescimento de fungos por análise de imagens, membro do LEMID (Laboratório de Epidemiologia e Manejo Integrado de Doenças) e orientado da Professora Larissa May De Mio.

#-----------------------------------------------------------------------------
# carrega o pacote

require(EBImage)

#-----------------------------------------------------------------------------
# lendo o arquivo

# lê imagem, imagem original com resolução de 2352x1568 foi recortada e
# reduzida para 400x400 pixel de dimensão

f0 <- readImage("http://www.leg.ufpr.br/~walmes/ridiculas/micelio.JPG")
str(f0)

display(f0) # vê a imagem
hist(f0)    # histograma dos componentes verde, vermelho e azul

par(mfrow=c(3,1)) # gráfico de densidade do vermelho, verde e azul
apply(f0, MARGIN=3, function(x) plot(density(x), xlim=c(0,1)))
layout(1)

#-----------------------------------------------------------------------------
# tratamento para escala cinza

f1 <- imageData(channel(f0, mode="red")) # vermelho parece separar melhor
f1 <- 1-f1                               # inverte as tonalidades
plot(density(f1), xlim=c(0,1))
b <- 0.47
abline(v=b)

filled.contour(f1, asp=1)
display(f1) # escala cinza com claro sendo a folha

#-----------------------------------------------------------------------------
# dicotomiza para branco e preto

f2 <- f1
f2[f1<b] <- 1
f2[f1>=b] <- 0
display(f2) # ops! temos a borda da placa presente, removê-la!
            # selecionar só pixels dentro de um círculo que exclui tal borda

# usado para selecionar pontos detro um círculo de certo raio
mx <- nrow(f1)/2; my <- ncol(f1)/2 # centro da imagem

# matriz de distâncias de cada pixel a partir do centro
M <- outer(1:nrow(f1), 1:ncol(f1),
           function(i,j) sqrt((i-mx)^2+(j-my)^2))
str(M)

f2 <- f1
f2[f1<b] <- 1
f2[f1>=b] <- 0
f2[M>155] <- 0
display(f2) # ok! região do micélio

f3 <- f1
f3[f1<0.65] <- 1
f3[f1>=0.65] <- 0
display(f3) # região do interior da placa de petri

#-----------------------------------------------------------------------------
# tratamento que remove pontos pretos dentro das regiões brancas

f2 <- bwlabel(f2)  # identifica os conjuntos brancos
f2 <- fillHull(f2) # elimina pontos pretos dentro do branco
f3 <- fillHull(f3)

display(f2)
display(f3)

kern <- makeBrush(3, shape="disc", step=FALSE) # remove 1 px

f2 <- erode(f2, kern) # remove alguns ruídos
f3 <- erode(f3, kern)

display(f2)
display(f3)

#-----------------------------------------------------------------------------
# calcula as dimensões e converte para cm pois o diâmetro da placa é 10cm

dimen <- c(area=pi*5^2, peri=2*pi*5, 5, 5, 5) # em cm
micel <- computeFeatures.shape(f2)
placa <- computeFeatures.shape(f3)

cm <- micel*dimen/placa
cm # área em cm² e demais em cm

#-----------------------------------------------------------------------------
# prepara para exportação (escolher preenchido ou borda)

f4 <- paintObjects(f2, f0, opac=c(NA, 0.45), col=c(NA, "red")) # preenchido
f4 <- paintObjects(f2, f0, opac=c(1, NA), col=c("black", NA))  # borda
display(f4)

xy <- computeFeatures.moment(f2, f0)[, c("m.cx", "m.cy")] # centro de massa
font <- drawfont(weight=600, size=15)
f5 <- drawtext(f4, xy=xy,
               labels=paste(format(cm[,"s.area"], digits=4), "cm²"),
               font=font, col="black")
display(f5)

writeImage(f5, "f042.jpg")

#-----------------------------------------------------------------------------
Categorias:dados Tags:

Legenda e contornos em gráficos de nível

Peso de 100 grãos em função do nível de saturação de água e potássio. Contornos e rótulo na escala de cores são os destaques desse gráfico.

Mais uma ridícula de lattice. Fiz duas coisas: adicionei rótulo à legenda de cores (colorkey) e coloquei contornos sobre o gráfico de níveis. Assim como a maioria das dicas gráficas, essa também é baseada nas mensagens da r-help. Os respectivos links estão no CMR. 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)

summary(influence.measures(m0))
soja[55,] # influente segundo dffits

m0 <- lm(pesograo~bloco+poly(agua,2)*poly(potassio,2),
         data=soja[-55,])
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[-55,])
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)

#-----------------------------------------------------------------------------
# representando com wireframe()

require(RColorBrewer)

display.brewer.all()
colr <- brewer.pal(11, "RdYlGn")
colr <- colorRampPalette(colr, space="rgb")

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

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),
          col.regions=colr(100),  drape=TRUE,
          screen=list(z=40, x=-70))

#-----------------------------------------------------------------------------
# representando em um levelplot()

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

levelplot(y~potassio*agua, data=pred, col.regions=colr(100),
          xlab=xlab, ylab=ylab)

#-----------------------------------------------------------------------------
# adicionando rotulo à legenda de cores, baseado nas mensagens da r-help
# http://r.789695.n4.nabble.com/Adding-title-to-colorkey-td4633584.html

library(grid)

# modo 1
levelplot(y~potassio*agua, data=pred, col.regions=colr(100),
          xlab=xlab, ylab=ylab,
          par.settings=list(layout.widths=list(axis.key.padding=4)))
grid.text(zlab, x=unit(0.88, "npc"), y=unit(0.5, "npc"), rot=90)

# modo 2
levelplot(y~potassio*agua, data=pred, col.regions=colr(100),
          xlab=xlab, ylab=ylab,
          ylab.right=zlab,
          par.settings=list(
            layout.widths=list(axis.key.padding=0, ylab.right=2)))

require(latticeExtra)

# modo 3
p <- levelplot(y~potassio*agua, data=pred, col.regions=colr(100),
               xlab=xlab, ylab=ylab,
               par.settings=list(
                 layout.widths=list(right.padding=4)))
p$legend$right <- list(fun=mergedTrellisLegendGrob(p$legend$right,
                         list(fun=textGrob, args=list(zlab, rot=-90, x=2)),
                         vertical=FALSE))
print(p)

#-----------------------------------------------------------------------------
# adicionando contornos, baseado em
# https://stat.ethz.ch/pipermail/r-help/2006-February/088166.html

#png(file="f037.png", width=500, height=400)
p <- levelplot(y~potassio*agua, data=pred, col.regions=colr(100),
               xlab=xlab, ylab=ylab,
               panel=function(..., at, contour=FALSE, labels=NULL){
                 panel.levelplot(..., at=at, contour=contour, labels=labels)
                 panel.contourplot(..., at=at, contour=TRUE,
                                   labels=list(labels=format(at, digits=4),
                                     cex=0.9))
               },
               par.settings=list(
                 layout.widths=list(right.padding=4)))
p$legend$right <- list(fun=mergedTrellisLegendGrob(p$legend$right,
                         list(fun=textGrob, args=list(zlab, rot=-90, x=2)),
                         vertical=FALSE))
print(p)
#dev.off()

#-----------------------------------------------------------------------------
Categorias:gráficos Tags:,

Cálculo de área foliar ocupada por cochonilhas

Área foliar (limão) ocupada por cochonilhas (%). As determinações de área foram feitas com funções disponíveis no pacote EBImage.

No post passado (Cálculo de área foliar) eu mostrei como calcular área foliar com o pacote EBImage do R, desenvolvido por Oleg Sklyar, Gregoire Pau, Mike Smith e Wolfgang Huber. O segundo mencionado é o mantenedor do pacote. Atualmente estou fazendo testes com o pacote porque considero promissor o seu uso em experimentos agronômicos, mais precisamente àqueles de fitopatologia, entomologia e solos. O R por ser livre pode ser usado em todas as etapas do experimento, a saber, a de aquisição dos dados considerando o processamento das imagens e de processamento de dados considerando a análise estatística. Dessa forma, alguém pode desenvolver procedimentos para análise de doenças em escalas diagramáticas avaliadas por computador, eliminando o traço de subjetividade e não reproducibilidade que existe nas notas dadas por um avaliador. E sei que os avaliadores são treinados para tal tarefa, não é qualidade desse trabalho que estou pontuando. O que considero desvantajoso nas escalas é que elas são sempre com poucos níveis (unitária de 0 à 5, por exemplo). Isso implica numa aproximação grosseira de métodos de análise de dados, principalmente porque os métodos de inferência frequentemente aplicados à esses casos supõem normalidade. A avaliação feita pelo computador gera medidas contínuas e isso é bom.

Aqui eu determino a área ocupada por cochonilha em duas folhas de limoeiro (uma planta que tenho no quintal de casa). Eu escolhi uma folha muito ocupada e outra pouco ocupada. O CMR abaixo apresenta todos os passos. Na figura que acompanha essa matéria você vê o resultado final das análises. Quando fiz minha avaliação visual das folhas eu chutei 5% e 40% de ocupação e confesso que fiquei surpreso com os resultados (errei feio). De fato, não temos boa capacidade de avaliar área em situações de distribuição pontual como essa. Por isso que acredito que uma avaliação de área pelo computador seria melhor que a de avaliadores treinados. Alguém se arrisca a submeter isso ao teste? Outro ponto importante nessa minha análise de brincadeira é que só foram consideradas as cochonilhas brancas, não me preocupei com as vermelhas. Algum esforço seria necessário para considerá-las, ou seja, para juntar o vermelho ao branco e separar do verde. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# carrega o pacote

require(EBImage)

#-----------------------------------------------------------------------------
# lê imagem com folhas digitalizadas

fol <- readImage("http://www.leg.ufpr.br/~walmes/ridiculas/folhas3.jpg")
str(fol)

display(fol) # vê a imagem
hist(fol)    # histograma dos componentes verde, vermelho e azul

#-----------------------------------------------------------------------------
# tratamento

fol2 <- imageData(channel(fol, mode="blue")) # seleciona um canal
fol2 <- 1-fol2                               # inverte as tonalidades
hist(fol2)                                   # histograma dos tons de cinza

display(fol2) # escala cinza com claro sendo a folha

#-----------------------------------------------------------------------------
# dicotomiza para branco e preto

fol2[fol2<0.5] <- 0
fol2[fol2>=0.5] <- 1

display(fol2)

#-----------------------------------------------------------------------------
# mais tratamento

fol3 <- bwlabel(fol2)     # coloca os rótulos nas regiões disjuntas
fol4 <- fillHull(fol3)    # remove os pontos pretos dentro do branco

display(fol3)
display(fol4)

#-----------------------------------------------------------------------------
# cálculos

f.ocu <- computeFeatures.shape(fol3) # folha ocupada
f.tot <- computeFeatures.shape(fol4) # folha total

p <- f.ocu[, "s.area"]/f.tot[, "s.area"]
p <- p[p<1]
p <- 1-p
100*p # áreas ocupadas nas folhas (%)

#-----------------------------------------------------------------------------
# prepara para exportação

fol5 <- paintObjects(fol3, fol, opac=c(NA, 0.45), col=c(NA, "blue"))
display(fol5)

xy <- computeFeatures.moment(fol3, fol)[, c("m.cx", "m.cy")] # centróides
font <- drawfont(weight=600, size=15)
fol5 <- drawtext(fol5, xy=xy, labels=paste(format(100*p, digits=3), "%"),
                 font=font, col="yellow")
display(fol5)

writeImage(fol5, "f039.jpg")

#-----------------------------------------------------------------------------
Categorias:dados Tags:,

Cálculo de área foliar

Áreas foliares (cm²)

Área foliar calculada com funções do pacote EBImage.

Sinceramente, eu vi no r-bloggers e não acreditei. É possível calcular área de figuras geométricas no R. Bem isso sempre foi possível, basta ter as coordenadas do polígono. Mas não estou falando de áreas geográficas das quais possuímos os mapas. Estou falando de áreas de figuras digitalizadas, como por exemplo a folha de uma planta, as asas de uma borboleta, um torrão de solo, uma pedra, a seção de de um tronco, uma semente. Enfim, as possibilidades são infinitas.

Em agronomia, a análise de imagens é algo que vem se tornando mais comum. Determinações de volume/comprimento/diâmetro de raízes, dimensão de agregados do solo (diâmetros, rendondezas, perimentros), porcentagem de área ocupada por doença ou atacada por inseto são alguns exemplos de aplicação. Alguns aplicativos para essas tarefas são disponíveis. Alguns deles são pagos. E todos eles não permitem automatizar o trabalho, pois requerem em alguma altura do processo, intervenção do usuário via mouse. Imagine ter que calcular a área de 3000 folhas fazendo o trabalho uma a uma? Não dá né?

Nessa matéria eu calculo área de algumas folhas que encontrei no chão, embaixo de algumas árvores. Para ter a medidas das folhas em cm² eu coloquei um quadrado de papel de área conhecida como referência. Nos vamos usar funções disponíveis no pacote EBImage. O meu tutorial não é muito diferente do original que me motivou (leaf area measuring — R package “EBImage”). O fato é que não me contentei em apenas ler mencionada matéria, tive que ver com meus próprios olhos, e já que foi assim, está qui o código que produzi. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# página de desenvolvimento do pacote
# http://www.bioconductor.org/packages/devel/bioc/html/EBImage.html

# instalação no linux, no terminal do linux, fazer
# sudo apt-get install libmagickcore-dev libmagickwand-dev

# instalando do bioconductor
source("http://bioconductor.org/biocLite.R")
biocLite("EBImage") # pacote EBImage, permite determinar área foliar

# instalando do tar.gz, pegar o link da página e rodar no terminal
# R CMD INSTALL EBImage_xxxxx.tar.gz # xxx representa a versão

#-----------------------------------------------------------------------------
# carrega o pacote

require(EBImage)

#-----------------------------------------------------------------------------
# lendo o arquivo

# lê imagem com folhas digitalizadas
fol <- readImage("http://www.leg.ufpr.br/~walmes/ridiculas/folhas.jpg")
str(fol)

display(fol) # vê a imagem
hist(fol)    # histograma dos componentes verde, vermelho e azul
             # picos de azul estão mais afastados, separam melhor

#-----------------------------------------------------------------------------
# tratamento

fol2 <- imageData(channel(fol, mode="blue")) # seleciona um canal
fol2 <- 1-fol2                               # inverte as tonalidades
hist(fol2)                                   # histigrama dos tons de cinza

display(fol2) # escala cinza com claro sendo a folha

#-----------------------------------------------------------------------------
# dicotomiza para branco e preto

fol2[fol2<0.5] <- 0
fol2[fol2>=0.5] <- 1

display(fol2)

#-----------------------------------------------------------------------------
# calcula atributos de cada região

fol3 <- bwlabel(fol2)     # coloca os rótulos nas regiões disjuntas
kern <- makeBrush(3, shape="disc", step=FALSE)
fol3 <- erode(fol3, kern) # remove alguns ruídos

display(fol3)

forma <- computeFeatures.shape(fol3)
area <- forma[, "s.area"]
area

# áreas foliares
areacm <- 25*area/min(area) # quadrado de área conhecida 25 cm²

#-----------------------------------------------------------------------------
# prepara para exportação

fol4 <- paintObjects(fol3, fol, opac=c(NA, 0.45), col=c(NA, "red"))
display(fol4)

xy <- computeFeatures.moment(fol3, fol)[, c("m.cx", "m.cy")] # centróides
font <- drawfont(weight=600, size=15)
fol5 <- drawtext(fol4, xy=xy, labels=paste(format(areacm, digits=4), "cm²"),
                 font=font, col="white")
display(fol5)

writeImage(fol5, "f038.jpg")

#-----------------------------------------------------------------------------
Categorias:dados Tags:,

Gráfico de densidade com destaque dos quantis

Gráficos de densidade com regiões destacando os quartis.

Gráficos de densidade são utilizados para se conhecer a forma da distribuição dos dados. A função lattice::densityplot() faz esse tipo de gráfico. Para aprimorar a visualização dos quartis dos dados, eu implementei uma função (my.densityplot()) para que destaque com regiões coloridas as porções de dados delimitadas pelos quartis. O código abaixo gera a figura no início desse post. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# densityplot() com regiões pintadas delimitadas pelos quantis

require(lattice)

my.densityplot <- function(x, cols="gray90", probs=c(1,3)/4, ...){
  panel.densityplot(x, ...)
  dx <- density(x)
  breaks <- sort(c(range(dx$x), quantile(x, probs=probs)))
  fx <- approxfun(dx$x, dx$y)
  do.polygon <- function(x){
    y <- fx(x)
    return(list(x=c(min(x), x, max(x)), y=c(0, y, 0)))
  }
  seqs <- lapply(1:(length(breaks)-1),
                 function(i){
                   x <- seq(breaks[i], breaks[i+1], l=20)
                   do.polygon(x)
                 })
  for(i in 1:length(seqs)){
    seqs[[i]]$col <- cols[i]
  }
  lapply(seqs, function(i) do.call(panel.polygon, i))
  panel.mathdensity(dmath=dnorm, col="black",
                    args=list(mean=mean(x),sd=sd(x)), lty=3)
}

require(RColorBrewer)
display.brewer.all()

cols <- brewer.pal(3, "BuPu") # Greens, Blues, BuPu, Reds
densityplot(~height|voice.part, data=singer,
            cols=cols, plot.points="rug", panel=my.densityplot)

#png("f036.png", width=500, height=320)
cols <- brewer.pal(4, "Greens") # Greens, Blues, BuPu, Reds
densityplot(~height|voice.part, data=singer, layout=c(4,2),
            strip=strip.custom(bg="gray90"),
            xlab="Altura", ylab="Densidade", col=1,
            cols=cols, probs=1:3/4, plot.points="rug", panel=my.densityplot)
#dev.off()

#-----------------------------------------------------------------------------
Categorias:gráficos Tags:,

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 (m) 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 m 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 m 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 m 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()

#-----------------------------------------------------------------------------
Categorias:gráficos Tags:, ,