Ajuste de modelos lineares e mistos no ambiente R

09 e 10 de Outubro de 2014 - Piracicaba - SP
Prof. Dr. Walmes M. Zeviani
Escola Superior de Agricultura “Luiz de Queiroz” - USP
Lab. de Estatística e Geoinformação - LEG
Pós Graduação em Genética e Melhoramento de Plantas Departamento de Estatística - UFPR

Análise exploratória

##-----------------------------------------------------------------------------
## Dados de carros Duster à venda no webmotors em 26/03/2014.

dus <-
    read.table("http://www.leg.ufpr.br/~walmes/data/duster_venda_260314.txt",
               header=TRUE, sep="\t", encoding="utf-8")
dus$ano <- factor(gsub(x=as.character(dus$ano), "/\\d{4}$", ""))
str(dus)

## Quantidade de NA em cada coluna.
apply(dus, MARGIN=2, function(x) sum(is.na(x)))

## Elimina registros com NA.
dus <- na.omit(dus)

##-----------------------------------------------------------------------------
## Gráfico de barras e setores.

x <- table(dus$cambio)
class(x)

## Se vem da xtabs() também tem classe `table`.
x <- xtabs(~cambio, data=dus)
class(x)

## barplot(x)
barplot(x,
        xlab="Tipo de câmbio",
        ylab="Frequência absoluta",
        col=c("seagreen", "yellowgreen"))

plot of chunk unnamed-chunk-2

barplot(x, horiz=TRUE,
        xlab="Tipo de câmbio",
        ylab="Frequência absoluta",
        col=c("seagreen", "yellowgreen"))
box(bty="L")

plot of chunk unnamed-chunk-2

## Cores com `green` no nome.
colors()
grep("green", colors(), value=TRUE)

## Gráfico de setores.
pie(x, col=c("seagreen", "yellowgreen"),
    main="Tipo de câmbio")

plot of chunk unnamed-chunk-2

## Para as cores do carro.
x <- xtabs(~cor, data=dus)
levels(dus$cor)

par(mar=c(4.1,7.1,2.1,2.1))
barplot(x, horiz=TRUE, las=1,
        col=c("blue", "white", "gray50", "Yellow", "gray90", "black",
            "green4", "red", "red4"))
mtext(side=2, text="Cor", line=5)
mtext(side=1, text="Frequência absoluta", line=2)
box(bty="L")

plot of chunk unnamed-chunk-2

##-----------------------------------------------------------------------------
## Gráficos de barras emplilhadas (stacked) e lado a lado.

x <- xtabs(~ano+cambio, data=dus)
x

## Barras empilhadas.
barplot(x, xlab="Câmbio", ylab="Frequência absoluta")

plot of chunk unnamed-chunk-2

colcamb <- c("seagreen", "yellowgreen")
barplot(t(x),
        xlab="Ano",
        ylab="Frequência absoluta",
        col=colcamb)
legend("topleft", legend=levels(dus$cambio),
       fill=colcamb, bty="n")

plot of chunk unnamed-chunk-2

## Barras lado a lado.
barplot(t(x), beside=TRUE,
        xlab="Ano", ylab="Frequência absoluta",
        col=colcamb)
legend("topleft", legend=levels(dus$cambio),
       fill=colcamb, bty="n")

plot of chunk unnamed-chunk-2

##-----------------------------------------------------------------------------
## Anotações nas barras.

x <- xtabs(~cambio+poten, data=dus); x

## Cores de preenchimento para as barras.
cols <- c("seagreen", "yellowgreen")

## Barras lado a lado.
bp <- barplot(t(x), beside=TRUE, col=cols,
              xlab="Tipo de câmbio", ylab="Frequência absoluta")

plot of chunk unnamed-chunk-2

bp

## Calcula a altura de uma palavra em termos da escala y do gráfico.
sh <- strheight("um texto qualquer"); sh
lim <- par()$usr[4]+3*sh

## Refaz o gráfico com espaço para o texto.
barplot(t(x), beside=TRUE, col=cols, ylim=c(0, lim),
              xlab="Tipo de câmbio", ylab="Frequência absoluta")
