Arquivo

Arquivo do Autor

Interface gráfica para ajuste de modelos de regressão não linear

Ajuste do modelo van Genuchten à curva de retenção de água do solo com interface gráfica desenvolvida com as ferramentas do pacote rpanel.

Modelos de regressão não linear são considerados quando alguma informação a priori existe sobre o fenômeno. Essa informação pode ser, por exemplo, que a curva seja sempre crescente, típico para curvas de crescimento/acúmulo. Em geral, esses modelos têm parâmetros com interpretação física/química/biológica e alguns parâmetros não têm, mas estão presentes para conferir flexibilidade.

Talvez um dos grandes gargalos no ajuste de modelos não lineares seja atribuição de valores iniciais para o método numérico de estimação. Para parâmetros com interpretação é possível obter valores ao ver gráficos de dispersão das variáveis ou pelo menos se tem idéia do intervalo em que o parâmetro se encontra. Para parâmetros sem interpretação a obtenção de valores iniciais é realmente trabalhosa.

Um procedimento muito útil para obter valores iniciais é fazer um grid de valores iniciais e plotar a curva corresponde para cada conjunto sobre o diagrama de dispersão dos dados. Você usa como valores iniciais os correspondentes à curva que melhor se aproxima do comportamento dos dados. Por outro lado, criar um grid de valores, sobrepor cada curva ao gráfico, escolher à curva que melhor se aproxima aos dados, fornecer os valores iniciais e estimar os parâmetros é algo de demanda tempo se o número de curvas à ajustar for grande. Seria de grande ajuda se esse processo pudesse ser automatizado ou facilitado, por exemplo, por meio de uma interface gráfica com botões de ação e controle de parâmetros por meio do mouse.

O pacote rpanel possui um conjunto de funções básicas que permitem construir interfaces gráficas para o usuário interagir com o R. O conjunto de funções compreende desde caixas de texto, caixas de seleção, deslizadores, até botões para disparo de alguma função. Usando essas ferramentas foi possível construir uma interface que permite selecionar o conjunto de dados, pré-ajustar a curva aos dados por meio de deslizadores e ajustar o modelo de regressão não linear. O código abaixo fornece as ferramentas para obter o que se vê nos gifs animados dessa matéria. Modificações devem ser feitas para outros dados/modelos. Minha intenção no futuro é aprimorar essas funções de forma que o usuário forneça o modelo e os limites dos intervalos para os parâmetros diretamente pela interface. Em caso de sucesso, isso pode virar um pacote ou complementar algum dos pacotes existentes dedicados à regressão não linear.

O primeiro exemplo usa gráficos do pacote graphics, básico do R. Os dados são de 10 curvas de retenção de água no solo onde ajustou-se o modelo van Genuchten (1980). O segundo conjunto de dados é de um experimento para avaliar o impacto da desfolha na produtividade do algodão. Para esse caso adotou-se o modelo potência e gráficos do pacote lattice. Clique na caixa abaixo para ver o código. Até à próxima ridícula.

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

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

#-----------------------------------------------------------------------------
# dados de curva de retenção de água no solo

cra <- read.table("http://www.leg.ufpr.br/~walmes/data/cra_manejo.txt",
                  header=TRUE, sep="\t")
str(cra)
cra$tens[cra$tens==0] <- 0.1

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

cras <- subset(cra, condi=="LVA3,5")
cras <- with(cras, aggregate(cbind(umid),
                             list(posi=posi, tens=tens, prof=prof), mean))
cras$caso <- with(cras, interaction(posi, prof))

xyplot(umid~log10(tens)|posi, groups=prof, data=cras, type=c("p","a"))

#-----------------------------------------------------------------------------
# ajuste do modelo van Genuchten de modo iterativo com rpanel

# nlsajust é excutada quando clicar no botão "Ajustar"
nlsajust <- function(panel){
  ## ajuste do modelo não linear
  n0 <- try(nls(umid~tr+(ts-tr)/(1+(a*tens)^n)^(1-1/n), # modelo
                data=da,                                # dados
                start=start))                           # valores iniciais
  ## em caso de não convergência imprime mensagem de erro
  if(class(n0)=="try-error"){
    par(usr=c(0, 1, 0, 1))
    text(0.5, 0.5, "Não convergiu!\nAproxime mais.", col="red", cex=2)
  } else {
    ## coloca à curva ajusta sobre os pontos
    with(as.list(coef(n0)), curve(tr+(ts-tr)/(1+(a*10^x)^n)^(1-1/n),
                                  add=TRUE, col=2))
    ## salva o ajuste numa lista
    aju[[i]] <<- n0
  }
  panel
}

vg <- function(panel){
  ## seleciona o subconjunto dos dados
  i <<- panel$caso
  da <<- subset(cras, caso==i)
  ## lista com valores iniciais vindos dos deslizadores
  start <<- panel[c("tr","ts","a","n")]
  ## diagrama de dispersão
  plot(umid~log10(tens), data=da,
       ylab=expression(Conteúdo~de~água~(g~g^{-1})),
       xlab=expression(log[10]~Tensão~matricial~(kPa)),
       main=expression(f(x)==tr+(ts-tr)/(1+(a*x)^n)^{1-1/n}))
  ## sobrepõe a curva controlada pelos deslizadores
  with(start, curve(tr+(ts-tr)/(1+(a*10^x)^n)^(1-1/n),
                    add=TRUE, col=2, lty=2))
  panel
}

par(mar=c(4.1,4.2,3.1,1))
# cria objetos vazios que serão preenchidos durante processo
da <- c(); start <- list(); aju <- list(); i <- c()
# abre interface gráfica
panel <- rp.control()
# controla parâmetros do modelo
rp.slider(panel, ts, 0.4, 0.8, initval=0.6, showvalue=TRUE, action=vg)
rp.slider(panel, tr, 0, 0.5, initval=0.3, showvalue=TRUE, action=vg)
rp.slider(panel, a, 0.1, 5, initval=1.3, showvalue=TRUE, action=vg)
rp.slider(panel, n, 1, 5, initval=1.6, showvalue=TRUE, action=vg)
# seleciona o conjunto de dados
rp.listbox(panel, caso, vals=levels(cras$caso),
           title="Subconjunto", action=vg)
