✎ ridiculas - Wiki do LEG

[ Por favor, preencha este local (:wiki:site_notice) ]

Você está aqui: start » ridiculas

Você carregou uma revisão antiga desse documento! Se você salvá-la, irá criar uma nova versão com esses dados.

===== R-idículas - dicas curtas sobre R =====

**//R-idículas//** é a página do LEG dedicada à fornecer //dicas curtas// sobre R, e.g. entre outras: otimização do ambiente, condução de análises, operação com dados e confecção de gráficos. As dicas estão organizadas pelo título, seguido de descrição, palavras-chave e CMR (código mínimo reproduzível). Se você deseja contribuir com a nossa página de R-idículas, envie e-mail para ''walmes@ufpr.br''.

----

==== Controlando gráficamente parâmetros de distribuições de probabilidade ====

Temporariamente sem descrição.\\
palavras-chave: #rpanel.

require(rpanel)

norm.panel <- function(panel){
  ##------------------------------------------------------------------
  ## panel$interval: vetor com domínio de plotagem da função
  ## penel$...: serão parâmetros da distribuição de probabilidades
  curve(dnorm(x, mean=panel$mean, sd=panel$sd),
        from=panel$interval[1], to=panel$interval[2])
  panel
  ##------------------------------------------------------------------
}

# passar os argumentos que serão fixos, abre a janelinha
panel <- rp.control(interval=c(-4,4))

# controla a média
rp.slider(panel, mean, -4, 4, initval=0, showvalue=TRUE, action=norm.panel)

# controla o desvio-padrão
rp.slider(panel, sd, 0.001, 10, initval=1, showvalue=TRUE, action=norm.panel)

#------------------------------------------------------------------------------------------
</code>

----

==== .Rprofile no Linux ====
<code R>
#------------------------------------------------------------------------------------------
#                                                                               por JCFaria
#------------------------------------------------------------------------------------------
</code>

Esse post tem a finalidade de compartilhar algumas coisas que considero importantes na inicialização do R no Linux!

Muitas das opções importantes, do ponto de vista funcional (não relativos à aparência), podem ser feitas no arquivo .Rprofile.
Esse arquivo deve ficar localizado no home do usuário (~/.Rprofile) e é um dos primeiros a ser lido quando uma sessão do
R é iniciada.

Tenho uma função (bem simples) que uso bastante em meu dia a dia: "cv" para calcular o coef. de variação de uma ANOVA:

Pois bem, a medida que vamos aumentando nossa intimidade com o R (inevitavelmente) iremos desenvolvendo nossas próprias
funções (o R foi projetado para isso).

<fc #000080>Ai vem o problema: ter que sempre carregar a função quando for usar, o que pode se tornar uma chatisse!
Pior ainda, ao limpar o workspace do usuário (.GlobalEnv) ela é removida e precisa ser recarregada. Imagine o contexto da necessidade de várias funções, que é o que geralmente acontece com o passar do tempo ...</fc>

<fc #000080>**Tem como contornar? Sim! De várias formas:**</fc>

**//1. Opção muito pouco prática://**
<code R>
oldp <- getwd()
setwd('/home/jcfaria/dados/r/funcoes/') 
source('cv.r')
setwd(oldp)
</code>

A função "cv" ficará disponível no meu workspace mas será removida com a instrução:

muito usada por várias GUIs.

**//2. Opção mais prática://**
A mesma que a anterior, contudo, a função não deverá se chamar "cv", mas sim ".cv".
Nesse caso ela permanecerá como um objeto oculto no meu workspace e não será removida com a intrução:

Contudo, poderá ser removida com a intrução:

**//3. Colocando suas funções em algum ambiente (environment) do R (optei pelo base)://**
<code R>
oldp <- getwd()
setwd('/home/jcfaria/dados/r/funcoes/') 
source('cv.r', local=baseenv())
setwd(oldp)
</code>

Ela não ficará no meu workspace, mas sim no base.
Como tal, poderá ser usada com qualquer outra função desse pacote.

**//4. Criando seu próprio ambiente (acho a solução mais elegante)://**
<code R>
oldp <- getwd()
setwd('/home/jcfaria/dados/r/funcoes/') 
.jcf <- new.env()
source('cv.r', local=.jcf)
setwd(oldp)
</code>

