===== R-idículas - dicas curtas sobre R ===== {{:nuvem.png?550 |}} **//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 <>= #------------------------------------------------------------------------------------------ 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) #------------------------------------------------------------------------------------------ # para enviar para outra pessoa pode-se enviar o resultado do comando dput(), assim dput(da) # cole na mensagem o resulado que aparece no console ao executar essa linha # assim a pessoa pode ler fazendo assim da <- structure(list(meses = structure(c(5L, 4L, 9L, 1L, 8L, 7L, 6L, 2L, 12L, 11L, 10L, 3L), .Label = c("Abr", "Ago", "Dez", "Fev", "Jan", "Jul", "Jun", "Mai", "Mar", "Nov", "Out", "Set"), class = "factor"), temp = c(26.49, 26.65, 27.19, 26.28, 26.62, 26.13, 25.83, 27.05, 27.6, 27.44, 26.56, 25.87), umidade = c(86.58, 88.49, 86.16, 89.75, 89.22, 87.83, 86.57, 83.14, 83.02, 85.16, 88.18, 90.63), rad = c(795.88, 710.24, 772.99, 574.88, 614.02, 680.08, 675.97, 756.44, 925.57, 927.71, 788.87, 703.94), chuva = c(0.36, 0.34, 0.21, 0.67, 0.31, 0.26, 0.15, 0.07, 0.14, 0.17, 0.19, 0.33)), .Names = c("meses", "temp", "umidade", "rad", "chuva" ), class = "data.frame", row.names = c(NA, -12L)) #------------------------------------------------------------------------------------------ # 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::glht() ==== 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 ==== 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 I (Pimentel Gomes p.185) 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 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 #------------------------------------------------------------------------------------------ # 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) #------------------------------------------------------------------------------------------ ---- ==== 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 ## ## 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)