# cria botão "Ajustar"
rp.button(panel, action=nlsajust, title="Ajustar")

sapply(aju, coef)    # estimativas dos coeficientes
lapply(aju, summary) # summary dos modelos ajustados

#-----------------------------------------------------------------------------
# dados de peso de capulhos em função da desfolha do algodão

cap <- read.table("http://www.leg.ufpr.br/~walmes/data/algodão.txt",
                  header=TRUE, sep="\t", encoding="latin1")
str(cap)
cap$desf <- cap$desf/100
cap <- subset(cap, select=c(estag, desf, pcapu))
cap$estag <- factor(cap$estag, labels=c("vegetativo","botão floral",
                                 "florescimento","maçã","capulho"))
str(cap)

xyplot(pcapu~desf|estag, data=cap, layout=c(5,1),
       xlab="Nível de desfolha artificial", ylab="Peso de capulhos")

#-----------------------------------------------------------------------------
# o exemplo a seguir mostra como usar com gráficos do pacote lattice

# função adapatada para outro modelo de regressão não linear
nlsajust <- function(panel){
  n0 <- try(nls(pcapu~f0-f1*desf^exp(C), data=da, start=start))
  if(class(n0)=="try-error"){
    trellis.focus("panel", nivel, 1, highlight=FALSE)
    grid.text(x=0.5, y=0.5, label="Não convergiu!\nAproxime mais.",
              gp=gpar(col="red"))
    trellis.unfocus()
  }
  trellis.focus("panel", nivel, 1, highlight=FALSE)
  with(as.list(coef(n0)), panel.curve(f0-f1*x^exp(C), add=TRUE, col=2))
  trellis.unfocus()
  print(coef(n0))
  aju[[i]] <<- n0
  panel
}

# função adapatada para outro modelo de regressão não linear
ptn <- function(panel){
  nivel <<- as.numeric(panel$nivel)
  i <<- levels(cap$estag)[nivel]
  da <<- subset(cap, estag==i)
  start <<- panel[c("f0","f1","C")]
  print({
    xyplot(pcapu~desf|estag, data=cap, layout=c(5,1), col=1,
           xlab="Nível de desfolha artificial", ylab="Peso de capulhos",
           main=expression(f(x)==f[0]-f[1]*x^exp(C)),
           strip=strip.custom(bg="gray90"))
  })
    trellis.focus("panel", nivel, 1, highlight=FALSE)
    with(start, panel.curve(f0-f1*x^exp(C), add=TRUE, col=2, lty=2))
    trellis.unfocus()
  panel
}

da <- c(); start <- list(); aju <- list(); i <- c(); nivel <- c()
panel <- rp.control()
rp.slider(panel, f0, 10, 50, initval=30, showvalue=TRUE, action=ptn)
rp.slider(panel, f1, 0, 30, initval=10, showvalue=TRUE, action=ptn)
rp.slider(panel, C, -5, 5, initval=0, showvalue=TRUE, action=ptn)
rp.listbox(panel, nivel, vals=1:nlevels(cap$estag),
           title="Painel", action=ptn)
rp.button(panel, action=nlsajust, title="Ajustar")

sapply(aju, coef)
lapply(aju, summary)

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

Ajuste do modelo potência à dados de redução na produção de algodão devido à desfolha com interface gráfica desenvolvida com as ferramentas do pacote rpanel.

Anúncios

Fatorial duplo com dois tratamentos adicionais

Número de sementes em 40 infectadas por Fusarium (tranformado) em função da aplicação de fungicida em combinação com polímero.

Experimentos fatoriais são aqueles em que dois ou mais fatores estão presentes. Nos fatoriais incompletos nem todas as combinações possíveis estão presentes. Algumas das combinações podem não ser de interesse e por isso não são realizadas, alguma combinação é acidentalmente perdida durante o experimento ou certas combinações não identificam novos níveis, como por exemplo aplicação do produto A e B na dose 0 que, apesar dos rótulos diferentes, são o mesmo tratamento.

Um caso comum de fatorial incompleto são os fatoriais (completos) com a presença de tratamentos adicionais. Experimentos fatorais são realizados para verificar se os fatores interagem e então escolher combinação de níveis mais promissora para o processo estudado. Caso não exista interação, entende-se que os fatores têm efeito aditivos. Em qualquer caso, adicionam-se tratamentos adicionais ao experimento quando se deseja comparar o promissor aos de referência ou extremos.

Esse é um tipo de ensaio típico na avaliação do efeito de insumos agrícolas (fungicidas, insecticidas e herbicidas). Os ensaios testam produtos e suas doses e precisam ter as chamadas testemunhas para avaliar a eficiência das combinações geradas. Devido não ser um delineamento de níveis completamente cruzados a análise do experimento fatorial com tratamentos adicionais têm alguns detalhes que exigem mais atenção.

Nessa matéria faço análise de um fatorial com dois tratamentos adicionais. Os dados analisados foram retirados do Zimmermann (2004) de um experimento que estudou três fungicidas (benlate, captam, derosal) em três formas de aplicação (puro, antes de polímero, misturado com polímero). Adicionou-se dois tratamentos que foram a testemunha e aplicação do polímero puro. A resposta
avaliada foi o número de sementes infectadas com Fusarium em um total de 40 sementes.

O código está comentado para informar o leitor. São realizados a importação e reordenação dos níveis, a análise de variância com decomposição das somas de quadrados, contrastes entre efeitos, e testes sobre as médias ajustadas. Foi dada prioridade para os procedimentos matriciais de análise que podem ser adaptados para situações de desbalanceamento. Até a próxima ridícula.

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