legend("topleft", title="Potência",
       legend=c("1.6","2.0"), fill=cols, bty="n")
text(x=c(bp), y=t(x), labels=x, pos=3)
box()

plot of chunk unnamed-chunk-2

##-----------------------------------------------------------------------------
## Gráficos de mosaico.

x <- xtabs(~ano+cambio, data=dus)
x

mosaicplot(x, ylab="Tipo de câmbio", xlab="Ano")

plot of chunk unnamed-chunk-2

mosaicplot(t(x), xlab="Tipo de câmbio", ylab="Ano")

plot of chunk unnamed-chunk-2

x <- xtabs(~novo+poten, data=dus); x

## Não dependência entre as variáveis.
mosaicplot(x, xlab="Condição", ylab="Potência",
           col=c("#009054","#900039"))

plot of chunk unnamed-chunk-2

## Pode-se especificar cores com a trinca RGB (red, green, blue),
## pode-se usar o padrão hexadecimal html para cores.

## Visite estes sites para pegar cores.
## browseURL("http://www.w3schools.com/html/html_colors.asp")
## browseURL("http://html-color-codes.info/")

##-----------------------------------------------------------------------------
## Histograma.

hist(dus$valor)

plot of chunk unnamed-chunk-2

hist(dus$valor, xlab="Preço de venda (R$)",
     ylab="Frequência absoluta", col="orange")
rug(dus$valor)

plot of chunk unnamed-chunk-2

## Se breaks é um escalar então entende-se que é uma *sugestão* para o
## número de clases.
hist(dus$valor,
     breaks=15,
     xlab="Preço de venda (R$)",
     ylab="Frequência absoluta",
     col="orange")
rug(dus$valor)

plot of chunk unnamed-chunk-2

## Se breaks é um vetor então entende-se que são os limites para
## classificação dos valores.
hist(dus$valor, breaks=seq(35000, 75000, 2500),
     xlab="Preço de venda (R$)",
     ylab="Frequência absoluta", col="#7700B7",
     sub="Amplitude de classe de R$ 2500", main=NULL)

plot of chunk unnamed-chunk-2

## Gráfico onde a altura é a densidade e não a frequência.
hist(dus$valor, prob=TRUE, breaks=seq(35000, 75000, 2500),
     xlab="Preço de venda (R$)",
     ylab="Frequência absoluta", col="#7700B7",
     sub="Amplitude de classe de R$ 2500", main=NULL)

plot of chunk unnamed-chunk-2

## Esse gráfico tem que a soma da área dos retângulos somam 1 pois o
## produto da amplitude pela densidade é a frequência relativa e a soma
## das frequência relativas é 1.

hist(dus$valor, prob=TRUE, seq(35000, 75000, 2000),
     xlab="Preço de venda (R$)", 
     ylab="Frequência absoluta", col="#6E0039",
     sub="Amplitude de classe de R$ 2500", main=NULL)
rug(dus$valor) ## Faz risquinhos no eixo x.

plot of chunk unnamed-chunk-2

##-----------------------------------------------------------------------------
## Anotações sobre um histograma.

## Com domínio do R se pode fazer gráficos espetaculares, como por
## exemplo esse com variação da tonalidade.

ht <- hist(dus$valor, prob=TRUE, breaks=seq(35000, 75000, 2000),
           xlab="Preço de venda (R$)", 
           ylab="Frequência absoluta", col="#6E0039",
           sub="Amplitude de classe de R$ 2500")
rug(dus$valor) ## Faz risquinhos no eixo x.

plot of chunk unnamed-chunk-2

## Destacar a barra da classe modal usando outra cor.
wm <- which.max(ht$counts)
cols <- rep("yellow", length(ht$counts))
cols[wm] <- "red"
cols

plot(ht, col=cols)

## Traçar os segmentos que indicam o valor interpolado para a moda.
ycoor <- with(ht, counts[wm+0:1])
xcoor <- with(ht, breaks[wm+0:1])
segments(xcoor[1], ycoor[1], xcoor[2], ycoor[2])

