Não foi possível enviar o arquivo. Será algum problema com as permissões?
Diferenças
Aqui você vê as diferenças entre duas revisões dessa página.
| - | Ambos lados da revisão anterior Revisão anterior Próxima revisão |
| - | ridiculas [2011/02/15 23:26] walmes |
|---|---|
| + | Revisão anterior |
| + | ridiculas [2012/10/25 15:45] (atual) walmes [section 8] |
| @@ Linha -1,9 +1,837 @@ | |
| {{: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''. | |
| ---- | |
| ==== | |
| Temporariamente sem descrição.\\ | |
| palavras-chave: #rpanel. | |
| <code R> | |
| #------------------------------------------------------------------------------------------ | |
| # 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) | |
| #------------------------------------------------------------------------------------------ | |
| </code> | |
| ---- | |
| ==== .Rprofile no Linux ==== | |
| <code R> | |
| #------------------------------------------------------------------------------------------ | |
| # por JCFaria | |
| #------------------------------------------------------------------------------------------ | |
| </code> | |
| Esse post tem a finalidade de compartilhar algumas coisas que considero importantes na inicialização do R no Linux! | |
| Muitas das opções importantes, do ponto de vista funcional (não relativos à aparência), podem ser feitas no arquivo .Rprofile. | |
| Esse arquivo deve ficar localizado no home do usuário (~/.Rprofile) e é um dos primeiros a ser lido quando uma sessão do | |
| R é iniciada. | |
| Tenho uma função (bem simples) que uso bastante em meu dia a dia: "cv" para calcular o coef. de variação de uma ANOVA: | |
| <code R> | |
| 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)) | |
| } | |
| </code> | |
| Pois bem, a medida que vamos aumentando nossa intimidade com o R (inevitavelmente) iremos desenvolvendo nossas próprias | |
| funções (o R foi projetado para isso). | |
| <fc #000080>Ai vem o problema: ter que sempre carregar a função quando for usar, o que pode se tornar uma chatisse! | |
| Pior ainda, ao limpar o workspace do usuário (.GlobalEnv) ela é removida e precisa ser recarregada. Imagine o contexto da necessidade de várias funções, que é o que geralmente acontece com o passar do tempo ...</fc> | |
| <fc #000080>**Tem como contornar? Sim! De várias formas:**</fc> | |
| **//1. Opção muito pouco prática://** | |
| <code R> | |
| oldp <- getwd() | |
| setwd('/home/jcfaria/dados/r/funcoes/') | |
| source('cv.r') | |
| setwd(oldp) | |
| </code> | |
| A função "cv" ficará disponível no meu workspace mas será removida com a instrução: | |
| <code R> | |
| > rm(list=ls()) | |
| </code> | |
| 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: | |
| <code R> | |
| > rm(list=ls()) | |
| </code> | |
| Contudo, poderá ser removida com a intrução: | |
| <code R> | |
| > rm(list=ls(all=TRUE)) | |
| </code> | |
| **//3. Colocando suas funções em algum ambiente (environment) do R (optei pelo base)://** | |
| <code R> | |
| oldp <- getwd() | |
| setwd('/home/jcfaria/dados/r/funcoes/') | |
| source('cv.r', local=baseenv()) | |
| setwd(oldp) | |
| </code> | |
| Ela não ficará no meu workspace, mas sim no base. | |
| Como tal, poderá ser usada com qualquer outra função desse pacote. | |
| **//4. Criando seu próprio ambiente (acho a solução mais elegante)://** | |
| <code R> | |
| oldp <- getwd() | |
| setwd('/home/jcfaria/dados/r/funcoes/') | |
| .jcf <- new.env() | |
| source('cv.r', local=.jcf) | |
| setwd(oldp) | |
| </code> | |
| Nesse último caso [**//4//**]: | |
| * O objeto ".jcf" ficará oculto no meu workspace evitando ser deletado com: rm(list=ls()) | |
| * Parar acessar a função "cv" é necessário: | |
| <code R> | |
| > .jcf$cv | |
| # ou | |
| > with(.jcf, cv) | |
| </code> | |
| Por exemplo: | |
| <code R> | |
| > av <- aov(Sepal.Length ~ Species, data=iris) | |
| > .jcf$cv(av) | |
| [1] 8.81 | |
| # ou | |
| > with(.jcf, cv(av)) | |
| [1] 8.81 | |
| </code> | |
| Esta forma de carregar funções de forma permanente no R pode ser usado para qualquer outro objeto! | |
| No Windows bastava usar no /etc/Rprofile.site: | |
| <code R> | |
| source('cv', local=TRUE) | |
| </code> | |
| que ela ficava disponível no pacote base. | |
| No Windows não testei na versão em desenvolvimento (instável) que uso no Linux, | |
| mas deve funcionar, pois usei por muitos anos em várias versões. | |
| Abaixo meu .Rprofile: | |
| <code R> | |
| ## José Cláudio Faria/UESC/DCET | |
| ##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬ | |
| ## General options | |
| library(utils) | |
| options(list(repos='http://cran-r.c3sl.ufpr.br/', | |
| width=80, | |
| editor='vim', | |
| browser='chromium')) | |
| ##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬ | |
| ## Package debug: options | |
| options(list(debug.height=10, | |
| debug.width=125, | |
| debug.screen.pos='-1680+00' # dois monitores! | |
| )) | |
| ##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬ | |
| ## Funtions: myself | |
| oldp <- getwd() | |
| setwd('/home/jcfaria/dados/r/funcoes/') | |
| #source('fitreg.r', local=baseenv()) | |
| #source('fitregl.r', local=baseenv()) | |
| #source('fitrlcor.r', local=baseenv()) | |
| #source('plotreg.r', local=baseenv()) | |
| #source('plotrl.r', local=baseenv()) | |
| #source('plotlm3d.r', local=baseenv()) | |
| #source('cv.r', local=baseenv()) | |
| #source('rs.r', local=baseenv()) | |
| #source('sourcedir.r', local=baseenv()) | |
| .jcf <- new.env() | |
| source('fitreg.r', local=.jcf) | |
| source('fitregl.r', local=.jcf) | |
| source('fitrlcor.r', local=.jcf) | |
| source('plotreg.r', local=.jcf) | |
| source('plotrl.r', local=.jcf) | |
| source('plotlm3d.r', local=.jcf) | |
| source('cv.r', local=.jcf) | |
| source('rs.r', local=.jcf) | |
| source('sourcedir.r', local=.jcf) | |
| #attach(.jcf) | |
| ##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬ | |
| ## Removing variables | |
| setwd(oldp) | |
| rm(oldp) | |
| ##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬ | |
| ## Packages: loading | |
| #library(fdth) | |
| #library(ScottKnott) | |
| #library(bpca) | |
| #library(TukeyC) | |
| #library(reshape) | |
| #library(debug) | |
| </code> | |
| ---- | |
| ==== Como fazer a justaposição de vários data.frames ==== | |
| Temporariamente sem descrição.\\ | |
| palavras-chave: #merge, #Reduce. | |
| <code R> | |
| #------------------------------------------------------------------------------------------ | |
| # 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() | |
| #------------------------------------------------------------------------------------------ | |
| </code> | |
| ---- | |
| ==== Gráfico de valores observados e curva de valores preditos para modelo linear generalizado ==== | |
| Temporariamente sem descrição.\\ | |
| palavras-chave: #glm, #poisson, #predict. | |
| <code R> | |
| #------------------------------------------------------------------------------------------ | |
| # 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)) | |
| #------------------------------------------------------------------------------------------ | |
| </code> | |
| ---- | |
| ==== Gráficos de barras com intervalos de confiança para as médias ==== | |
| Temporariamente sem descrição.\\ | |
| palavras-chave: #intervalo, #erro_padrão, #barras. | |
| <code R> | |
| #------------------------------------------------------------------------------------------ | |
| # 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") | |
| #------------------------------------------------------------------------------------------ | |
| </code> | |
| ---- | |
| ==== Gráficos do R exportados pelo dispositivo tikz ==== | |
| Temporariamente sem descrição.\\ | |
| palavras-chave: #tikz, #sweave, #latex. | |
| <code> | |
| \documentclass{article} | |
| \usepackage{Sweave} | |
| \usepackage{tikz} | |
| \SweaveOpts{keep.source=true} | |
| \title{Usando ti\textit{k}z no Sweave} | |
| \author{Walmes Zeviani\\ LEG/UFPR} | |
| \begin{document} | |
| \maketitle | |
| C\'{o}digo m\'{i}nimo reproduz\'{i}vel usando \texttt{tikzDevice} para exportar gr\'{a}ficos feitos no R | |
| para codifica\c{c}\~{a}o ti\textit{k}z. | |
| {\footnotesize | |
| <<results=hide>>= | |
| #------------------------------------------------------------------------------------------ | |
| require(tikzDevice) | |
| set.seed(2011); x <- rnorm(200) | |
| tikz("plot.tex", w=5, h=3) | |
| hist(x, freq=FALSE, ylab="Densidade", | |
| main="Histograma de uma amostra de $X \\sim N(\\mu=0, \\sigma^2=1)$") | |
| curve(dnorm(x), col=2, add=TRUE, lwd=2); rug(x) | |
| legend("topleft", col=2, lty=1, lwd=2, bty="n", | |
| legend="$\\displaystyle \\frac{1}{\\sqrt{2\\pi\\sigma^2}}\\cdot e^{-\\frac{(x-\\mu)^2}{2\\sigma^2}}$") | |
| box(); dev.off() | |
| #------------------------------------------------------------------------------------------ | |
| @ | |
| } | |
| \input{plot.tex} | |
| \end{document} | |
| </code> | |
| ---- | |
| ==== Gráfico com dois eixos coordenados ==== | |
| Temporariamente sem descrição.\\ | |
| palavras-chave: #série, #eixo. | |
| <code R> | |
| #------------------------------------------------------------------------------------------ | |
| # por Ivan, Benilton e Walmes | |
| #------------------------------------------------------------------------------------------ | |
| # dados de séries de dados indexadas no tempo (meses) | |
| lines <- 'meses temp umidade rad chuva | |
| Jan 26.49 86.58 795.88 0.36 | |
| Fev 26.65 88.49 710.24 0.34 | |
| Mar 27.19 86.16 772.99 0.21 | |
| Abr 26.28 89.75 574.88 0.67 | |
| Mai 26.62 89.22 614.02 0.31 | |
| Jun 26.13 87.83 680.08 0.26 | |
| Jul 25.83 86.57 675.97 0.15 | |
| Ago 27.05 83.14 756.44 0.07 | |
| Set 27.60 83.02 925.57 0.14 | |
| Out 27.44 85.16 927.71 0.17 | |
| Nov 26.56 88.18 788.87 0.19 | |
| Dez 25.87 90.63 703.94 0.33' | |
| da <- read.table(textConnection(lines), header=TRUE) | |
| str(da) | |
| #------------------------------------------------------------------------------------------ | |
| # 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) | |
| #------------------------------------------------------------------------------------------ | |
| </code> | |
| ---- | |
| ==== Análise de dados de proporção usando modelo linear generalizado ==== | |
| Temporariamente sem descrição.\\ | |
| palavras-chave: #binomial, #sucessos, #deviance, #wireframe, #fatorial, #superfície. | |
| <code R> | |
| #------------------------------------------------------------------------------------------ | |
| # 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) | |
| #------------------------------------------------------------------------------------------ | |
| </code> | |
| ---- | |
| ==== Desdobramento de interação em experimento fatorial ==== | |
| Temporariamente sem descrição.\\ | |
| palavras-chave: #desdobramento, #comparações, #teste_de_médias, #fatorial. | |
| <code R> | |
| #------------------------------------------------------------------------------------------ | |
| # 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) | |
| }) | |
| #------------------------------------------------------------------------------------------ | |
| </code> | |
| ---- | |
| ==== Gráfico para comportamento do desempenho em avaliações ==== | |
| Temporariamente sem descrição.\\ | |
| palavras-chave: . | |
| <code R> | |
| #------------------------------------------------------------------------------------------ | |
| # 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() | |
| #------------------------------------------------------------------------------------------ | |
| </code> | |
| ---- | |
| ==== Desdobramento de interação usando a multcomp::glht() ==== | |
| Temporariamente sem descrição.\\ | |
| palavras-chave: #desdobramento, #glht, #comparações, #teste_de_médias, #fatorial, #contrastes. | |
| <code R> | |
| #------------------------------------------------------------------------------------------ | |
| # 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 | |
| #------------------------------------------------------------------------------------------ | |
| </code> | |
| ---- | |
| ==== Regressão na análise de variância ==== | |
| No exemplo abaixo fazemos a análise de um experimento fatorial. São combinados três níveis de um fator qualitativo (''cultivar'' de sorgo) com 6 níveis de um fator quantitativo (''dose'' de nitrogênio). A resposta avaliada nas parcelas é o índice agronômico. O experimento foi instalado em blocos. A interação entre as fontes de variação foi significativa. O desdobramento dessa interação foi feita nos dois sentidos, ou seja, estudo do fator dose fixando nível do fator cultivar e estudo do fator cultivar fixando o nível do fator dose. Foi feito a ajuste de modelos de regressão polinomial para o estudo do fator dose. São ilustrados passos para teste de hipótese sobre os efeitos linear, quadrático, etc, e falta de ajuste do modelo ajustado, além de como obter as estimativas dos termos do modelo e o coeficiente de determinação (R²).\\ | |
| palavras-chave: #regressão, #polinômio, #anova. | |
| @@ Linha -116,7 +944,9 @@ | |
| </code> | |
| 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. | |
| @@ Linha -173,5 +1003,7 @@ | |
| </code> | |
| ==== Experimento em parcelas subdivididas com resposta do tipo binomial | |
| Temporariamente sem descrição.\\ | |
| @@ Linha -221,5 +1053,7 @@ | |
| </code> | |
| ---- | |
| ==== | |
| Temporariamente sem descrição.\\ | |
| @@ Linha -239,2 +1073,595 @@ | |
| #------------------------------------------------------------------------------------------ | |
| </code> | |
| ---- | |
| ==== Representação de ajuste de curvas de regressão no mesmo gráfico ==== | |
| Faremos o ajuste de três curvas de regressão à dados de índice agronômico de sorgo em função da dose de nitrogênio aplicada. Os dados estão classificados para três cultivares de sorgo. Apesar de se tratar de um experimento em blocos, isso será ignorado no exemplo que tem o objetivo de apresentar os comandos para ajustar curvas e representar o ajuste no mesmo gráfico. | |
| palavras-chave: #regressão, #lapply, #curvas, #ajuste, #lines, #subset. | |
| <code R> | |
| #------------------------------------------------------------------------------------------ | |
| # por Walmes | |
| #------------------------------------------------------------------------------------------ | |
| # dados. | |
| sorgo <- read.table("http://www.leg.ufpr.br/~walmes/docs/anovareg.txt", header=TRUE) | |
| sorgo <- transform(sorgo, bloco=factor(bloco), cultivar=factor(cultivar)) | |
| str(sorgo) | |
| #------------------------------------------------------------------------------------------ | |
| # usaremos esse banco de dados para ajustarmos curvas de regressão separadas por cultivar. | |
| # iremos ignorar a presença de blocos assumindo amostragem aleatória simples dos dados. | |
| #------------------------------------------------------------------------------------------ | |
| # gráfico de dispersão dos dados. | |
| levels(sorgo$cultivar) | |
| par(mfrow=c(1,3)) | |
| with(sorgo, sapply(levels(sorgo$cultivar), | |
| function(x){ | |
| plot(indice~dose, subset(sorgo, cultivar==x), main=x) | |
| } | |
| ) | |
| ) | |
| #------------------------------------------------------------------------------------------ | |
| # ajustaremos uma curva para cada cultivar. pode-se usar a lapply o tempo todo para isso. | |
| levels(sorgo$cultivar) | |
| mAG <- lm(indice~dose+I(dose^2), data=subset(sorgo, cultivar==levels(sorgo$cultivar)[1])) | |
| mBR <- lm(indice~dose+I(dose^2), data=subset(sorgo, cultivar==levels(sorgo$cultivar)[2])) | |
| mPi <- lm(indice~dose+I(dose^2), data=subset(sorgo, cultivar==levels(sorgo$cultivar)[3])) | |
| lapply(list(mAG, mBR, mPi), summary) | |
| lapply(list(mAG, mBR, mPi), anova) | |
| #------------------------------------------------------------------------------------------ | |
| # agora faremos os gráficos dos valores observados com as retas dentro | |
| # jitter() é para "agitar" os valores de dose antes de plotar, apenas para melhorar aspecto | |
| layout(1) | |
| plot(indice~jitter(dose), data=sorgo, col=sorgo$cultivar) | |
| with(subset(sorgo, cultivar==levels(sorgo$cultivar)[1]), lines(predict(mAG)~dose, col=1)) | |
| with(subset(sorgo, cultivar==levels(sorgo$cultivar)[2]), lines(predict(mBR)~dose, col=2)) | |
| with(subset(sorgo, cultivar==levels(sorgo$cultivar)[3]), lines(predict(mPi)~dose, col=3)) | |
| #------------------------------------------------------------------------------------------ | |
| </code> | |
| Pode-se usar as funções ''sapply()'' e/ou ''lapply()'' para tornar o código mais eficiente/enxuto. | |
| ---- | |
| ==== Análise intrabloco e interbloco de experimento em blocos incompletos tipo III ==== | |
| Para análise intrabloco pode-se usar a função ''agricolae::BIB.test()''.\\ | |
| palavras-chave: #bloco_incompleto, #intrabloco, #interbloco, #efeito_aleatório, #Tukey, #média_ajustada. | |
| <code R> | |
| #------------------------------------------------------------------------------------------ | |
| # 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) | |
| #------------------------------------------------------------------------------------------ | |
| </code> | |
| ---- | |
| ==== Análise intrabloco e interbloco de experimento em blocos incompletos tipo II ==== | |
| Temporariamente sem descrição.\\ | |
| palavras-chave: #bloco_incompleto, #intrabloco, #interbloco, #efeito_aleatório, #Tukey, #média_ajustada. | |
| <code R> | |
| #------------------------------------------------------------------------------------------ | |
| # 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) | |
| #------------------------------------------------------------------------------------------ | |
| </code> | |
| ---- | |
| ==== Análise intrabloco e interbloco de experimento em blocos incompletos tipo I ==== | |
| Temporariamente sem descrição.\\ | |
| palavras-chave: #bloco_incompleto, #intrabloco, #interbloco, #efeito_aleatório, #Tukey, #média_ajustada. | |
| <code R> | |
| #------------------------------------------------------------------------------------------ | |
| # 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) | |
| #------------------------------------------------------------------------------------------ | |
| </code> | |
| ---- | |
| ==== Sobreposição de Gráficos Usando polygon() ==== | |
| O exemplo abaixo usa a função polygon() para colorir uma região da densidade normal. \\ | |
| palavras-chave: #graficos #sobreposicao #polygon #plot | |
| <code R> | |
| ## | |
| ## Por Benilton Carvalho | |
| ## | |
| ## Determinando os parâmetros da Normal | |
| media <- 170 | |
| stder <- 8 | |
| ## Criando o grid no qual deseja-se visualizar a densidade | |
| x <- seq(media-3*stder, media+3*stder, .01) | |
| ## Determinando a densidade em cada ponto do grid | |
| y <- dnorm(x, media, stder) | |
| ## Gráfico da densidade no grid de interesse | |
| plot(x, y, type='l', xlab='Eixo X', ylab='Eixo Y') | |
| abline(h=0, v=media, lty=2) | |
| ## Determinando o sub-grid de interesse | |
| ## e a densidade em cada ponto da regiao | |
| ## Região inferior apenas. | |
| hx <- seq(media-3*stder, media-2*stder, .01) | |
| hy <- dnorm(hx, media, stder) | |
| ## A funcao polygon() necessita do caminho inteiro | |
| ## da regiao a ser colorida. Então, para cada hx[i], | |
| ## existem 2 pontos determinando o polígono: a) na | |
| ## altura da densidade e b) em zero (para fechar o | |
| ## polígono). | |
| n <- length(hy) | |
| polygon(c(hx, rev(hx)), c(hy, rep(0, n)), col=2) | |
| ## Repetindo para a região superior | |
| hx <- seq(media+2*stder, media+3*stder, .01) | |
| hy <- dnorm(hx, media, stder) | |
| n <- | |
| length(hy) | |
| polygon(c(hx, rev(hx)), c(hy, rep(0, n)), col=2) | |
| </code> | |