require(lattice)
require(latticeExtra)
require(grid)
require(MASS)
require(doBy)
require(gmodels)
require(multcomp)

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

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

#-----------------------------------------------------------------------------
# modificar ordem dos níveis, primeiro níveis fatoriais, depois adicionais

levels(fu$fungicida)               # ordem correta
levels(fu$aplicacao)               # precisa alterar
levels(fu$aplicacao)[c(4,1,2,3,5)] # para ficar assim

fu$aplicacao <- factor(fu$aplicacao,
                       levels=levels(fu$aplicacao)[c(4,1,2,3,5)])
levels(fu$aplicacao) # agora ordem correta

#-----------------------------------------------------------------------------
# análise exploratória

xtabs(~fungicida+aplicacao, fu)

# os níveis não se cruzam completamente, portanto é um fatorial incompleto
# no caso 3x3+2, fatorial com 2 tratamentos adicionais

xyplot(infec40~tratamentos, fu) # resposta do tipo contagem de sucessos
xyplot(infec40~tratamentos, fu, jitter.x=TRUE)

#-----------------------------------------------------------------------------
# Vamos primeito analisar fora da estrutura fatorial, ou seja, considerando
# 3x3+2 = 11 níveis de um único fator.
# Vamos considerar distribuição para reproduzir os resultados do livro
# do Zimmermann (2004). Todavia, não é a análise mais correta devido os dados
# serem discretos. A distribuição binomial parece ser mais indicada.

m0 <- lm(infec40~tratamentos, data=fu)
par(mfrow=c(2,2)); plot(m0); layout(1) # não tão ruim, tentar boxcox?

m0 <- lm(infec40+1~tratamentos, data=fu) # soma 1 para eliminar os 0
boxcox(m0) # aponta log

m1 <- lm(log(infec40+1)~tratamentos, data=fu)
par(mfrow=c(2,2)); plot(m1); layout(1) # alguma melhoria?

# Zimmermann usa asin(sqrt(x/40)).
m2 <- lm(asin(sqrt(infec40/40))~tratamentos, data=fu)
par(mfrow=c(2,2)); plot(m2); layout(1) # alguma melhoria?

#-----------------------------------------------------------------------------
# Visualmente todos os resíduos são semelhantes, apresentam fugas.

lapply(list(m0, m1, m2), function(i) shapiro.test(residuals(i)))
lapply(list(m0, m1, m2), function(i) ks.test(scale(residuals(i)), "pnorm"))

# Pelos p-valores, a tranformação Box-Cox foi melhor, embora os resíduos ainda
# não apresentem um comportamento satísfatório.
# O foco não é testar transformações, então vamos seguir o que Zimmermann fez
# e trabalhar com asin(sqrt(x/40)).

#-----------------------------------------------------------------------------
# análise exploratória com a variável transformada

fu$y <- asin(sqrt(fu$infec40/40)) # cria a transformada

xyplot(y~tratamentos, fu, jitter.x=TRUE)

xyplot(y~fungicida, groups=aplicacao, fu, jitter.x=TRUE,
       type=c("p","a"), auto.key=TRUE)

# Apresenta indicativos de efeito de interação e que os adicionais
# são piores que o os fatoriais.

#-----------------------------------------------------------------------------
# análise na estrutura fatorial 3x3+2

m3 <- lm(y~fungicida*aplicacao, data=fu)
anova(m3)

# Os graus de liberdade são o que são pois:
# fungicida 4: são 5 niveis (3 fatoriais, 2 adicionais), 5-1=4;
# aplicacao 4: são 5, mas 2 estão confundidos (os adicionais), 5-2-1=2;
# fun:apli  4: (3-1)*(3-1)=4;

# O que queremos é a anova com somas de quadrados particionadas da seguinte
# forma
# fonte de variação     gl
# fungicidas             2
# aplicacao              2
# fun:apli               4
# fatorial vs adicionais 1
# polimero vs testemunha 1

#-----------------------------------------------------------------------------
# Na minha opinião, obter esse quadro de análise de variância é desperdício
# de tempo pois:
# * não tem sentido prático comparar o efeito médio dos fatoriais contra o
#   efeito médio dos adicionais, principalmente sabendo que a interação
#   fungicida:aplicacao foi significativa;
# * a interação fungicida:aplicacao significativa torna não interpretável
#   o teste para os efeitos principais;
# * polímero vs testemunha tem só um grau de liberdade, posso comprar os
#   efeitos por um testes t;

coef(m3) # tem NA para os efeitos não estimados (não presentes)
efstm <- colnames(vcov(m3)) # nomes dos efeitos que foram estimados
b <- coef(m3)[efstm] # só os estimados
b

# matriz de coeficientes para estimação das médias ajustadas
X <- popMatrix(m3, effect=c("fungicida", "aplicacao"))
str(X)

Xn <- attr(X, "grid") # 25 linhas, 5x5=25 (faz pensando que é completo)
Xo <- unique(fu[,c("fungicida","aplicacao")]) # que ocorrem (3x3+2)
Xo

#-----------------------------------------------------------------------------
# calculando as médias ajustadas matricialmente

rownames(X) <- apply(Xn, 1, paste, collapse=":")
on <- apply(Xo, 1, paste, collapse=":")
X <- X[on,efstm] # só linhas e colunas que correspondem à efeitos que ocorrem

X%*%b                                             # são as médias ajustadas
with(fu,
     tapply(y, list(fungicida, aplicacao), mean)) # são médias amostrais

# elas coincidem por ser um experimento regular (balanceado/ortogonal)

#-----------------------------------------------------------------------------
# fazendo contrastes
# não funcionam com presenta de NA
# contrast::contrast()
# multcomp::glht()
# para isso usar a gmodels::estimable()

m4 <- aov(formula(m3), data=fu) # aov não NA em coef