Nesse último caso [**//4//**]:
  * O objeto ".jcf" ficará oculto no meu workspace evitando ser deletado com: rm(list=ls())
  * Parar acessar a função "cv" é necessário:

<code R>
> .jcf$cv

# ou

> with(.jcf, cv)
</code>

Por exemplo:
<code R>
> av <- aov(Sepal.Length ~ Species, data=iris)
> .jcf$cv(av)
[1] 8.81

# ou

> with(.jcf, cv(av))
[1] 8.81
</code>

Esta forma de carregar funções de forma permanente no R pode ser usado para qualquer outro objeto!

No Windows bastava usar no /etc/Rprofile.site:

<code R>
source('cv', local=TRUE)
</code>

que ela ficava disponível no pacote base.

No Windows não testei na versão em desenvolvimento (instável) que uso no Linux,
mas deve funcionar, pois usei por muitos anos em várias versões.

Abaixo meu .Rprofile:
<code R>
## José Cláudio Faria/UESC/DCET

##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬
## General options
library(utils)
options(list(repos='http://cran-r.c3sl.ufpr.br/',
             width=80,
             editor='vim',
             browser='chromium'))

##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬
## Package debug: options
options(list(debug.height=10,
             debug.width=125,
             debug.screen.pos='-1680+00' # dois monitores!
             ))

##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬
## Funtions: myself
oldp <- getwd()
setwd('/home/jcfaria/dados/r/funcoes/')
#source('fitreg.r',    local=baseenv())
#source('fitregl.r',   local=baseenv())
#source('fitrlcor.r',  local=baseenv())
#source('plotreg.r',   local=baseenv())
#source('plotrl.r',    local=baseenv())
#source('plotlm3d.r',  local=baseenv())
#source('cv.r',        local=baseenv())
#source('rs.r',        local=baseenv())
#source('sourcedir.r', local=baseenv())

.jcf <- new.env()
source('fitreg.r',    local=.jcf)
source('fitregl.r',   local=.jcf)
source('fitrlcor.r',  local=.jcf)
source('plotreg.r',   local=.jcf)
source('plotrl.r',    local=.jcf)
source('plotlm3d.r',  local=.jcf)
source('cv.r',        local=.jcf)
source('rs.r',        local=.jcf)
source('sourcedir.r', local=.jcf)
#attach(.jcf)

##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬
## Removing variables
setwd(oldp)
rm(oldp)

##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬
## Packages: loading
#library(fdth)
#library(ScottKnott)
#library(bpca)
#library(TukeyC)
#library(reshape)
#library(debug)
</code>

----

==== Como fazer a justaposição de vários data.frames ====

Temporariamente sem descrição.\\
palavras-chave: #merge, #Reduce.

id <- 1:30  # número único que identifica os registros
n <- 20     # número de registros por data.frame

a1 <- data.frame(id=sample(id, n), v1=rnorm(n))     # resposta 1
a2 <- data.frame(id=sample(id, n), v2=rpois(n,10))  # resposta 2
a3 <- data.frame(id=sample(id, n), v3=runif(n))     # resposta 3

merge(a1, a2, by="id") # justapõe 2 data.frames de cada vez

a0 <- list(a1, a2, a3) # cria uma lista com todos os data.frames

Reduce(function(x, y) merge(x, y, by="id"), a0, accumulate=FALSE) # justapõe todos

#------------------------------------------------------------------------------------------
# font: http://rwiki.sciviews.org/doku.php?id=tips:data-frames:merge
#       https://stat.ethz.ch/pipermail/r-help/2008-April/160836.html
#       http://econometricsense.blogspot.com/2011/01/merging-multiple-data-frames-in-r.html
#       http://www.youtube.com/watch?v=E4uR5I1uLFM

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

mergeAll <- function(..., by="date", all=TRUE){
  dotArgs <- list(...)
  Reduce(
         function(x, y){
           merge(x, y, by = by, all = all, suffixes=paste(".", names(dotArgs), sep = ""))
         },
         dotArgs)
}

mergeAll(a = a1, b = a2, c = a3, by="id")

str(.Last.value)

You also might be able to set it up to capture names without you having to put "a = a" etc. using substitute.

# http://r.789695.n4.nabble.com/merge-multiple-data-frames-td4331089.html
# ver reshape::merge_all()
#------------------------------------------------------------------------------------------
</code>

----

==== Gráfico de valores observados e curva de valores preditos para modelo linear generalizado ====

Temporariamente sem descrição.\\
palavras-chave: #glm, #poisson, #predict.

da <- expand.grid(trat=gl(2,1), tempo=1:20)
y <- rpois(nrow(da), lambda=da$tempo/5)

g0 <- glm(y~trat/(tempo+I(tempo^2)), data=da, family=poisson)

new <- expand.grid(trat=gl(2,1), tempo=seq(1,20,l=50))
new$p0 <- predict(g0, newdata=new, type="response")

plot(y~tempo, da, col=da$trat)
with(subset(new, trat=="1"), lines(p0~tempo, col=1))
with(subset(new, trat=="2"), lines(p0~tempo, col=2))

#------------------------------------------------------------------------------------------
</code>

----

==== Gráficos de barras com intervalos de confiança para as médias ====

Temporariamente sem descrição.\\
palavras-chave: #intervalo, #erro_padrão, #barras.

da <- data.frame(trat=gl(5,8))
da$y <- as.numeric(da$trat)+rnorm(nrow(da))

m0 <- lm(y~trat, da)

new <- data.frame(trat=levels(da$trat))
new$pred <- predict(m0, newdata=new, interval="confidence")
str(new)

ylim <- c(0, max(new$pred)*1.05)
bp <- barplot(new$pred[,1], ylim=ylim)
arrows(bp, new$pred[,2], bp, new$pred[,3], code=3, angle=90)
box()

#------------------------------------------------------------------------------------------
# outras referências

browseURL("http://addictedtor.free.fr/graphiques/graphcode.php?graph=54")
browseURL("http://addictedtor.free.fr/graphiques/RGraphGallery.php?graph=72")

#------------------------------------------------------------------------------------------
</code>

----

==== Gráficos do R exportados pelo dispositivo tikz ====

Temporariamente sem descrição.\\
palavras-chave: #tikz, #sweave, #latex.

<code>
\documentclass{article}

\usepackage{Sweave}
\usepackage{tikz}

\SweaveOpts{keep.source=true}

\title{Usando ti\textit{k}z no Sweave}
\author{Walmes Zeviani\\ LEG/UFPR}

\begin{document}

\maketitle

C\'{o}digo m\'{i}nimo reproduz\'{i}vel usando \texttt{tikzDevice} para exportar gr\'{a}ficos feitos no R
para codifica\c{c}\~{a}o ti\textit{k}z.

{\footnotesize
<<results=hide>>=
#------------------------------------------------------------------------------------------
require(tikzDevice)
set.seed(2011); x <- rnorm(200)
tikz("plot.tex", w=5, h=3)
hist(x, freq=FALSE, ylab="Densidade",
     main="Histograma de uma amostra de $X \\sim N(\\mu=0, \\sigma^2=1)$")
curve(dnorm(x), col=2, add=TRUE, lwd=2); rug(x)
legend("topleft", col=2, lty=1, lwd=2, bty="n",
       legend="$\\displaystyle \\frac{1}{\\sqrt{2\\pi\\sigma^2}}\\cdot e^{-\\frac{(x-\\mu)^2}{2\\sigma^2}}$")
box(); dev.off()
#------------------------------------------------------------------------------------------
@
}

\input{plot.tex}

\end{document}
</code>

----

==== Gráfico com dois eixos coordenados ====

Temporariamente sem descrição.\\
palavras-chave: #série, #eixo.

<code R>
#------------------------------------------------------------------------------------------
#                                                               por Ivan, Benilton e Walmes
#------------------------------------------------------------------------------------------
# dados de séries de dados indexadas no tempo (meses)

lines <- 'meses  temp umidade    rad   chuva
    Jan    26.49   86.58    795.88  0.36
    Fev    26.65   88.49    710.24  0.34
    Mar    27.19   86.16    772.99  0.21
    Abr    26.28   89.75    574.88  0.67
    Mai    26.62   89.22    614.02  0.31
    Jun    26.13   87.83    680.08  0.26
    Jul    25.83   86.57    675.97  0.15
    Ago    27.05   83.14    756.44  0.07
    Set    27.60   83.02    925.57  0.14
    Out    27.44   85.16    927.71  0.17
    Nov    26.56   88.18    788.87  0.19
    Dez    25.87   90.63    703.94  0.33'
da <- read.table(textConnection(lines), header=TRUE)
str(da)

#------------------------------------------------------------------------------------------
# gráfico de duas séries no mesmo gráfico com um eixo coordenado para cada

par(mar=c(5,5,4,5))
plot(da$temp, type="b", pch=15, lwd=1.5,
     xlab="Meses (Anos 2010)", main="Temperatura e umidade em 2010",
     ylab=expression(Temperatura~group("(", degree*C, ")")),
     ylim=c(24,30), xlim=c(1,12), axes=FALSE)
axis(1, at=1:12, labels=da$meses)
axis(2)
par(new=TRUE)
plot(da$umidade, type="b", pch=14, axes=FALSE, frame=TRUE, ann=FALSE)
axis(4)
mtext(text="Umidade (%)", 4, line=3)
legend("topleft", bty="n", seg.len=3, lty=1, lwd=c(1.5,1), pch=c(15,14),
       legend=c("Temperatura","Umidade"), merge=TRUE, trace=FALSE)

#------------------------------------------------------------------------------------------
# funções sugeridas

apropos("month")
demo(plotmath)

#------------------------------------------------------------------------------------------
</code>

----

==== Análise de dados de proporção usando modelo linear generalizado ====

Temporariamente sem descrição.\\
palavras-chave: #binomial, #sucessos, #deviance, #wireframe, #fatorial, #superfície.

rend <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/rendimento.txt", header=TRUE)
rend <- transform(rend, k=factor(K), a=factor(A), bloc=factor(bloc))
str(rend)

#------------------------------------------------------------------------------------------
# ajuste modelo de caselas aos dados assumindo distribuição binomial (link=logit)

g0 <- glm(cbind(nv, nvi)~bloc+k*a, data=rend, family=binomial)

#------------------------------------------------------------------------------------------
# análise de resíduos usual para verificar anomalias

par(mfrow=c(2,2))
plot(g0)
layout(1)

#------------------------------------------------------------------------------------------
# quadro de estimativas e quadro de análise de deviance, faz a vez da anova

summary(g0)
anova(g0, test="Chisq")
  
#------------------------------------------------------------------------------------------
# obter modelo mais parcimonioso, usar fatores na forma contínua

g1 <- glm(cbind(nv, nvi)~bloc+K+A+I(K^2)+I(A^2)+K:A, data=rend, family=binomial)
summary(g1) # comparar deviance residual com grau de liberdade residual
g1 <- update(g1, formula=.~.-K:A, family=quasibinomial)
summary(g1)
anova(g1, test="F")
  
#------------------------------------------------------------------------------------------
# faz a predição dos valores (usando o apenas o bloco 1)

pred <- with(rend,
             expand.grid(A=seq(min(A),max(A),l=20),
                         K=seq(min(K),max(K),l=20),
                         bloc="1"))
pred$prob <- predict(g1, newdata=pred, type="response")
 
#------------------------------------------------------------------------------------------
# gráfico

require(lattice)
wireframe(prob~A+K, data=pred,
          zlab=list("Probabilidade de germinação", rot=90),
          scales=list(arrows=FALSE),
          screen=list(z=-50, x=-60), drape=TRUE)

#------------------------------------------------------------------------------------------
</code>

----

==== Desdobramento de interação em experimento fatorial ====

Temporariamente sem descrição.\\
palavras-chave: #desdobramento, #comparações, #teste_de_médias, #fatorial.

rend <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/rendimento.txt", header=TRUE)
str(rend)
rend <- transform(rend, K=factor(K), A=factor(A), bloc=factor(bloc))
str(rend)
  
#------------------------------------------------------------------------------------------
# análise gráfica

require(lattice)
xyplot(rg~K|A, groups=bloc, data=rend, type="b", auto.key=TRUE)
  
#------------------------------------------------------------------------------------------
# ajuste do modelo

m0 <- aov(rg~bloc+A*K, data=rend)
summary(m0)
  
#------------------------------------------------------------------------------------------
# checagem

par(mfrow=c(2,2)); plot(m0); layout(1)
  
#------------------------------------------------------------------------------------------
# desdobrando somas de quadrados para a variação de K dentro de A

m1 <- aov(rg~bloc+A/K, data=rend)
summary(m1)
summary(m1, split=list("A:K"=list(
                         "A-37.5"=c(1,4,7,10),
                         "A-50.0"=c(2,5,8,11),
                         "A-62.5"=c(3,6,9,12)
                         )))
  
#------------------------------------------------------------------------------------------
# para facilitar encontrar as posições pode-se fazer a busca por expessões regulares

names(coef(m1))[8:19]
grep("A37.5", names(coef(m1))[8:19])
grep("A50", names(coef(m1))[8:19])
grep("A62.5", names(coef(m1))[8:19])
  
#------------------------------------------------------------------------------------------
# usando as expressões regulares vamos desdobrar A dentro de K

m2 <- aov(rg~bloc+K/A, data=rend)
summary(m2)
names(coef(m2))
  
#------------------------------------------------------------------------------------------
# buscando pela expressão regular

grep("K0", names(coef(m2))[10:19])
desAinK <- sapply(paste("K", levels(rend$K),sep=""), simplify=FALSE,
                  grep, x=names(coef(m2)[10:19]))

#------------------------------------------------------------------------------------------
# decomposição das somas de quadrados

summary(m2, split=list("K:A"=desAinK))

#------------------------------------------------------------------------------------------
# usando o ExpDes (inglês) (https://sites.google.com/site/ericbferreira/unifal/downloads-1)

require(ExpDes)
help(package="ExpDes")
help(fat2.rbd, help_type="html")
  
#------------------------------------------------------------------------------------------
# aplicando a função do Eric fat2.rbd

with(rend, fat2.rbd(A, K, bloc, rg, mcomp="sk", quali=c(TRUE, TRUE)))
  
#------------------------------------------------------------------------------------------
# desdobrando a interação em testes de Tukey para níveis de K fixando os níveis de A

require(agricolae)

with(subset(rend, A=="37.5"),
     HSD.test(rg, K, DFerror=df.residual(m0), MSerror=deviance(m0)/df.residual(m0)))
with(subset(rend, A=="50"),
     HSD.test(rg, K, DFerror=df.residual(m0), MSerror=deviance(m0)/df.residual(m0)))
with(subset(rend, A=="62.5"),
     HSD.test(rg, K, DFerror=df.residual(m0), MSerror=deviance(m0)/df.residual(m0)))
  
#------------------------------------------------------------------------------------------
# usando funções para fazer o desdobramento (lapply)

levels(rend$A)
sapply(levels(rend$A), simplify=FALSE,
       function(a){
         with(subset(rend, A==a),
              HSD.test(rg, K,
                       DFerror=df.residual(m0),
                       MSerror=deviance(m0)/df.residual(m0)))
       })
  
#------------------------------------------------------------------------------------------
# fazendo o mesmo para o teste ScottKnott (a ordem A*K e K*A é importante!)

require(ScottKnott)
sk <- SK.nest(x=rend, y=rend$rg, model="y~bloc+A*K", which="A:K", fl2=1)
summary(sk)
sk <- SK.nest(x=rend, y=rend$rg, model="y~bloc+K*A", which="K:A", fl2=1)
summary(sk)
  
#------------------------------------------------------------------------------------------
# fazer o teste de ScottKnott com um comando apenas (lapply)

levels(rend$A)
lapply(seq_along(levels(rend$A)),
       function(a){
         sk <- SK.nest(x=rend, y=rend$rg, model="y~bloc+K*A", which="K:A", fl2=a)
         summary(sk)
       })
  
#------------------------------------------------------------------------------------------

levels(rend$K)
lapply(seq_along(levels(rend$K)),
       function(a){
         sk <- SK.nest(x=rend, y=rend$rg, model="y~bloc+A*K", which="A:K", fl2=a)
         summary(sk)
       })

#------------------------------------------------------------------------------------------
# uma forma mais interessante de ajeitar o resultado

require(plyr)

aux <- sapply(levels(rend$A), simplify=FALSE,
              function(a){
                with(subset(rend, A==a),
                     HSD.test(rg, K,
                              DFerror=df.residual(m0),
                              MSerror=deviance(m0)/df.residual(m0)))
              })

tTukey <- ldply(aux, NULL)
tTukey$M <- gsub(" ", "", tTukey$M, fixed=TRUE)
tTukey$trt <- as.factor(as.numeric(as.character(tTukey$trt)))
str(tTukey)

#------------------------------------------------------------------------------------------
# a cereja em cima do bolo como diz o PJ

barchart(means~trt|.id, data=tTukey, horiz=FALSE, layout=c(3,1),
         panel=function(x, y, subscripts, ...){
           panel.barchart(x, y, subscripts=subscripts, ...)
           panel.text(x, y, label=tTukey[subscripts,"M"], pos=3)
         })

#------------------------------------------------------------------------------------------
</code>

----

==== Gráfico para comportamento do desempenho em avaliações ====

Temporariamente sem descrição.\\
palavras-chave: .

p1 <- 10*rbeta(60, 7, 4) # notas na primeira avaliação

stati <- function(x){    # função para obter algumas estatísticas
  x <- na.omit(x); nobs <- length(x); xbar <- mean(x); md <- median(x)
  max <- max(x); min <- min(x); sd <- sd(x); fv <- fivenum(x)
  M <- data.frame(n=nobs, média=xbar, mediana=md, máximo=max, mínimo=min,
                  q1=fv[2], q3=fv[4], des.pad=sd)
  round(M,2)
}

texto <- paste(capture.output(t(as.matrix(stati(p1))))[-1], collapse="\n")
rgb.palette <- colorRampPalette(c("red","yellow"), space="rgb") # escala de cores
op <- par()

par(mar=c(5.1,4.1,0,0), fig=c(0,0.7,0,0.7))
plot(ecdf(p1), main=NULL, xlim=c(0,10), col=1, lwd=2,
     col.01line=NULL, xaxt="n", verticals=TRUE, cex=0, las=2,
     ylab="Distribuição acumulada empírica", xlab="Nota")
axis(1, at=0:10)
rug(p1)
abline(v=0:10, h=seq(0,1,by=0.1), col="gray75")
abline(v=7, lty=3)
par(family="mono")
text(0,1, label=texto, adj=c(0,1))

par(op)
par(mar=c(5.1,0,0,3.1), new=TRUE, fig=c(0.7,1,0,0.7))
boxplot(p1, notch=TRUE, yaxt="n", ylim=c(0,10), col="orange")
rug(p1, side=4, ticksize=0.07)
abline(h=7, lty=3)
axis(4, at=0:10, las=2)
points(1, mean(p1))

par(mar=c(0,4.1,1.5,0), new=TRUE, fig=c(0,0.7,0.7,1))
ht <- hist(p1, freq=FALSE, xaxt="n", xlim=c(0,10), main=NULL,
           ylab="Densidade", col=rgb.palette(10), las=2)
text(ht$mids, ht$density, ht$counts, pos=1)
rug(p1)
abline(v=7, lty=3)
lines(density(p1, from=0, to=10), col=1, lwd=2)
box()

#------------------------------------------------------------------------------------------
</code>

----

==== Desdobramento de interação usando a multcomp::ghlt() ====

Temporariamente sem descrição.\\
palavras-chave: #desdobramento, #glht, #comparações, #teste_de_médias, #fatorial, #contrastes.

da <- expand.grid(A=gl(3,4, labels="a"), B=gl(2,1, labels="b"))
da$y <- rnorm(nrow(da))
da

m0 <- lm(y~A*B, da)
summary(m0)

#------------------------------------------------------------------------------------------
# pacotes necessários

require(contrast)
require(multcomp)

#------------------------------------------------------------------------------------------
# níveis dos fatores

ii <- levels(da$A)
jj <- levels(da$B)

#------------------------------------------------------------------------------------------
# preparar a matriz contrastes para A dentro de B, esse esquema foi um artifício que criei

m.contr <- outer(ii, ii,
                 function(x,y){ paste(x, y, sep="-") })
m.contr <- m.contr[upper.tri(m.contr)]
p.contr <- do.call(rbind, strsplit(m.contr, "-"))
p.contr

#------------------------------------------------------------------------------------------
# lista com as matrizes dos contrastes que específicam o desdobramento A dentro de B

c.X <- sapply(jj, simplify=FALSE,
              function(j){
                sapply(1:nrow(p.contr),
                       function(i){
                         c.contr <- contrast(m0,
                                             list(A=p.contr[i,1], B=j),
                                             list(A=p.contr[i,2], B=j))
                         c.contr$X
                       })})
c.X

#------------------------------------------------------------------------------------------
# usando a glht para desdobrar, cada entrada da lista é um nível de B

mc.contr <- sapply(c.X, simplify=FALSE,
                   function(X){
                     summary(glht(m0, linfct=t(X)))
                   })
mc.contr

#------------------------------------------------------------------------------------------
</code>

----

==== Regressão na análise de variância ====

No exemplo abaixo fazemos a análise de um experimento fatorial. São combinados três níveis de um fator qualitativo (''cultivar'' de sorgo) com 6 níveis de um fator quantitativo (''dose'' de nitrogênio). A resposta avaliada nas parcelas é o índice agronômico. O experimento foi instalado em blocos. A interação entre as fontes de variação foi significativa. O desdobramento dessa interação foi feita nos dois sentidos, ou seja, estudo do fator dose fixando nível do fator cultivar e estudo do fator cultivar fixando o nível do fator dose. Foi feito a ajuste de modelos de regressão polinomial para o estudo do fator dose. São ilustrados passos para teste de hipótese sobre os efeitos linear, quadrático, etc, e falta de ajuste do modelo ajustado, além de como obter as estimativas dos termos do modelo e o coeficiente de determinação (R²).\\
palavras-chave: #regressão, #polinômio, #anova.

<code R>
#------------------------------------------------------------------------------------------
#                                                                                por Walmes
#------------------------------------------------------------------------------------------
# dados
sorgo <- read.table("http://www.leg.ufpr.br/~walmes/docs/anovareg.txt", header=TRUE)
sorgo <- transform(sorgo, bloco=factor(bloco), cultivar=factor(cultivar))
str(sorgo)
                                                                                          #
#------------------------------------------------------------------------------------------
# gráficos exploratórios
require(lattice)
xyplot(indice~dose|cultivar, groups=bloco, data=sorgo,
       jitter.x=TRUE, type=c("p","l"), layout=c(3,1))
xyplot(indice~dose, groups=cultivar, data=sorgo, jitter.x=TRUE, type=c("p","a"))
                                                                                          #
#------------------------------------------------------------------------------------------
# análise de variância do modelo de fatores
m0 <- aov(indice~bloco+cultivar*ordered(dose), data=sorgo)
summary(m0)
                                                                                          #
#------------------------------------------------------------------------------------------
# checagem
par(mfrow=c(2,2))
plot(m0)
layout(1)
                                                                                          #
#------------------------------------------------------------------------------------------
# desdobrando as somas de quadrados de doses dentro de cultivar
# dicas: forneça para ’by’ o número de níveis de cultivar (3)
# forneça para ’length.out’ os graus de liberdade de dose (6-1)
m1 <- aov(indice~bloco+cultivar/ordered(dose), data=sorgo)
summary(m1)
coef(m1)
summary(m1, split=list("cultivar:ordered(dose)"=list(
                         "Ag-1002"=seq(1, by=3, length.out=5),
                         "BR-300"=seq(2, by=3, length.out=5),
                         "Pioneer-B815"=seq(3, by=3, length.out=5)
                         )))
                                                                                          #
#------------------------------------------------------------------------------------------
# desdobrando somas de quadrados de cultivar dentro das doses
# dicas: forneça para ’by’ o número de níveis de dose (6)
# forneça para ’length.out’ os graus de liberdade de cultivar (3-1)
m2 <- aov(indice~bloco+ordered(dose)/cultivar, data=sorgo)
coef(m2)
summary(m2, split=list("ordered(dose):cultivar"=list(
                         "N.0"=seq(1, by=6, length.out=2),
                         "N.60"=seq(2, by=6, length.out=2),
                         "N.120"=seq(3, by=6, length.out=2),
                         "N.180"=seq(4, by=6, length.out=2),
                         "N.240"=seq(5, by=6, length.out=2),
                         "N.300"=seq(6, by=6, length.out=2)
                         )))
                                                                                          #
#------------------------------------------------------------------------------------------
# desdobrando efeitos dos graus polinômio dentro de dose dentro de cultivar
# lof é falta de ajuste (lack of fit)
summary(m1, split=list("cultivar:ordered(dose)"=list(
                         "Ag-1002.L"=1,
                         "Ag-1002.Q"=4,
                         "Ag-1002.C"=7,
                         "Ag-1002.lof"=c(10,13),
                         "BR-300.L"=2,
                         "BR-300.Q"=5,
                         "BR-300.C"=8,
                         "BR-300.lof"=c(11,14),
                         "Pioneer-B815.L"=3,
                         "Pioneer-B815.Q"=6,
                         "Pioneer-B815.C"=9,
                         "Pioneer-B815.lof"=c(12,15)
                         )))
                                                                                          #
#------------------------------------------------------------------------------------------
# obter as equações de regressão e R^2 para os modelos linear, quadrático e cúbico
# dica: usar contraste tipo soma zero para blocos para se anularem na fórmula
# e remover o intercepto especificando o ’-1’, trocar a ordem dos termos no modelo
# linear (estimativas corretas mas erros padrões e p-valores precisam de correção)
m3 <- aov(indice~-1+cultivar/dose+bloco, data=sorgo,
          contrast=list(bloco=contr.sum))
summary.lm(m3)
                                                                                          #
#------------------------------------------------------------------------------------------
# quadrático (estimativas corretas mas erros padrões e p-valores precisam de correção)
m4 <- aov(indice~-1+cultivar/(dose+I(dose^2))+bloco, data=sorgo,
          contrast=list(bloco=contr.sum))
summary.lm(m4)
                                                                                          #
#------------------------------------------------------------------------------------------
# cúbico (estimativas corretas mas erros padrões e p-valores precisam de correção)
m5 <- aov(indice~-1+cultivar/(dose+I(dose^2)+I(dose^3))+bloco, data=sorgo,
          contrast=list(bloco=contr.sum))
summary.lm(m5)
                                                                                          #
#------------------------------------------------------------------------------------------
# calcular os R^2
sapply(c(linear=1, quadrático=2, cúbico=3),
       function(degree){
         sapply(levels(sorgo$cultivar),
                function(i){
                  da <- with(subset(sorgo, cultivar==i),
                             aggregate(indice, list(dose=dose), mean))
                  summary(lm(x~poly(dose, degree, raw=TRUE), da))$r.squared
                })})
                                                                                          #
#------------------------------------------------------------------------------------------
</code>

----

==== Experimento com dois fatores de efeito aditivo e perda de muitas parcelas ====

Apresenta a análise de um experimento de dois fatores qualitativos com a perda de muitas parcelas. Nesse caso, a perda de parcelas não compromete a estimabilidade dos efeitos dos níveis dos fatores, apenas implica que eles sejam estimados com precisões diferentes. Além disso, a ortogonalidade entre os efeitos é perdida e a interpretação dos testes de hipótese pela análise de variância requer cuidados. Finalmente, é ilustrado a obtenção das estimativas das médias marginais para os níveis dos fatores (conhecidas por lsmeans na documentação de outro aplicativo) e realizadas as comparações múltiplas dessas médias com correção para o nível nominal de significância.\\
palavras-chave: #parcela_perdida, #desbalanceamento, #médias_ajustadas, #aditivo.

da <- expand.grid(rept=1:5, ep=factor(1:5), tr=factor(1:4))
da$y <- c(58.4, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
          68.4, NA, NA, NA, NA, 258.8, 265.6, NA, NA, NA, NA, NA, 250, NA, 278.8,
          268.8, NA, NA, NA, 309.6, NA, NA, NA, NA, NA, NA, NA, NA, NA, 254, 598.8,
          NA, NA, NA, NA, 250, 399.6, 260, NA, NA, NA, 288.4, NA, NA, NA, 397.2, NA,
          NA, 337.6, NA, 415.2, NA, 450.8, NA, NA, NA, NA, 393.2, NA, NA, NA, NA,
          NA, NA, NA, 380.4, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 634, 417.2,
          NA, NA, NA, NA, NA)

#------------------------------------------------------------------------------------------
# ajuste do modelo aditivo com teste F marginal

m0 <- lm(y~ep+tr, data=da)
drop1(m0, test="F")

#------------------------------------------------------------------------------------------
# análise gráfica dos resíduos

par(mfrow=c(2,2))
plot(m0)
layout(1)

#------------------------------------------------------------------------------------------
# estimativas dos efeitos (solução) sob a restrição do R

summary(m0)

#------------------------------------------------------------------------------------------
# obtenção das médias ajustadas dos níveis de tratamento

require(contrast)
lapply(levels(da$tr),
       function(i){
         contrast(m0, type="average", list(tr=i, ep=levels(da$ep)))
       }
       )

#------------------------------------------------------------------------------------------
# comparação múltipla de efeitos

require(multcomp)
summary(glht(m0, linfct=mcp(tr="Tukey")))

#------------------------------------------------------------------------------------------
</code>

----

==== Experimento em parcelas subdivididas com resposta do tipo binomial ====

Temporariamente sem descrição.\\
palavras-chave: #binomial, #subdividida, #verossimilhança, #experimento, #glm.

require(MASS)

#------------------------------------------------------------------------------------------
# simulandos com a estrutura de um modelo de parcelas subdivididas e resposta binomial
# melhor forma de entender como são gerados os valores observados num experimento

nA <- 4;                              # número de níveis de um fator A de efeito fixo
nB <- 8;                              # número de níveis de um fator B de efeito aleatório
nC <- 2;                              # número de níveis de um fator C de efeito fixo
da <- expand.grid(A=factor(1:nA), B=factor(1:nB), C=factor(1:nC))
da$parcela <- with(da, factor(paste(A,B,sep="-"))) # níveis de parcela aleatório
XA <- model.matrix(~A-1, data=da)        # matriz dos efeitos fixos para A
XB <- model.matrix(~B-1, data=da)        # matriz dos efeitos aleatórios de B
XC <- model.matrix(~C-1, data=da)        # matriz dos efeitos fixos para C
Xp <- model.matrix(~parcela-1, data=da)  # matriz dos efeitos aleatórios de parcela
bA <- c(-1,0,1,1)                        # vetor de efeitos fixos A
bB <- rnorm(nlevels(da$B),0,1)           # vetor de efeitos aleatórios B
bC <- c(0,-1)                            # vetor de efeitos fixos C
bp <- rnorm(nlevels(da$parcela),0,1)     # vetor de efeitos aleatórios parcela

da$eta <- XA%*%bA+XB%*%bB+XC%*%bC+Xp%*%bp   # vetor de valores preditos
summary(exp(da$eta)/(1+exp(da$eta)))        # vetor das probabilidades preditas
da$yobs <- rbinom(nrow(da),
                  prob=exp(da$eta)/(1+exp(da$eta)),
                  size=100) # vetor de dados observados
summary(da$yobs)
str(da)

#------------------------------------------------------------------------------------------
# estimação de parâmetros via máxima verossimilhança (usando Penalized Quasi Likelihood)

m0 <- glmmPQL(cbind(yobs, 100-yobs)~A*C, random=~1|B/parcela, data=da, family=binomial)
summary(m0)

#------------------------------------------------------------------------------------------
</code>

----

==== Duas ou mais funções em um mesmo gráfico ====

Temporariamente sem descrição.\\
palavras-chave: #gráfico, #função, #curva.

x <- seq(0, 10, by=0.1)
y1 <- 4-0.3*x+0.02*x^2
y2 <- 3+exp(-0.9*x)
matplot(x, cbind(y1, y2), lty=1, col=1:2, type="l")

#------------------------------------------------------------------------------------------
</code>

----

==== Representação de ajuste de curvas de regressão no mesmo gráfico ====

Faremos o ajuste de três curvas de regressão à dados de índice agronômico de sorgo em função da dose de nitrogênio aplicada. Os dados estão classificados para três cultivares de sorgo. Apesar de se tratar de um experimento em blocos, isso será ignorado no exemplo que tem o objetivo de apresentar os comandos para ajustar curvas e representar o ajuste no mesmo gráfico.
palavras-chave: #regressão, #lapply, #curvas, #ajuste, #lines, #subset.

sorgo <- read.table("http://www.leg.ufpr.br/~walmes/docs/anovareg.txt", header=TRUE)
sorgo <- transform(sorgo, bloco=factor(bloco), cultivar=factor(cultivar))
str(sorgo)

#------------------------------------------------------------------------------------------
# usaremos esse banco de dados para ajustarmos curvas de regressão separadas por cultivar.
# iremos ignorar a presença de blocos assumindo amostragem aleatória simples dos dados.

#------------------------------------------------------------------------------------------
# gráfico de dispersão dos dados.

levels(sorgo$cultivar)
par(mfrow=c(1,3))
with(sorgo, sapply(levels(sorgo$cultivar),
                   function(x){
                     plot(indice~dose, subset(sorgo, cultivar==x), main=x)
                   }
                   )
     )

#------------------------------------------------------------------------------------------
# ajustaremos uma curva para cada cultivar. pode-se usar a lapply o tempo todo para isso.

levels(sorgo$cultivar)
mAG <- lm(indice~dose+I(dose^2), data=subset(sorgo, cultivar==levels(sorgo$cultivar)[1]))
mBR <- lm(indice~dose+I(dose^2), data=subset(sorgo, cultivar==levels(sorgo$cultivar)[2]))
mPi <- lm(indice~dose+I(dose^2), data=subset(sorgo, cultivar==levels(sorgo$cultivar)[3]))

lapply(list(mAG, mBR, mPi), summary)
lapply(list(mAG, mBR, mPi), anova)

#------------------------------------------------------------------------------------------
# agora faremos os gráficos dos valores observados com as retas dentro
# jitter() é para "agitar" os valores de dose antes de plotar, apenas para melhorar aspecto

layout(1)
plot(indice~jitter(dose), data=sorgo, col=sorgo$cultivar)
with(subset(sorgo, cultivar==levels(sorgo$cultivar)[1]), lines(predict(mAG)~dose, col=1))
with(subset(sorgo, cultivar==levels(sorgo$cultivar)[2]), lines(predict(mBR)~dose, col=2))
with(subset(sorgo, cultivar==levels(sorgo$cultivar)[3]), lines(predict(mPi)~dose, col=3))

#------------------------------------------------------------------------------------------
</code>
Pode-se usar as funções ''sapply()'' e/ou ''lapply()'' para tornar o código mais eficiente/enxuto.

----

==== Análise intrabloco e interbloco de experimento em blocos incompletos tipo III ====

Para análise intrabloco pode-se usar a função ''agricolae::BIB.test()''.\\
palavras-chave: #bloco_incompleto, #intrabloco, #interbloco, #efeito_aleatório, #Tukey, #média_ajustada.

bib3 <- c(1, 1, 35, 1, 2, 28, 1, 3, 27, 2, 1, 30, 2, 2, 20, 2, 4, 22,
          3, 1, 28, 3, 2, 16, 3, 5, 18, 4, 1, 36, 4, 3, 29, 4, 4, 30,
          5, 1, 29, 5, 3, 19, 5, 5, 22, 6, 1, 25, 6, 4, 16, 6, 5, 19,
          7, 2, 26, 7, 3, 30, 7, 4, 28, 8, 2, 27, 8, 3, 29, 8, 5, 27,
          9, 2, 29, 9, 4, 29, 9, 5, 27, 10, 3, 27, 10, 4, 26, 10, 5, 29)
bib3 <- as.data.frame(matrix(bib3, ncol=3, byrow=TRUE))
names(bib3) <- c("bloc","trat","resp") 
bib3 <- transform(bib3, bloc=factor(bloc), trat=factor(trat))
str(bib3)

#------------------------------------------------------------------------------------------
# número de ocorrência de tratamento com blocos e número de níveis dos fatores

with(bib3, table(trat, bloc))
levels(bib3$trat)
levels(bib3$bloc)

#------------------------------------------------------------------------------------------
# análise intrabloco
#------------------------------------------------------------------------------------------
# soma de quadrados sequencial com tratamentos ajustados aos blocos (intrabloco)

m0 <- lm(resp~bloc+trat, data=bib3)
anova(m0)

#------------------------------------------------------------------------------------------
# checagem das pressuposições

par(mfrow=c(2,2))
plot(m0)
layout(1)

#------------------------------------------------------------------------------------------
# soma de quadrados sequencial com blocos ajustados aos tratamentos

m1 <- lm(resp~trat+bloc, data=bib3)
anova(m1)

#------------------------------------------------------------------------------------------
# soma de quadrados marginal, ambos fatores ajustados um ao outro

drop1(m1, scope=.~., test="F")
car::Anova(m1, type="III")

#------------------------------------------------------------------------------------------
# obtenção das médias ajustadas
# uso da opção de contraste tipo soma zero torna mais simples as operações matriciais
# nesse caso (intercept) representa a média geral

m0 <- lm(resp~bloc+trat, data=bib3,
         contrast=list(bloc=contr.sum, trat=contr.sum))

#------------------------------------------------------------------------------------------
# vamos usar a matriz de contraste dos efeitos dos tratamentos e posições dos efeitos no
# vetor de estimativas

m0$contrast$trat
m0$assign

#------------------------------------------------------------------------------------------
# obtenção das médias ajustadas é o produto matricial abaixo somado à média geral

maju <- c(m0$contrast$trat%*%coef(m0)[m0$assign==2]+coef(m0)[1])
maju

#------------------------------------------------------------------------------------------
# vetor de coeficientes do contraste tr1 vs tr2 e sua estimativa

contr <- rep(c(0,0,1,-1,0), c(1, nlevels(bib3$bloc)-1, 1, 1, nlevels(bib3$trat)-3)); contr
sum(contr*coef(m0))

#------------------------------------------------------------------------------------------
# variância do contraste, como tem equilibrio, todos os contrastes tem mesma variância

v.dif <- contr%*%vcov(m0)%*%contr; v.dif

#------------------------------------------------------------------------------------------
# diferença mínima significativa pelo teste de Tukey

delta <- qtukey(0.95, nlevels(bib3$trat), df=df.residual(m0))*sqrt(0.5*v.dif); delta

#------------------------------------------------------------------------------------------
# aplicação do teste de Tukey

taju <- agricolae::order.stat(levels(bib3$trat), maju, delta)

#------------------------------------------------------------------------------------------
# gráfico de barras com as médias

taju2 <- taju$means
names(taju2) <- taju$trt

bp <- barplot(taju2,
              xlab="Tratamentos", ylab="Variável resposta",
              ylim=c(0, 1.2*max(taju2)))
text(bp, taju2,
     labels=paste(format(taju$means,dig=3), gsub(" ","",taju$M), sep="\n"), pos=3)
mtext(3, line=2,
      text="Comparação múltipla de médias", cex=1.5)
mtext(3, line=0.5,
      text="Médias seguidas de mesma letra não diferem entre si pelo teste de Tukey (5%)")
box()

#------------------------------------------------------------------------------------------
# análise interbloco
#------------------------------------------------------------------------------------------
# ajuste do modelo, anova() não faz anova e sim teste de Wald sequencial para efeitos fixos

require(nlme)

mm0 <- lme(resp~trat, random=~1|bloc, data=bib3)
anova(mm0)

#------------------------------------------------------------------------------------------
# checagem da normalidade efeitos/erros aleatórios

par(mfrow=c(1,2))
qqnorm(residuals(mm0)); qqline(residuals(mm0))
qqnorm(unlist(ranef(mm0))); qqline(unlist(ranef(mm0)))
layout(1)

#------------------------------------------------------------------------------------------
# estimativa dos componetes de variância

VarCorr(mm0)

#------------------------------------------------------------------------------------------
# matriz de contraste dos tratamentos

contr <- mm0$contrasts$trat
idx <- grep("trat", names(fixef(mm0)))

#------------------------------------------------------------------------------------------
# médias ajustadas dos tratamentos

maju <- c(contr%*%fixef(mm0)[idx]+fixef(mm0)[1])
maju

#------------------------------------------------------------------------------------------
# vetor de coeficientes do contraste tr1 vs tr2 e sua estimativa

contr <- rep(c(0, 1, -1, 0), c(1, 1, 1, nlevels(bib3$trat)-3)); contr
contr%*%maju

#------------------------------------------------------------------------------------------
# variância do contraste, como tem equilibrio, todos os contrastes tem mesma variância

v.dif <- contr%*%vcov(mm0)%*%contr; v.dif

#------------------------------------------------------------------------------------------
# diferença mínima significativa pelo teste de Tukey

delta <- c(qtukey(0.95, nlevels(bib3$trat), anova(mm0)["trat","denDF"])*sqrt(0.5*v.dif))
delta

#------------------------------------------------------------------------------------------
# aplicação do teste de Tukey

agricolae::order.stat(levels(bib3$trat), maju, delta)

#------------------------------------------------------------------------------------------
</code>

----

==== Análise intrabloco e interbloco de experimento em blocos incompletos tipo II ====

Temporariamente sem descrição.\\
palavras-chave: #bloco_incompleto, #intrabloco, #interbloco, #efeito_aleatório, #Tukey, #média_ajustada.

bib2 <- c(35,1, 28,2, 32,2, 37,3, 35,3, 25,4, 28,4, 27,5, 30,5, 32,6, 24,6, 26,7,
          31,7, 27,1, 38,1, 40,3, 36,3, 27,5, 23,5, 30,7, 28,7, 25,2, 26,2, 28,4,
          23,4, 24,6, 28,6, 33,1, 30,1, 22,4, 27,4, 34,7, 32,7, 39,3, 33,3, 24,6,
          28,6, 34,2, 29,2, 26,5, 23,5, 33,1)
bib2 <- matrix(bib2, ncol=2, byrow=TRUE)
bib2 <- as.data.frame(bib2)
bib2$grup <- gl(3,14)
bib2$bloc <- gl(7,2)
names(bib2)[1:2] <- c("resp","trat")
bib2 <- transform(bib2, trat=factor(trat), resp=resp/10)
str(bib2)

m0 <- lm(terms(resp~grup/bloc+trat, keep.order=TRUE), data=bib2,
         contrast=list(grup=contr.sum, bloc=contr.sum, trat=contr.sum))
anova(m0)
car::Anova(m0, type="III")
drop1(m0, scope=.~., test="F")

#------------------------------------------------------------------------------------------
# checagem das pressuposições
 
par(mfrow=c(2,2))
plot(m0)
layout(1)

#------------------------------------------------------------------------------------------
# médias ajustadas

Xtrat <- m0$contrast$trat
assi <- m0$assign

maju <- c(Xtrat%*%coef(m0)[assi==3]+coef(m0)[1])
maju

#------------------------------------------------------------------------------------------
# vetor de coeficientes do contraste tr1 vs tr2 e sua estimativa

contr <- rep(c(0,0,0,1,-1,0), c(1,sum(assi==1),sum(assi==2),1,1,sum(assi==3)-2)); contr
sum(contr*coef(m0))

#------------------------------------------------------------------------------------------
# variância do contraste, como tem equilíbrio, todos os contrastes tem mesma variância

v.dif <- contr%*%vcov(m0)%*%contr; v.dif

#------------------------------------------------------------------------------------------
# diferença mínima significativa pelo teste de Tukey

delta <- qtukey(0.95, nlevels(bib2$trat), df=df.residual(m0))*sqrt(0.5*v.dif); delta

#------------------------------------------------------------------------------------------
# aplicação do teste de Tukey

taju <- agricolae::order.stat(levels(bib2$trat), maju, delta)
taju

#------------------------------------------------------------------------------------------
# gráfico de barras com as médias

taju2 <- taju$means
names(taju2) <- taju$trt

bp <- barplot(taju2,
              xlab="Tratamentos", ylab="Variável resposta",
              ylim=c(0, 1.2*max(taju2)))
text(bp, taju2,
     labels=paste(format(taju$means,dig=2), gsub(" ","",taju$M), sep="\n"), pos=3)
mtext(3, line=2,
      text="Comparação múltipla de médias", cex=1.5)
mtext(3, line=0.5,
      text="Médias seguidas de mesma letra não diferem entre si pelo teste de Tukey (5%)")
box()

#------------------------------------------------------------------------------------------
# análise interbloco
#------------------------------------------------------------------------------------------
# criar o fator bloco dentro de grupo para associar efeito aleatório

bib2$bloc.grup <- factor(paste(bib2$grup, bib2$bloc, sep="-"))
str(bib2)

#------------------------------------------------------------------------------------------
# ajuste do modelo e teste de Wald sequencial para os efeitos fixos

require(nlme)

mm0 <- lme(resp~grup+trat, random=~1|bloc.grup, data=bib2,
           contrast=list(grup=contr.sum, trat=contr.sum))
anova(mm0)

#------------------------------------------------------------------------------------------
# checagem da normalidade efeitos/erros aleatórios

par(mfrow=c(1,2))
qqnorm(residuals(mm0)); qqline(residuals(mm0))
qqnorm(unlist(ranef(mm0))); qqline(unlist(ranef(mm0)))
layout(1)

#------------------------------------------------------------------------------------------
# estimativas dos componentes de variância

VarCorr(mm0)

#------------------------------------------------------------------------------------------
# médias ajustadas

Xtrat <- mm0$contrast$trat
assi <- lapply(list(grup="grup",trat="trat"), function(x){ grep(x, names(fixef(mm0))) })

maju <- c(Xtrat%*%fixef(mm0)[assi$trat]+fixef(mm0)[1])
maju

#------------------------------------------------------------------------------------------
# vetor de coeficientes do contraste tr1 vs tr2 e sua estimativa

contr <- rep(c(0,0,1,-1,0), c(1,length(assi$grup),1,1,length(assi$trat)-2)); contr
sum(contr*fixef(mm0))

#------------------------------------------------------------------------------------------
# variância do contraste, como tem equilíbrio, todos os contrastes tem mesma variância

v.dif <- contr%*%vcov(mm0)%*%contr; v.dif

#------------------------------------------------------------------------------------------
# diferença mínima significativa pelo teste de Tukey

delta <- qtukey(0.95, nlevels(bib2$trat), df=anova(mm0)["trat","denDF"])*sqrt(0.5*v.dif)

#------------------------------------------------------------------------------------------
# aplicação do teste de Tukey

agricolae::order.stat(levels(bib2$trat), maju, delta)

#------------------------------------------------------------------------------------------
</code>

----

==== Análise intrabloco e interbloco de experimento em blocos incompletos tipo I ====

Temporariamente sem descrição.\\
palavras-chave: #bloco_incompleto, #intrabloco, #interbloco, #efeito_aleatório, #Tukey, #média_ajustada.

bib1 <- c(20,1,18,2,15,3,16,4,14,5,15,6,16,7,18,8,
          24,1,18,3,25,2,19,8,13,4,16,5,12,6,16,7,
          23,1,17,4,26,2,18,7,15,3,17,6,13,5,16,8,
          21,1,13,5,23,2,16,3,10,4,12,7,13,6,11,8,
          28,1,14,6,27,2,18,4,18,3,15,8,16,5,17,7,
          22,1,17,7,24,2,16,6,18,3,14,5,15,4,17,8,
          23,1,15,8,21,2,13,5,15,3,12,7,13,4,16,6)
bib1 <- matrix(bib1, ncol=2, byrow=TRUE)
bib1 <- as.data.frame(bib1)
bib1 <- cbind(rept=gl(7,8), bib1)
bib1$bloc <- gl(4,2)
names(bib1) <- c("rept","resp","trat","bloc")
bib1$trat <- factor(bib1$trat)
str(bib1)

#------------------------------------------------------------------------------------------
# análise intrabloco
#------------------------------------------------------------------------------------------
# ajuste do modelo, anova sequencial e marginal, só o teste F para trat é válido

m0 <- lm(terms(resp~rept/bloc+trat, keep.order=TRUE), data=bib1,
         contrast=list(rept=contr.sum, bloc=contr.sum, trat=contr.sum))
anova(m0)
car::Anova(m0, type="III")
drop1(m0, scope=.~., test="F")

#------------------------------------------------------------------------------------------
# checagem das pressuposições

par(mfrow=c(2,2))
plot(m0)
layout(1)

#------------------------------------------------------------------------------------------
# médias ajustadas
 
Xtrat <- m0$contrast$trat
assi <- m0$assign
 
maju <- c(Xtrat%*%coef(m0)[assi==3]+coef(m0)[1])
maju

#------------------------------------------------------------------------------------------
# vetor de coeficientes do contraste tr1 vs tr2 e sua estimativa
 
contr <- rep(c(0,0,0,1,-1,0), c(1,sum(assi==1),sum(assi==2),1,1,sum(assi==3)-2)); contr
sum(contr*coef(m0))
 
#------------------------------------------------------------------------------------------
# variância do contraste, como tem equilíbrio, todos os contrastes tem mesma variância
 
v.dif <- contr%*%vcov(m0)%*%contr; v.dif
 
#------------------------------------------------------------------------------------------
# diferença mínima significativa pelo teste de Tukey
 
delta <- qtukey(0.95, nlevels(bib1$trat), df=df.residual(m0))*sqrt(0.5*v.dif); delta
 
#------------------------------------------------------------------------------------------
# aplicação do teste de Tukey

taju <- agricolae::order.stat(levels(bib1$trat), maju, delta)
taju

#------------------------------------------------------------------------------------------
# gráfico de barras com as médias

taju2 <- taju$means
names(taju2) <- taju$trt

#------------------------------------------------------------------------------------------
# análise interbloco
#------------------------------------------------------------------------------------------
# criar o fator bloco dentro de repetição para associar efeito aleatório

bib1$bloc.rept <- factor(paste(bib1$rept, bib1$bloc, sep="-"))
str(bib1)

#------------------------------------------------------------------------------------------
# ajuste do modelo e teste de Wald sequencial para os efeitos fixos

mm0 <- lme(resp~rept+trat, random=~1|bloc.rept, data=bib1,
           contrast=list(rept=contr.sum, trat=contr.sum))
anova(mm0)

#------------------------------------------------------------------------------------------
# checagem da normalidade efeitos/erros aleatórios

par(mfrow=c(1,2))
qqnorm(residuals(mm0)); qqline(residuals(mm0))
qqnorm(unlist(ranef(mm0))); qqline(unlist(ranef(mm0)))
layout(1)

#------------------------------------------------------------------------------------------
# estimativas dos componentes de variância

VarCorr(mm0)

#------------------------------------------------------------------------------------------
# médias ajustadas
 
Xtrat <- mm0$contrast$trat
assi <- lapply(list(rept="rept",trat="trat"), function(x){ grep(x, names(fixef(mm0))) })
 
maju <- c(Xtrat%*%fixef(mm0)[assi$trat]+fixef(mm0)[1])
maju

#------------------------------------------------------------------------------------------
# vetor de coeficientes do contraste tr1 vs tr2 e sua estimativa
 
contr <- rep(c(0,0,1,-1,0), c(1,length(assi$rept),1,1,length(assi$trat)-2)); contr
sum(contr*fixef(mm0))
 
#------------------------------------------------------------------------------------------
# variância do contraste, como tem equilíbrio, todos os contrastes tem mesma variância
 
v.dif <- contr%*%vcov(mm0)%*%contr; v.dif
 
#------------------------------------------------------------------------------------------
# diferença mínima significativa pelo teste de Tukey
 
delta <- qtukey(0.95, nlevels(bib1$trat), df=anova(mm0)["trat","denDF"])*sqrt(0.5*v.dif)
 
#------------------------------------------------------------------------------------------
# aplicação do teste de Tukey
 
agricolae::order.stat(levels(bib1$trat), maju, delta)
 
#------------------------------------------------------------------------------------------
</code>

----

==== Sobreposição de Gráficos Usando polygon() ====

O exemplo abaixo usa a função polygon() para colorir uma região da densidade normal. \\

palavras-chave: #graficos #sobreposicao #polygon #plot

<code R>
##
## Por Benilton Carvalho
##
## Determinando os parâmetros da Normal
media <- 170
stder <- 8

## Criando o grid no qual deseja-se visualizar a densidade
x <- seq(media-3*stder, media+3*stder, .01)

## Determinando a densidade em cada ponto do grid
y <- dnorm(x, media, stder)

## Gráfico da densidade no grid de interesse
plot(x, y, type='l', xlab='Eixo X', ylab='Eixo Y')
abline(h=0, v=media, lty=2)

## Determinando o sub-grid de interesse
##   e a densidade em cada ponto da regiao
## Região inferior apenas.
hx <- seq(media-3*stder, media-2*stder, .01)
hy <- dnorm(hx, media, stder)

## A funcao polygon() necessita do caminho inteiro
##   da regiao a ser colorida. Então, para cada hx[i],
##   existem 2 pontos determinando o polígono: a) na
##   altura da densidade e b) em zero (para fechar o
##   polígono).
n <- length(hy)
polygon(c(hx, rev(hx)), c(hy, rep(0, n)), col=2)

## Repetindo para a região superior
hx <- seq(media+2*stder, media+3*stder, .01)
hy <- dnorm(hx, media, stder)
n <-
 length(hy)
polygon(c(hx, rev(hx)), c(hy, rep(0, n)), col=2)
</code>

Por favor digite todas as letras dentro da caixa para provar que você é humano. I L᠎ T B A Por favor deixe esse campo em branco:

Resumo da edição

Observe: editando esta página você aceita disponibilizar o seu conteúdo sob a seguinte licença: GNU Free Documentation License 1.3