ycoor <- with(ht, counts[wm-1:0])
xcoor <- with(ht, breaks[wm+0:1])
segments(xcoor[1], ycoor[1], xcoor[2], ycoor[2])

## Por semelhança de triangulos a moda obtida é:
ac <- with(ht, diff(breaks[1:2]))
d <- with(ht, abs(diff(counts[wm+(-1:1)])))
xmoda <- with(ht, breaks[wm]+(ac*d[1])/sum(d)); xmoda

abline(v=mean(dus$valor))
abline(v=xmoda, col="yellow")

plot of chunk unnamed-chunk-2

## Como aprimorar um histograma.
plot(ht, col=NULL, lty=0, ann=FALSE, axes=FALSE)
abline(h=seq(0, 100, by=10), lty=2)
plot(ht, col=cols, ann=FALSE, axes=FALSE, add=TRUE)
rug(dus$valor)
axis(side=1, at=seq(35000, 75000, 5000))
axis(side=2, at=seq(0, 100, by=10))
box(bty="L")
title(main="Histograma do valor (R$)",
      sub="Dados retirados do webmotors.com",
      xlab="Valor (R$)",
      ylab="Frequência absoluta")
mtext(side=3, line=0,
      text=paste("Amostra de tamanho", length(dus$valor)))
mtext(side=4, line=-1, col="gray70", outer=TRUE, adj=0,
      text="Feito por Walmes Zeviani")
legend("topright", fill="red", legend="Classe modal", bty="n")

plot of chunk unnamed-chunk-2

## Outra variação de um histograma.
ht <- hist(dus$valor, seq(35000, 75000, 2000), plot=FALSE)
nc <- length(ht$mids)             ## Número de classes.
ac <- diff(ht$breaks[1:2])        ## Amplitude de classe.
ma <- mean(dus$valor)             ## Média da amostra.
md <- median(dus$valor)           ## Mediana da amostra.
qts <- fivenum(dus$valor)[c(2,4)] ## 1Q e 3Q da amostra.
modal <- which.max(ht$counts)     ## Classe modal.
modal <- list(x=ht$mids[modal], y=ht$counts[modal])
colseq <- rgb(red=0.25, blue=0.7,
              green=seq(0.1, 0.9, length.out=nc))

plot(ht, col=colseq, ylim=c(0, modal$y+strheight("1")),
     xlab="Preço de venda (R$)",
     ylab="Frequência absoluta",
     sub=paste("Amplitude de classe de R$", ac),
     main=NULL, border="gray50")
text(x=modal$x, y=modal$y, labels=modal$y, pos=3)
rug(dus$valor)
arrows(ma, 0, ma, modal$y/3, code=1, length=0.15)
text(ma, modal$y/3, labels=paste("Média:", round(ma,2)), pos=3)
arrows(md, 0, md, modal$y/6, code=1, length=0.15)
text(ma, modal$y/6, labels=paste("Mediana:", round(md,1)),
     pos=ifelse(md<ma, 2, 4))
box()

plot of chunk unnamed-chunk-2

## Responda: o que de informação foi acrescentado com as barras mudando
## de cor? Alguns não gostam de distração ou poluição visual, outros
## acreditam que isso atrai o leitor.

##-----------------------------------------------------------------------------
## Gráficos de densidade.

den <- density(dus$valor)
str(den)

plot(den)

plot of chunk unnamed-chunk-2

modal <- which.max(den$y)
modal <- list(x=den$x[modal], y=den$y[modal])

plot(den, type="n", xlab="Preço de venda (R$)", ylab="Densidade",
     ylim=c(0, modal$y+strheight("1")), main="",
     sub=paste("Bandwidth:", round(den$bw,3)))
with(den, polygon(x, y, col="gray90"))
with(modal, segments(x, 0, x, y, col=2))
with(modal, text(x, y, labels=round(x, 2), pos=3))
arrows(ma, 0, ma, modal$y/3, code=1, length=0.15)
text(ma, modal$y/3, labels=paste("Média:", round(ma,2)), pos=3)
arrows(md, 0, md, modal$y/6, code=1, length=0.15)
text(ma, modal$y/6, labels=paste("Mediana:", round(md,1)),
     pos=ifelse(md<ma, 2, 4))