# polímero contra a testemunha
p.vs.t <- matrix(X[10,]-X[11,], nrow=1)
rownames(p.vs.t) <- "p.vs.t"
estimable(m4, cm=p.vs.t)

# fatorial contra adicional
f.vs.a <- matrix(apply(X[1:9,], 2, mean)-apply(X[10:11,], 2, mean), 1)
rownames(f.vs.a) <- "f.vs.a"
estimable(m4, cm=f.vs.a)

estimable(m4, cm=rbind(p.vs.t, f.vs.a)) # os dois de uma só vez

# E aí, complicado fazer tudo isso, né? Qual a sensação em saber de que
# pouco disso vai ser útil? Se deu interação teremos que estudar níveis
# de um fator fixando os níveis do outro.

#-----------------------------------------------------------------------------
# fazendo a decomposição das somas de quadrados dentro da anova

cf <- contrasts(C(fu$fungicida, helmert))
cf # o segredo está em montar os contrastes corretamente
cf[,3] <- c(-1,-1,-1,1.5,1.5) # fatorial contra adicional
cf[,4] <- c(0,0,0,-1,1)       # polímero vs testemunha
cf
contrasts(fu$fungicida) <- cf

m5 <- aov(y~fungicida*aplicacao, data=fu)
summary(m5, expand=FALSE,
        split=list("fungicida"=list(
                     "fungicida"=1:2,
                     "fat.vs.adi"=3,
                     "pol.vs.tes"=4)))

# Esse procedimento depende de balanceamento pois as somas de quadrados são
# sequênciais.

#-----------------------------------------------------------------------------
# Uma vez que deu interação é de prática fazer o desdobramento que consiste
# em estudar o efeito de aplicação dentro dos níveis de fungicida e
# vice-versa. Além do mais, por ter tratamentos adicionais, pode-se comparar
# com os níveis dos fatoriais. Vamos listas o número de hipóteses que teremos
# que fazer:
# * 3 níveis de aplicação em cada fungicida: 3*choose(3,2)=9
# * 3 níveis de fungicida em cada aplicação: 3*choose(3,2)=9
# * testemunha contra cada um dos 9 níveis do 3x3: 9
# * polimero contra cada um dos 9 níveis do 3x3: 9
# Com isso são 4*9=36 comparações. Já que é assim, vou fazer todos
# contra todos, choose(11, 2)=55, que o usuário pode olhar a que lhe
# interessa, afinal de contas, de 36 para 55 não têm diferença

mm <- model.matrix(m4)
dim(mm)
mm <- mm[,efstm]

# modelo de classe lm equivalente à m5, sem NA
m6 <- lm(y~-1+mm, data=fu)
anova(m6)

cbind(coef(m4), coef(m6)) # são o mesmo modelo

#-----------------------------------------------------------------------------
# médias ajustadas (ma)

cima <- confint(glht(m6, linfct=X)) # intervalos de confiança para ma
cima # tem cobertura *global* de 95% (e não indivídual)

str(cima)
cima$confint

result <- cbind(Xo, cima$confint)
str(result)

segplot(fungicida:aplicacao~lwr+upr, data=result,
        xlab="Plantas infectadas (asin(sqrt(x/40))",
        ylab="Fungicida:aplicação",
        centers=Estimate, draw.bands=FALSE,
        segments.fun=panel.arrows, ends="both",
        angle=90, length=1, unit="mm")

#-----------------------------------------------------------------------------
# carrega algumas funções que serão úteis para comparar médias e fazer
# gráficos com intervalos de confiança

objs <- ls()
objs <- c(objs, "apc", "prepanel.cbH", "panel.cbH")
source("http://www.leg.ufpr.br/~walmes/ridiculas/ridiculas_functionsi.R",
       encoding="latin1")

rm(list=ls()[!(ls()%in%objs)])
ls()

#-----------------------------------------------------------------------------
# comparando as médias ajustadas, choose(11, 2), contrastes de Tukey

str(X)
Xc <- apc(X)
str(Xc)

#-----------------------------------------------------------------------------
# fazendo as comparações com a glht, método para correção de p-valor
# fdr, o single-step demora muito

c0 <- summary(glht(m6, linfct=Xc), test=adjusted(type="fdr"))
c0

c0$focus <- "comparacoes" # para poder usar a cld()
cld(c0)

result$cld <- cld(c0)$mcletters$Letters

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

segplot(fungicida:aplicacao~lwr+upr, data=result,
        xlab="Plantas infectadas (asin(sqrt(x/40))",
        ylab="Fungicida:aplicação",
        centers=Estimate, draw.bands=FALSE,
        segments.fun=panel.arrows, ends="both",
        angle=90, length=1, unit="mm",
        panel=function(x, y, z, centers, subscripts, ...){
          panel.segplot(x, y, z, centers=centers, subscripts=subscripts, ...)
          panel.text(centers, z,
                     paste(format(centers, digits=2),
                           result$cld[subscripts], sep=" "), pos=3)
        })

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

result$desloc <- 0.05*c(0,-1,1,0,0)[match(result$aplicacao,
                                          levels(fu$aplicacao))]

#png("f048.png", 500, 400)
xyplot(y~fungicida, groups=aplicacao, data=fu, jitter.x=TRUE, amount=0.05,
       xlab="Fungicida", ylab=expression(asin(sqrt(x/40))),
       prepanel=prepanel.cbH, ly=result$lwr, uy=result$upr,
       auto.key=list(title="Aplicação", cex.title=1.2, columns=2,
         type="o", divide=1, lines=TRUE, points=FALSE))+
  as.layer(xyplot(Estimate~fungicida, groups=aplicacao, data=result, type="l",
                  ly=result$lwr, uy=result$upr,
                  desloc=result$desloc,
                  cty="bars",
                  panel.groups=panel.cbH,
                  panel=panel.superpose))
