ridiculas - Wiki do LEG

Não foi possível enviar o arquivo. Será algum problema com as permissões?

Você está aqui: start » ridiculas

Essa é uma revisão anterior do documento!

Tabela de conteúdos

Controlando gráficamente parâmetros de distribuições de probabilidade
.Rprofile no Linux
Como fazer a justaposição de vários data.frames
Gráfico de valores observados e curva de valores preditos para modelo linear generalizado
Gráficos de barras com intervalos de confiança para as médias
Gráficos do R exportados pelo dispositivo tikz
Gráfico com dois eixos coordenados
Análise de dados de proporção usando modelo linear generalizado
Desdobramento de interação em experimento fatorial
Gráfico para comportamento do desempenho em avaliações
Desdobramento de interação usando a multcomp::ghlt()
Regressão na análise de variância
Experimento com dois fatores de efeito aditivo e perda de muitas parcelas
Experimento em parcelas subdivididas com resposta do tipo binomial
Duas ou mais funções em um mesmo gráfico
Representação de ajuste de curvas de regressão no mesmo gráfico
Análise intrabloco e interbloco de experimento em blocos incompletos tipo III
Análise intrabloco e interbloco de experimento em blocos incompletos tipo II
Análise intrabloco e interbloco de experimento em blocos incompletos tipo I
Sobreposição de Gráficos Usando polygon()

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.

#------------------------------------------------------------------------------------------
#                                                                                por Walmes
#------------------------------------------------------------------------------------------
 
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)
 
#------------------------------------------------------------------------------------------

.Rprofile no Linux

#------------------------------------------------------------------------------------------
#                                                                               por JCFaria
#------------------------------------------------------------------------------------------

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:

cv <- function(av)
{
  if(is.null(av) || !inherits(av, 'aov'))
    stop('Please, check the parameter!')
  qmee <- with(av, sum(residuals^2) / df.residual)
  cv   <- 100 * sqrt(qmee) / mean(av$fitted.values)
  return(round(cv, 2))
}

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

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 …

Tem como contornar? Sim! De várias formas:

1. Opção muito pouco prática:

oldp <- getwd()
setwd('/home/jcfaria/dados/r/funcoes/') 
source('cv.r')
setwd(oldp)

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

> rm(list=ls())

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:

> rm(list=ls())

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

> rm(list=ls(all=TRUE))

3. Colocando suas funções em algum ambiente (environment) do R (optei pelo base):

oldp <- getwd()
setwd('/home/jcfaria/dados/r/funcoes/') 
source('cv.r', local=baseenv())
setwd(oldp)

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

oldp <- getwd()
setwd('/home/jcfaria/dados/r/funcoes/') 
.jcf <- new.env()
source('cv.r', local=.jcf)
setwd(oldp)

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:

> .jcf$cv
 
# ou
 
> with(.jcf, cv)

Por exemplo:

> av <- aov(Sepal.Length ~ Species, data=iris)
> .jcf$cv(av)
[1] 8.81
 
# ou
 
> with(.jcf, cv(av))
[1] 8.81

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:

source('cv', local=TRUE)

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:

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

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

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

#------------------------------------------------------------------------------------------
#                                                                                por Walmes
#------------------------------------------------------------------------------------------
 
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()
#------------------------------------------------------------------------------------------

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

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

#------------------------------------------------------------------------------------------
#                                                                                por Walmes
#------------------------------------------------------------------------------------------
 
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))
 
#------------------------------------------------------------------------------------------

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

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

#------------------------------------------------------------------------------------------
#                                                                                por Walmes
#------------------------------------------------------------------------------------------
# cmr para colocar intervalos de confiança para as medias num gráfico de 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")
 
#------------------------------------------------------------------------------------------

Gráficos do R exportados pelo dispositivo tikz

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

\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}

Gráfico com dois eixos coordenados

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

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

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

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

#------------------------------------------------------------------------------------------
#                                                                                por Walmes
#------------------------------------------------------------------------------------------
# dados de número de sementes viáveis de soja
 
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)
 
#------------------------------------------------------------------------------------------

Desdobramento de interação em experimento fatorial

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

#------------------------------------------------------------------------------------------
#                                                                                por Walmes
#------------------------------------------------------------------------------------------
# importa e prepara os dados
 
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)
         })
 
#------------------------------------------------------------------------------------------

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

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

#------------------------------------------------------------------------------------------
#                                                                                por Walmes
#------------------------------------------------------------------------------------------
# gráficos para apresentar comportamento de notas de alunos em provas
 
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()
 
#------------------------------------------------------------------------------------------

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

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

#------------------------------------------------------------------------------------------
#                                                                                por Walmes
#------------------------------------------------------------------------------------------
# dados e ajuste
 
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
 
#------------------------------------------------------------------------------------------

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.

#------------------------------------------------------------------------------------------
#                                                                                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
                })})
                                                                                          #
#------------------------------------------------------------------------------------------

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.

#------------------------------------------------------------------------------------------
#                                                                                por Walmes
#------------------------------------------------------------------------------------------
# dados
 
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")))
 
#------------------------------------------------------------------------------------------

Experimento em parcelas subdivididas com resposta do tipo binomial

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

#------------------------------------------------------------------------------------------
#                                                                                por Walmes
#------------------------------------------------------------------------------------------
# pacote que contem a função glmmPQL. pode-se usar a lme4::glmer()
 
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)
 
#------------------------------------------------------------------------------------------

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

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

#------------------------------------------------------------------------------------------
#                                                                                por Walmes
#------------------------------------------------------------------------------------------
# cria os vetores usando as funções abaixo e plota 2 num mesmo gráfico
 
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")
 
#------------------------------------------------------------------------------------------

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.

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

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.

#------------------------------------------------------------------------------------------
#                                                                                por Walmes
#------------------------------------------------------------------------------------------
# dados de um experimento em blocos incompletos tipo III (Pimentel Gomes, p.179)
 
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)
 
#------------------------------------------------------------------------------------------

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.

#------------------------------------------------------------------------------------------
#                                                                                por Walmes
#------------------------------------------------------------------------------------------
# dados de um experimento em blocos incompletos tipo II (Pimentel Gomes p.188)
 
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)
 
#------------------------------------------------------------------------------------------
# análise intrabloco
#------------------------------------------------------------------------------------------
# ajuste do modelo, anova sequencial e marginal
 
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)
 
#------------------------------------------------------------------------------------------

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