rug(dus$valor)

plot of chunk unnamed-chunk-2

## Frequência acumulada empírica.
y <- ecdf(dus$valor)
plot(y)

plot of chunk unnamed-chunk-2

plot(y, xlab="Preço de venda (R$)",
     ylab="Frequência relativa acumulada",
     cex=NA, verticals=TRUE, main=NULL)

## Destacando a frequência de veículos com preço de 50 à 60 mil.
lim <- c(50000,60000)
ptbl <- prop.table(table(cut(dus$valor,
                             breaks=c(-Inf,lim,Inf))))
cs <- cumsum(ptbl)[seq_along(lim)]

plot(y, xlab="Preço de venda (R$)",
     ylab="Frequência relativa acumulada",
     cex=NA, verticals=TRUE, main=NULL)
segments(lim, 0, lim, cs, lty=2)
segments(lim, cs, par()$usr[3], cs, lty=2)
arrows(lim[1], cs[1], lim[1], cs[2], code=3, length=0.15)
text(lim[1], median(cs), labels=round(ptbl[2], 4),
     srt=90, adj=c(0.5,-0.5))
rug(dus$valor)

plot of chunk unnamed-chunk-2

##-----------------------------------------------------------------------------
## Diagrama de dispersão.

plot(valor~km, data=dus)

plot of chunk unnamed-chunk-2

plot(valor~km, data=dus,
     xlab="Distância percorrida (km)",
     ylab="Preço de venda (R$)")

## Adicionar uma linha de tendência suave.
plot(valor~km, data=dus,
     xlab="Distância percorrida (km)",
     ylab="Preço de venda (R$)")
with(dus, lines(lowess(x=km, y=valor), lwd=2))

plot of chunk unnamed-chunk-2

## Usar cores diferentes para identificar o tipo de câmbio, com linhas
## de tendência e grid.
plot(valor~km, data=dus, col=c(2,4)[dus$cambio],
     xlab="Distância percorrida (km)",
     ylab="Preço de venda (R$)")
with(subset(dus, cambio=="AUTOMÁTICO"),
     lines(lowess(x=km, y=valor), col=2, lwd=1.5,))
with(subset(dus, cambio=="MANUAL"),
     lines(lowess(x=km, y=valor), col=4, lwd=1.5,))
legend("topright", lty=1, col=c(2,4), lwd=1.5,
       legend=levels(dus$cambio), bty="n")
grid()

plot of chunk unnamed-chunk-2

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

boxplot(valor~ano, data=dus)

plot of chunk unnamed-chunk-2

boxplot(valor~cat, data=dus)

plot of chunk unnamed-chunk-2

levels(dus$cat)

levels(dus$cat) <- c("?", "Dynamique", "Expression",
                     "Tech Road I", "Tech Road II")

boxplot(valor~cat, data=dus,
        xlab="Modelo", ylab="Preço de venda (R$)")

plot of chunk unnamed-chunk-2

## Larguras proporcionais à raiz da quantidade em cada grupo.
pal <- c("#583882","#35165F","#43256C","#705199","#8A71AA")
boxplot(valor~cat, data=dus, varwidth=TRUE, pars=list(boxwex=1.25),
        col=pal, xlab="Modelo", ylab="Preço de venda (R$)")

plot of chunk unnamed-chunk-2

## Indicação do valor da média.
mds <- with(dus, tapply(valor, cat, mean))

boxplot(valor~cat, data=dus,
        xlab="Modelo", ylab="Preço de venda (R$)")
points(x=1:nlevels(dus$cat), y=mds, pch=15, cex=1.5)

plot of chunk unnamed-chunk-2

##-----------------------------------------------------------------------------
## Gráficos com o valor para a média e barra de erro para o
## desvio-padrão.

res <- aggregate(valor~cat, data=dus,
                 FUN=function(x){ c(m=mean(x), s=sd(x)) })

## Criando os limites superior e inferior.
res <- transform(res,
                 lwr=valor[,1]-valor[,2],
                 upr=valor[,1]+valor[,2],
                 catf=as.integer(cat))
res