trellis.focus("panel", 1, 1, h=FALSE)
x <- as.numeric(result$fungicida)+result$desloc
y <- result$Estimate
lab <- paste(format(y, digits=2), result$cld, sep=" ")
a <- paste(rep("0", max(nchar(lab))), collapse="")
grid.rect(x=unit(x, "native"), y=unit(y, "native"),
          width=unit(1, data=a, "strwidth"),
          height=unit(1.2, data=a, "strheight"),
          hjust=-0.1, gp=gpar(fill="gray90", col=NA, fontsize=10))
grid.text(lab, x=unit(x, "native"), y=unit(y, "native"),
          hjust=-0.2, gp=gpar(col="black", fontsize=10))
trellis.unfocus()
#dev.off()

#-----------------------------------------------------------------------------
Categorias:delineamento Tags:, ,

Como fazer e interpretar o gráfico quantil-quantil

Gráfico quantil-quatil da distribuição normal com envelopes de confiança simulados.

O gráfico quantil-quantil (q-q) é uma ferramenta muito útil para checar adequação de distribuição de frequência dos dados à uma distribuição de probabilidades. Situações como essa ocorrem principalmente na análise de resíduos de modelos de regressão onde o gráfico q-q é usado para verificar se os resíduos apresentam distribuição normal. O gráfico q-q é melhor que o histograma e o gráfico de distribuição acumulada empírica porque nós temos mais habilidade para verificar se uma reta se ajusta aos pontos do que se uma curva de densidade se ajusta a um histograma ou uma curva de probabilidade acumulada se ajusta à acumulada empírica. Compare às três visualizações com o código a seguir.

#-----------------------------------------------------------------------------
# q-q vs histograma vs ecdf

y <- rnorm(50)
par(mfrow=c(1,3))
qqnorm(y); qqline(y)
plot(density(y))
curve(dnorm(x, mean(y), sd(y)), add=TRUE, col=2)
plot(ecdf(y))
curve(pnorm(x, mean(y), sd(y)), add=TRUE, col=2)

# qual você sente mais segurança para verificar adequação?
# nossos olhos têm mais habilidade para comparar alinhamentos
#-----------------------------------------------------------------------------

Apesar de muito usado, poucos usuários conhecem o procedimento para fazê-lo e interpretá-lo. O procedimento é simples e pode ser estendido para outras distribuições de probabilidade, não apenas para a distribuição normal como muitos podem pensar. Além do mais, alguns padrões do gráfico q-q obdecem à certas características dados dados, como assimetria, curtose e discreticidade. Saber identificar essas características é fundamental para indicar uma transformação aos dados. Abaixo o código para o gráfico q-q para distribuição normal, q-q para qualquer distribuição e o q-q com envelope de confiança obtido por simulação. A execução e estudo do código esclarece o procedimento. O envelope de confiança por simulação torna-se proibitivo para grandes amostras pelo tempo gasto na simulação.

#-----------------------------------------------------------------------------
# função para fazer o gráfico quantil-quantil da normal

qqn <- function(x, ref.line=TRUE){
  x <- na.omit(x)               # remove NA
  xo <- sort(x)                 # ordena a amostra
  n <- length(x)                # número de elementos
  i <- seq_along(x)             # índices posicionais
  pteo <- (i-0.5)/n             # probabilidades teóricas
  qteo <- qnorm(pteo)           # quantis teóricos sob a normal padrão
  plot(xo~qteo)                 # quantis observados ~ quantis teóricos
  if(ref.line){
    qrto <- quantile(x, c(1,3)/4) # 1º e 3º quartis observados
    qrtt <- qnorm(c(1,3)/4)       # 1º e 3º quartis teóricos
    points(qrtt, qrto, pch=3)     # quartis, passa uma reta de referência
    b <- diff(qrto)/diff(qrtt)    # coeficiente de inclinação da reta
    a <- b*(0-qrtt[1])+qrto[1]    # intercepto da reta
    abline(a=a, b=b)              # reta de referência
  }
}

x <- rnorm(20)
par(mfrow=c(1,2))
qqn(x)
qqnorm(x); qqline(x)
layout(1)

#-----------------------------------------------------------------------------
# função para fazer o gráfico quantil-quantil de qualquer distribuição

qqq <- function(x, ref.line=TRUE, distr=qnorm, param=list(mean=0, sd=1)){
  x <- na.omit(x)               # remove NA
  xo <- sort(x)                 # ordena a amostra
  n <- length(x)                # número de elementos
  i <- seq_along(x)             # índices posicionais
  pteo <- (i-0.5)/n             # probabilidades teóricas
  qteo <- do.call(distr,        # quantis teóricos sob a distribuição
                  c(list(p=pteo), param))
  plot(xo~qteo)                 # quantis observados ~ quantis teóricos
  if(ref.line){
    qrto <- quantile(x, c(1,3)/4) # 1º e 3º quartis observados
    qrtt <- do.call(distr,        # 1º e 3º quartis teóricos
                    c(list(p=c(1,3)/4), param))
    points(qrtt, qrto, pch=3)     # quartis, por eles passa uma reta de referência
    b <- diff(qrto)/diff(qrtt)    # coeficiente de inclinação da reta
    a <- b*(0-qrtt[1])+qrto[1]    # intercepto da reta
    abline(a=a, b=b)              # reta de referência
  }
}

x <- rnorm(20)
par(mfrow=c(1,2))
qqq(x)
qqnorm(x); qqline(x)
layout(1)

x <- runif(20)
qqq(x, ref.line=TRUE, distr=qunif, param=list(min=0, max=1))

x <- rgamma(20, shape=4, rate=1/2)
qqq(x, ref.line=TRUE, distr=qgamma, param=list(shape=4, rate=1/2))

#-----------------------------------------------------------------------------
# envelope para o gráfico de quantis (simulated bands)

