Arquivo

Archive for the ‘gráficos’ Category

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.

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

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

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

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

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

Média e desvio-padrão de muitas variáveis separado por grupos

Gráfico de barras com as médias e desvios-padrões para variáveis resposta de um experimento avaliando cultivares de banana e meios de cultura em cultura de tecidos.

Em estatística descritiva é comum a síntese dos dados em medidas de posição e dispersão. As medidas mais frequentemente usadas são a média e o desvio-padrão. Não raramente, existe a necessidade de obter essas medidas separando as observações por alguma variável categórica do conjunto de dados. Nesse post apresento algumas formas de obter isso. Estou certo de que as soluções que apresento abaixo não cobrem 5% das disponíveis considerando a extensa quantidade de pacotes/funções do R ligadas à análise descritiva, portanto, considere isso apenas como um ponto de partida. Por fim, sintetizo em um gráfico os resultamos aplicados em um conjunto de dados reais. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# gerando um conjunto de dados artificial

da <- expand.grid(A=1:3, B=1:2, r=1:15)
da <- cbind(da, as.data.frame(matrix(rnorm(nrow(da)*3), nrow(da))))
str(da)

#-----------------------------------------------------------------------------
# usando funções do reshape e plyr

require(reshape)
require(plyr)

ddply(da[,-c(1:3)], .(da$A, da$B), mean) # as médias
ddply(da[,-c(1:3)], .(da$A, da$B), sd)   # os desvios-padrões

db <- melt(da, id.vars=1:3, id.measures=4:6)
str(db)

#-----------------------------------------------------------------------------
# fazendo empilhamento das observações

ddply(db[,-3], .(A, B, variable), summarise, m=mean(value), sd=sd(value))

#-----------------------------------------------------------------------------
# usando funções do pacote doBy

library(doBy)

summaryBy(V1+V2+V3~A+B, data=da,
          FUN=function(x){ c(m=mean(x), sd=sd(x)) })

#-----------------------------------------------------------------------------
# usando funções do Hmisc

require(Hmisc)

with(da, bystats(cbind(V1, V2, V3), A, B,
                 fun=function(x){ c(m=colMeans(x), s=sd(x)) }))

#-----------------------------------------------------------------------------
# dados reais

bnn <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/cnpaf/banana.txt",
                  header=TRUE, sep="\t")
str(bnn)
bnn <- melt(bnn, id.measures=4:12)
bnn <- ddply(bnn, .(cult, meio, variable),
             summarise, m=mean(value), sd=sd(value))
str(bnn)

#-----------------------------------------------------------------------------
# gráfico de barras com desvio-padrão

require(lattice)
require(RColorBrewer)
display.brewer.all()
cols <- brewer.pal(nlevels(bnn$meio), "Set3")

png("f028.png", 500, 500)
trellis.par.set(superpose.polygon=list(col=cols),
                strip.background=list(col="gray90"))
barchart(m~cult|variable, groups=meio, data=bnn,
         sdv=bnn$sd, scales="free", layout=c(3,3),
         xlab="Cultivar", ylab=expression(bar(x)+sd),
         key=list(title="Meio de cultura", cex.title=1.2,
           columns=3, text=list(levels(bnn$meio)),
           rectangle=list(col=trellis.par.get("superpose.polygon")$col)),
         prepanel=function(y, sdv, subscripts=subscripts, ...){
           ly <- as.numeric(y+sdv[subscripts])
           list(ylim=range(y, ly, finite=TRUE))
         },
         panel=function(x, y, subscripts, groups, sdv, box.ratio, ...){
           panel.barchart(x, y, subscripts=subscripts,
                          groups=groups, box.ratio=box.ratio, ...)
           d <- 1/(nlevels(groups)+nlevels(groups)/box.ratio)
           g <- (as.numeric(groups[subscripts])-1); g <- (g-median(g))*d
           panel.arrows(as.numeric(x)+g, y, as.numeric(x)+g, y+sdv[subscripts],
                        angle=90, length=0.025)
         })
dev.off()

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

Moda de uma amostra: método do histograma vs método kernel

Moda de uma amostra aleatória de uma distribuição contínua obtida pelo método do histograma e da densidade kernel.

