ridiculas - Wiki do LEG

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

Você está aqui: start » ridiculas

Diferenças

Aqui você vê as diferenças entre duas revisões dessa página.

-	ridiculas [2011/03/24 19:23] walmes [section 7]
-	Ambos lados da revisão anterior Revisão anterior Próxima revisão
+	Revisão anterior
+	ridiculas [2012/10/25 15:45] (atual) walmes [section 8]
@@ Linha -1,8 +1,835 @@ ~~removida~~ criada

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

	//Ridículas// é a página do LEG dedicada à fornecer //dicas curtas// sobre R, e.g. 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 Ridículas, envie e-mail para ''walmes@ufpr.br''{{: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.

	<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 AK e KA é 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²).\\
@@ Linha -117,5 +944,7 @@ ~~removida~~ criada
	</code>

	=----

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

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

	=----

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

	Temporariamente sem descrição.\\
@@ Linha -222,5 +1053,7 @@ ~~removida~~ criada
	</code>

	=----

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

	Temporariamente sem descrição.\\
@@ Linha -241,6 +1074,7 @@ ~~removida~~ criada
	</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.
@@ Linha -298,8 +1132,9 @@ ~~removida~~ criada
	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 =====

	~~Temporariamente sem descrição~~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.

@@ Linha -332,11 +1167,18 @@ ~~removida~~ criada
	# soma de quadrados sequencial com tratamentos ajustados aos blocos (intrabloco)

	m0 <- lm(resp~bloc+trat, data=~~tipo3~~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=~~tipo3~~bib3)
	anova(m1)

@@ Linha -345,4 +1187,5 @@ ~~removida~~ criada

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

	#------------------------------------------------------------------------------------------
@@ Linha -386,5 +1229,22 @@ ~~removida~~ criada
	# 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()

	#------------------------------------------------------------------------------------------
@@ Linha -394,6 +1254,15 @@ ~~removida~~ criada

	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)

	#------------------------------------------------------------------------------------------
@@ Linha -438,3 +1307,361 @@ ~~removida~~ criada
	#------------------------------------------------------------------------------------------
	</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.5v.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.5v.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.5v.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.5v.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-3stder, media+3stder, .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-3stder, media-2stder, .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+2stder, media+3stder, .01)
	hy <- dnorm(hx, media, stder)
	n <-
	length(hy)
	polygon(c(hx, rev(hx)), c(hy, rep(0, n)), col=2)
	</code>

Diferenças

Navegação

Busca

Ferramentas

QR Code