qqqsb <- function(x, ref.line=TRUE, distr=qnorm, param=list(mean=0, sd=1),
                  sb=TRUE, nsim=500, alpha=0.95, ...){
  x <- na.omit(x)               # remove NA
  xo <- sort(x)                 # ordena a amostra
  n <- length(x)                # número de elementos
  i <- seq_along(x)             # índices posicionais
  pteo <- (i-0.5)/n             # probabilidades teóricas
  qteo <- do.call(distr,        # quantis teóricos sob a distribuição
                  c(list(p=pteo), param))
  plot(xo~qteo, ...)            # quantis observados ~ quantis teóricos
  if(ref.line){
    qrto <- quantile(x, c(1,3)/4) # 1º e 3º quartis observados
    qrtt <- do.call(distr,        # 1º e 3º quartis teóricos
                    c(list(p=c(1,3)/4), param))
    points(qrtt, qrto, pch=3)     # quartis, passa uma reta de referência
    b <- diff(qrto)/diff(qrtt)    # coeficiente de inclinação da reta
    a <- b*(0-qrtt[1])+qrto[1]    # intercepto da reta
    abline(a=a, b=b)              # reta de referência
  }
  if(sb){
    rdistr <- sub("q", "r",       # função que gera números aleatórios
                  deparse(substitute(distr)))
    aa <- replicate(nsim,         # amostra da distribuição de referência
                    sort(do.call(rdistr, c(list(n=n), param))))
    lim <- apply(aa, 1,           # limites das bandas 100*alpha%
                 quantile, probs=c((1-alpha)/2,(alpha+1)/2))
    matlines(qteo, t(lim),        # coloca as bandas do envelope simulado
             lty=2, col=1)
  }
}

x <- rnorm(20)

#png("f047.png", 400, 300)
qqqsb(x, xlab="Quantis teóricos da distribuição normal padrão",
      ylab="Quantis observados da amostra",
      main="Gráfico quantil-quantil da normal")
#dev.off()

x <- rpois(50, lambda=20)
qqqsb(x, distr=qpois, param=list(lambda=20))

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

A normalidade dos resíduos é um dos pressupostos da análise de regressão. São os resíduos e não os dados que devem apresentar normalidade. Se a distribuição dos dados, ou melhor, da sua variável resposta (Y) condicional ao efeito das suas variáveis explicativas for normal, os resíduos terão distribuição normal. Porém, se você aplica um teste de normalidade aos dados (Y) você não está considerando os efeitos das variáveis explicativas, ou seja, você está aplicando um teste na distribuição marginal de Y que não tem porque atender a normalidade. Todo teste de normalidade supõe que os dados têm uma média e uma variância e só os resíduos atendem essa premissa porque os dados (Y) têm média dependente do efeito das variáveis explicativas.

Distribuição condicional e marginal da variável resposta.

#-----------------------------------------------------------------------------
# distribuição condicional vs distribuição marginal

layout(1)
x <- rep(1:4, e=50)
y <- rnorm(x, mean=0+0.75*x, sd=0.1)
da <- data.frame(x, y)
plot(y~x, da)

db <- split(da, f=da$x)

#png("f046.png", 400, 300); par(mar=c(1,1,1,1))
plot(y~x, da, xlim=c(1,6), xaxt="n", yaxt="n", xlab="", ylab="")
abline(a=0, b=0.75)
lapply(db,
       function(d){
         dnst <- density(d$y)
         lines(d$x[1]+dnst$y*0.1, dnst$x)
         abline(v=d$x[1], lty=2)
       })
points(rep(5, length(y)), da$y, col=2)
abline(v=5, lty=2, col=2)
dnst <- density(da$y)
lines(5+dnst$y*2, dnst$x, col=2)
text(3, 1, "distribuição\ncondicional")
text(5.5, 1, "distribuição\nmarginal", col=2)
#dev.off()

par(mfrow=c(1,2))
qqnorm(da$y)
qqnorm(residuals(lm(y~x, da)))
layout(1)

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

Muitos usuários preferem aplicar um teste de normalidade do que olhar para o gráfico q-q. Isso têm duas razões: (1) costume, o usuário sempre usou aplicativos para análise de dados que não dispõem de recursos gráficos, eles conduzem toda análise sem sequer ver os dados em um gráfico, (2) consideram subjetiva a análise gráfica. Do meu ponto de vista, a subjetividade está presente ao aplicar um teste também pois o usuário é quem escolhe o teste e o nível de significância. Outro ponto é que os testes de normalidade assumem independência entre as observações e os resíduos não são independentes e foram obtidos após a estimação de p parâmetros, o que não é considerado por nenhum teste. Ou seja, qualquer teste aplicado fornece um resultado aproximado. Mas o que de fato eu defendo é que a análise gráfica é indiscutivelmente mais informativa. Veja, se o teste rejeita a normalidade é porque os dados não apresentam distribuição normal por algum motivo. Quando você visualiza o q-q é possível explicar a fuga de normalidade que pode ser sistemática: (1) devido à desvio de assimetria, (2) de curtose, (3) à mistura de distribuições, (4) à presença de um outlier e (5) ao dados serem discretos. Essas são as principais causas de afastamento. Cada uma delas sugere uma alternativa para corrigir a fuga: tranformação, modelagem da variância, deleção de outlier, etc. Nesse sentido, como identificar esses padrões de fuga? É o que os gráficos animados vão mostrar. Até a próxima ridícula.

Fuga de normalidade por assimetria representada pelo gráfico quantil-quantil. A assimetria no q-q aparece com pontos dispostos em forma de arco.

#-----------------------------------------------------------------------------
# assimetria

dir.create("frames")
setwd("frames")

png(file="assimet%03d.png", width=300, height=150)
par(mar=c(1,1,1,1))
par(mfrow=c(1,2))
for(i in 10*sin(seq(0.01, pi-0.01, l=100))){
  curve(dbeta(x, i, 10-i), 0, 1, xaxt="n", yaxt="n", xlab="", ylab="")
  y <- rbeta(100, i, 10-i)
  qqnorm(y, xaxt="n", yaxt="n", main=NULL, xlab="", ylab=""); qqline(y)
}
dev.off()