Uma pergunta bem comum nos fóruns de estatística, geralmente associado à alunos de graduação, é de alguma função no R para obter a moda de um conjunto de dados. Bem, como moda se refere ao valor de maior frequência na amostra, a tarefa é simples para dados discretos: basta obter a distribuição de frequência dos valores (função table()) e tomar aquele associado à maior frequência.

Para dados contínuos, teoricamente, não teríamos dois valores observados iguais em uma amostra. Em amostras reais às vezes temos, seja pelo fato da amostra ser grande e/ou pelo instrumento de medida não ser preciso para diferenciar duas observações. Por exemplo, uma balança que registra massa em quilos vai estar sempre arredondando (censurando) a parte decimal dos valores. Então na prática até que podemos ter, eventualmente, valores repetidos em amostras de variáveis contínuas. Porém, não é assim que calculamos moda nessas amostras, porque valores repetidos podem não ocorrer e quando ocorrem, raramente excede à duas ocorrência de um mesmo valor.

A obtenção da moda pelo método do histograma é simples e está bem documentada nos livros texto da disciplina de estatística básica. Consiste em agrupar os dados em classes e fazer uma interpolação linear que considera a frequência da classe modal e de suas vizinhas. A interpolação linear usa semelhança de triângulos para encontrar a moda da amostra. A fórmula é essa

mo= l_i+\displaystyle A_c \left(\frac{f_m-f_e}{2f_m-f_e-f_d}\right )

em que l_i é o limite inferior da classe modal, A_c é a amplitude da classe modal, f_m é a frequência da classe modal, f_e é a frequência da classe à esquerda da modal, f_d é a frequência da classe à direita da modal.

Outra forma de obter a classe modal é calculando o valor da amostra que corresponde ao máximo da função de densidade kernel. É um método mais moderno e computacionalmente mais trabalhoso, pois deve-se ajustar a densidade kernel, fazer aproximação da função kernel por uma função interpoladora e depois obter o máximo dessa função. Mas como somos usuários de R, já temos há um bom tempo as ferramentas necessárias para isso.

No CMR abaixo eu implementei os métodos para obtenção da moda de uma amostra pelos dois métodos: do histograma e da densidade kernel. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# amostra aleatória de distribuição discreta

x <- rpois(100, lambda=12)
plot(sort(x))
tb <- table(x); tb
tb[which.max(tb)]

#-----------------------------------------------------------------------------
# amostra aleatória de distribuição contínua

x <- rbeta(1000, 0.5,0.5)
ht <- hist(x) # histograma

#-----------------------------------------------------------------------------
# função para obter a moda pelo método do histograma, moda obtida por
# interpolação linear, baseado em semelhaça de triângulos

moda.hist <- function(ht, plotit=TRUE){
  ## ht: um objeto do uso da função hist()
  mcl <- which.max(ht$counts)
  li <- ht$breaks[mcl]
  width <- diff(ht$breaks[mcl+0:1])
  counts <- c(0,ht$counts,0)
  delta <- abs(diff(counts[1+mcl+(-1:1)]))
  moda <- li+width*delta[1]/sum(delta)
  cols <- rep(5, length(ht$counts))
  cols[mcl] <- 3
  if(plotit==TRUE){
    plot(ht, col=cols)
    abline(v=moda)
  }
  return(moda=moda)
}

ht <- hist(x)
moda.hist(ht)

#-----------------------------------------------------------------------------
# função para obter a moda da função densidade kernel

moda.dens <- function(dn, plotit=TRUE){
  ## dn: um objeto do uso da função density()
  ini <- dn$x[which.max(dn$y)]
  fx <- approxfun(dn$x, dn$y)
  op <- optim(c(ini), fx, method="BFGS", control=list(fnscale=-1))
  if(plotit==TRUE){
    plot(dn)
    abline(v=op$par, col=2)
  }
  return(moda=op$par)
}

dn <- density(x, kernel="triangular", bw=0.01)
moda.dens(dn)

dn <- density(x, kernel="gaussian", bw=0.05)
moda.dens(dn)

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

png("f027.png", 400, 300)
x <- rgamma(1000, 5, 5)
ht <- hist(x, freq=FALSE, col="gray86")
mo1 <- moda.hist(ht, plotit=FALSE)
abline(v=mo1, col=2, lwd=2)
dn <- density(x)
lines(dn, col=4, lwd=2)
mo2 <- moda.dens(dn, plotit=FALSE)
abline(v=mo2, col=3, lwd=2)
box()
dev.off()

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