with(res, matplot(x=catf, y=cbind(lwr, upr),
                  xlim=extendrange(x=catf, f=0.1),
                  ylim=extendrange(x=c(lwr, upr), f=0.1),
                  type="n", ann=FALSE, xaxt="n"))
grid()
with(res, points(x=catf, y=valor[,"m"], pch=19))
with(res, arrows(catf, lwr, catf, upr, code=3,
                 angle=90, length=0.1))
with(res, axis(side=1, at=catf, labels=as.character(cat)))
title(xlab="Categoria", ylab="Valor (R$)")
mtext(side=3, line=0,
      text=expression("Barras de erro representam "*bar(x) %+-% 1*s))

plot of chunk unnamed-chunk-2

## Com o boxplot ao lado.
with(res, matplot(x=catf, y=cbind(lwr, upr),
                  xlim=extendrange(x=catf, f=0.1),
                  ylim=extendrange(x=c(lwr, upr, dus$valor), f=0.1),
                  type="n", ann=FALSE, xaxt="n"))
grid()
with(res, points(x=catf, y=valor[,"m"], pch=19))
with(res, arrows(catf, lwr, catf, upr, code=3,
                 angle=90, length=0.1))
with(res, axis(side=1, at=catf, labels=as.character(cat)))
title(xlab="Categoria", ylab="Valor (R$)")
mtext(side=3, line=0,
      text=expression("Barras de erro representam "*bar(x) %+-% 1*s))
boxplot(valor~cat,
        at=1:nlevels(dus$cat)+0.2, col="gray45",
        data=dus, add=TRUE, ann=FALSE, axes=FALSE,
        pars=list(boxwex=0.1))

plot of chunk unnamed-chunk-2

## Com boxplot e pontos dispersos dos lados.

xlim <- extendrange(x=1:nlevels(dus$cat), f=0.1)
ylim <- extendrange(x=c(res$lwr, res$upr, dus$valor), f=0.1)

par(mar=c(5.1,4.1,4.1,0))
layout(matrix(c(1,2), ncol=2), widths=c(0.85,0.15))
with(dus,
     plot.default(
         x=jitter(as.integer(cat), factor=0.25)-0.2,
         y=valor, xaxt="n", ann=FALSE, col="gray50",
         xlim=xlim, ylim=ylim))
grid()
with(res, points(x=catf, y=valor[,"m"], pch=19))
with(res, arrows(catf, lwr, catf, upr, code=3,
                 angle=90, length=0.05))
with(res, axis(side=1, at=catf, labels=as.character(cat), cex.axis=0.95))
title(xlab="Categoria", ylab="Valor (R$)")
mtext(side=3, line=0,
      text=expression("Barras de erro representam "*bar(x) %+-% 1*s))
boxplot(valor~cat,
        at=1:nlevels(dus$cat)+0.2, col="gray45",
        data=dus, add=TRUE, ann=FALSE, axes=FALSE,
        pars=list(boxwex=0.1))
par(mar=c(5.1,0.1,4.1,1))
yhist <- hist(dus$valor, plot=FALSE, breaks=20)
str(yhist)
with(yhist,
     plot(x=NULL, y=NULL, ann=FALSE, axes=FALSE,
          ylim=ylim, xlim=c(0, max(density))))
rug(side=2, dus$valor)
snc <- 1:length(yhist$mids)
with(yhist,
     rect(0, breaks[snc],
          density[snc], breaks[snc+1],
          col="gray70"))
den <- density(dus$valor)
with(den, lines(x=y, y=x, col="red", lwd=2))

plot of chunk unnamed-chunk-2


print(sessionInfo(), locale=FALSE)
## R version 3.1.1 (2014-07-10)
## Platform: i686-pc-linux-gnu (32-bit)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  base     
## 
## other attached packages:
## [1] rmarkdown_0.3.3 knitr_1.6      
## 
## loaded via a namespace (and not attached):
## [1] digest_0.6.4    evaluate_0.5.5  formatR_0.10    htmltools_0.2.6 stringr_0.6.2  
## [6] tools_3.1.1
Sys.time()
## [1] "2014-10-08 15:16:42 BRT"