# converte os pngs para um gif usando ImageMagick
system("convert -delay 10 assimet*.png assimet.gif")

# remove os arquivos png
file.remove(list.files(pattern=".png"))

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

Fuga de normalidade por mistura de distribuições (alteração de curtose) representada pelo gráfico quantil-quantil. Mistura de distribuições geram fugas nos extremos.

#-----------------------------------------------------------------------------
# mistura de distribuições

png(file="mistura%03d.png", width=300, height=150)
par(mar=c(1,1,1,1))
par(mfrow=c(1,2))
for(i in sin(seq(0, pi, l=100))){
  curve(i*dnorm(x,0,1)+(1-i)*dnorm(x,0,6), -20, 20,
        xaxt="n", yaxt="n", xlab="", ylab="")
  y <- c(rnorm(ceiling(500*i),0,1), rnorm(floor(500*(1-i)),0,6))
  qqnorm(y, xaxt="n", yaxt="n", main=NULL, xlab="", ylab=""); qqline(y)
}
dev.off()

# converte os pngs para um gif usando ImageMagick
system("convert -delay 10 mistura*.png mistura.gif")

# remove os arquivos png
file.remove(list.files(pattern=".png"))

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

Fuga de normalidade devido dados serem discretos (padrão escada) representada pelo gráfico quantil-quantil.

#-----------------------------------------------------------------------------
# discreticidade

png(file="discret%03d.png", width=300, height=150)
par(mar=c(1,1,1,1))
par(mfrow=c(1,2))
for(i in c(1:100, 99:1)){
  x <- 0:i
  px <- dbinom(x, i, 0.5)
  plot(px~x, type="h", xaxt="n", yaxt="n", xlab="", ylab="")
  y <- rbinom(100, i, 0.5)
  qqnorm(y, xaxt="n", yaxt="n", main=NULL, xlab="", ylab=""); qqline(y)
}
dev.off()

# converte os pngs para um gif usando ImageMagick
system("convert -delay 10 discret*.png discret.gif")

# remove os arquivos png
file.remove(list.files(pattern=".png"))

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

Interpretação da matriz de covariância das estimativas

Quero falar de um erro comum de interpretação de resultados em análise de regressão através de um exemplo. Considere que você tem pessoas com diferentes pesos e diferentes alturas. Facilmente você aceita que uma pessoa mais alta tem maior peso, certo? Ou seja, existe uma correlação positiva entre peso e altura. Pois bem, vamos simular observações desse experimento e em seguida ajustar uma regressão linear simples para relação peso-altura.

require(MASS)

set.seed(12345)
da <- as.data.frame(mvrnorm(80, c(altura=30, peso=60),
                            matrix(c(2,1.9,1.9,3), 2,2)))
plot(peso~altura, da)

m0 <- lm(peso~altura, data=da) # ajusta a reta
summary(m0)              # estimativa dos parâmetros
abline(m0)               # adiciona uma reta ao gráfico
vcov(m0)                 # covariância das estimativas

Aqui está o ponto que eu quero comentar: a interpretação da matriz de covariância. Perceba que a covariância foi negativa. Já vi gente interpretando isso da seguinte forma: espera-se que pessoas mais altas tenham menor peso. Duplamente errado. Primeiro nos sabemos por experiência que a correlação de peso e altura é positiva. Segundo, a matriz de covariância se refere as estimativas dos parâmetros e não as variáveis envolvidas. Nunca a interprete dessa forma.

Então como interpretar? Bem, a matriz de covariância das estimativas, é um reflexo da função objetivo ao redor da solução. A função objetivo nesse caso é minimizar a soma de quadrados (SQ). Então se eu aumento o valor do parâmetro b0, o parâmetro b1 tem que diminuir para compensar, ou seja, para diminuir a SQ. Ocorre um efeito compensatório aqui. Além do mais, esse efeito pode ser eliminado facilmente aplicando uma reparametrização, como por exemplo, centrar a covariável na média.

da$alturac <- with(da, altura-mean(altura)) # centrando a covariável
par(mfrow=c(1,2))
plot(peso~altura, da)
abline(m0)
plot(peso~alturac, da)

m1 <- lm(peso~alturac, data=da) # ajusta a reta
summary(m1)              # estimativa dos parâmetros
abline(m1)               # adiciona uma reta ao gráfico
vcov(m1)                 # covariância das estimativas

require(ellipse)
par(mfrow=c(1,2)) # regiões de confiança conjunta
plot(ellipse(vcov(m0)), type="l")
plot(ellipse(vcov(m1)), type="l")

cov2cor(vcov(m0))
cov2cor(vcov(m1))

Perceba que os resultados são os mesmos em termos de estatísticas de teste e medidas de ajuste. Isso é esperado por eu só fiz uma traslação da altura. Mas o importante é que agora a matriz de covariância tem covariâncias praticamente nulas, que é resultado da translação. O intercepto é então o valor esperado para alturac igual a zero (centro dos dados). Veja porque essa covariância é zero: se eu aumentar b0, não adiante eu alterar b1 que não vai diminuir a SQ, e vice-versa, porque agora eles são ortogonais. Ortogonalidade entre parâmetros é uma característica desejável pois permite você inferir sobre um deles sem considerar o outro. Além disso, tem vantagens do ponto de vista de estimação por métodos numéricos de otimização. Em outras palavras, pegando conceitos de probabilidade, se a distribuição amostral de duas variáveis aleatórias é normal bivariada com covariância nula, a distribuição condicional de A|B é igual a distribuição marginal de A pela independência. Reforçando, eu não preciso conhecer valores de B para descrever A. Considerando tudo que foi argumentado, é sempre preferível que você adote o modelo que apresente menor covariância entre os parâmetros. Até a próxima ridícula.

Categorias:inferência Tags:, ,

Semântica para descrever modelos

A função lm() do R resolve 90% dos casos de análise de experimentos. Não são para lm() modelos não lineares, de efeito aleatório, de distribuição diferente da gaussiana ou dados não independentes (correlacionados). Para uso eficiente da função é importante saber declarar os modelos. Para isso existe uma semântica (termo bem apropriado usado pelo Fernando Toledo) para declarar os modelos. Aqui eu vou listar alguns dos casos possíveis.

fórmula significado
A+B efeito simples aditivo entre fatores
A:B efeito do produto cartesiano de níveis dos fatores
A*B*C efeitos simples e todas interações possíveis
(A+B+C)^2 efeitos simples e interações de até 2 termos (duplas)
(A+B+C+C+E)^3 efeitos simples e interações de até 3 termos (triplas)
A/B efeito simples de A e de B aninhado em A, A+A:B
I(x^2) operador I() torna x^2 literal (remove o poder semântico de ^2)
-A:B remove o termo do modelo
A*B*C-A:B:C declara modelo triplo mas remove interação tripla, (A+B+C)^2
0+A remove intercepto (coluna de 1) do modelo
-1+A idem ao anterior
poly(x, degree=2) polinômio ortogonal de grau 2 em x
poly(x, degree=2, raw=TRUE) polinômio cru de grau 2 em x
A*(B+C) propriedade distributiva, A*B+A*C
A/(B+C) idem ao anterior, A/B+A/C
A/B/C propriedade recursiva, A+A:B+A:B:C
. representa todas as colunas, exceto a da resposta
(.)^2 inclui todas as interações duplas das variáveis

Além dessas opções de sintaxe, ainda temos a possibilidade de usar o operador semântico Error() para declarar termos de efeito aleatório. Esse operador só é interpretado pela função aov(). Com ele fazemos a análise de experimentos em parcelas subdivididas, por exemplo, de tal forma que a análise de variância apresenta os testes F considerando os quadrados médios corretos.
Além disso podermos usar a função terms(..., keep.order=TRUE) para fixar a ordem dos termos na análise de variância, pois os termos são ordenados pela ordem de magnitude, assim, sempre interações vem depois de efeitos simples, algumas vezes não desejamos isso, como é o caso da análise de experimentos em blocos incompletos. Os exemplos vão tornar esses apontamentos mais claros. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# observe as matrizes de delineamento geradas nessas situações

da <- replicate(5, gl(3,1), simplify=FALSE)
da <- data.frame(do.call(expand.grid, da))
names(da) <- LETTERS[1:ncol(da)]
da$x <- runif(nrow(da))

formulas <- list(~A+B, ~A:B, ~A*B, ~(A+B)^2, ~(A+B+C+D+E)^3,
                 ~A/B, ~A/B/C, ~A/(B+C), ~A*(B+C), ~A*(B+C)^2,
                 ~(A+B+C)^2-A:B, ~x, ~0+A, ~-1+A, ~0+x, ~A*x, ~A:x,
                 ~A/x, ~x+I(x^2), ~poly(x, degree=2),
                 ~poly(x, degree=2, raw=TRUE), ~A+A:B+C,
                 terms(~A+A:B+C, keep.order=TRUE))
names(formulas) <- apply(sapply(formulas, paste), 2, paste, collapse="")

X <- lapply(formulas, model.matrix, data=da)
lapply(X, head)

#-----------------------------------------------------------------------------
# para obter a equação de um polinômio para os níveis

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

soja$A <- factor(soja$agua)
soja$K <- soja$potassio

m0 <- lm(rengrao~bloco+A*(K+I(K^2)), soja) # I() torna literal
m1 <- lm(rengrao~bloco+A*poly(K, degree=2, raw=TRUE), soja)  # original
m2 <- lm(rengrao~bloco+A*poly(K, degree=2, raw=FALSE), soja) # ortogonal

cbind(coef(m0), coef(m1), coef(m2)) # 2 primeiras colunas são iguais

m3 <- lm(rengrao~0+A/(K+I(K^2))+bloco, soja) # para ter os coeficientes
coef(m3) # os interceptos, termos lineares e termos quadráticos
c0 <- coef(m3)[grep("bloco", names(coef(m3)), invert=TRUE)]
matrix(c0, ncol=3, dimnames=list(levels(soja$A), c("b0","b1","b2")))

#-----------------------------------------------------------------------------
# tranformações dentro da fórmula

m0 <- lm(rengrao~bloco+A*(K+sqrt(K)), soja)  # não requer I()
m0 <- lm(rengrao~bloco+A*I(K-mean(K)), soja) # centrar na média

#-----------------------------------------------------------------------------
# fazendo de conta que o experimento é em parcelas subdividas para ilustrar
# o uso do operador semântico Error()

soja$K <- factor(soja$K)

# combinar os níveis de bloco e A nos identificamos as parcelas
m4 <- aov(rengrao~bloco+A*K+Error(bloco:A), soja) # dá mensagem de aviso
summary(m4) # anova separa pelos estratos

soja <- transform(soja, parcela=interaction(bloco, A))

# é o mesmo modelo
m5 <- aov(rengrao~bloco+A*K+Error(parcela), soja) # não dá mensagem de aviso
summary(m5) # anova separa pelos estratos

#-----------------------------------------------------------------------------
# fazendo uso da função terms() para fixar a ordem dos termos na anova()

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

bib1 <- transform(bib1,
                  repetição=factor(repetição),
                  variedade=factor(variedade),
                  bloco=factor(bloco))
str(bib1)

# repetição + bloco dentro de repetição + variedade
m0 <- lm(y~repetição/bloco+variedade, bib1)
anova(m0) # por ser de maior ordem o termo repetição:bloco é o último

m1 <- lm(terms(y~repetição/bloco+variedade, keep.order=TRUE), bib1)
anova(m1) # ordem obedeceu a ordem dentro da fórmula

#-----------------------------------------------------------------------------
Categorias:delineamento Tags:, ,

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:,