##-----------------------------------------------------------------------------
## Definições da sessão, pacotes a serem usados.
pkg <- c("lattice", "latticeExtra", "doBy", "multcomp", "reshape",
"plyr", "wzRfun")
sapply(pkg, library, character.only=TRUE, logical.return=TRUE)
## lattice latticeExtra doBy multcomp reshape plyr
## TRUE TRUE TRUE TRUE TRUE TRUE
## wzRfun
## TRUE
source("lattice_setup.R")
##-----------------------------------------------------------------------------
## Informações sobre as versões dos pacotes.
sessionInfo()
## R version 3.1.1 (2014-07-10)
## Platform: i686-pc-linux-gnu (32-bit)
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=pt_BR.UTF-8
## [4] LC_COLLATE=en_US.UTF-8 LC_MONETARY=pt_BR.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=pt_BR.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] splines stats graphics grDevices utils datasets base
##
## other attached packages:
## [1] wzRfun_0.1 plyr_1.8.1 reshape_0.8.5 multcomp_1.3-3
## [5] TH.data_1.0-3 mvtnorm_0.9-99992 doBy_4.5-10 MASS_7.3-33
## [9] survival_2.37-7 latticeExtra_0.6-26 RColorBrewer_1.0-5 lattice_0.20-29
## [13] knitr_1.5
##
## loaded via a namespace (and not attached):
## [1] evaluate_0.5.5 formatR_0.10 grid_3.1.1 lme4_1.1-6
## [5] Matrix_1.1-4 methods_3.1.1 minqa_1.2.3 nlme_3.1-117
## [9] Rcpp_0.11.1 RcppEigen_0.3.2.1.2 sandwich_2.3-0 stringr_0.6.2
## [13] tools_3.1.1 zoo_1.7-11
## obs: Para instalar um pacote faça:
## install.packages("nome_do_pacote", dependencies=TRUE)
##-----------------------------------------------------------------------------
## Lendo arquivos de dados.
## Diretório de trabalho.
getwd()
## [1] "/home/walmes/Dropbox/CursoR/CIAEEAR/scripts"
## Arquivos no diretório.
list.files()
## [1] "as.lm.R" "ciaeear_01dic.html" "ciaeear_01dic.R" "ciaeear_01dic.Rmd"
## [5] "ciaeear_02dic.R" "ciaeear_02dic.Rmd" "ciaeear_03dbc.R" "ciaeear_03dbc.Rmd"
## [9] "ciaeear_03.R" "ciaeear_04dbc.R" "ciaeear_04dbc.Rmd" "ciaeear_05dql.R"
## [13] "ciaeear_05dql.Rmd" "ciaeear_06dql.R" "ciaeear_06dql.Rmd" "ciaeear_07dic.R"
## [17] "ciaeear_07dic.Rmd" "ciaeear_08fat.R" "ciaeear_08fat.Rmd" "ciaeear_09anc.R"
## [21] "ciaeear_09anc.Rmd" "ciaeear_10fat.R" "ciaeear_10fat.Rmd" "ciaeear_11reg.R"
## [25] "ciaeear_11reg.Rmd" "ciaeear_12reg.R" "ciaeear_12reg.Rmd" "ciaeear_13reg.R"
## [29] "ciaeear_13reg.Rmd" "ciaeear_14nls.Rmd" "ciaeear_init.Rmd" "fig"
## [33] "knitr_setup.R" "lattice_setup.R" "ordinal.Rmd"
## Url de um arquivo com dados.
url <- "http://www.leg.ufpr.br/~walmes/data/pimentel_racoes.txt"
## Opções:
## 1. Fazer download do arquivo para sua máquina em seguida ler.
## 2. Fazer a leitura do arquivo direto do link, sem ter cópia na sua
## máquina.
## Faz o download do arquivo a partir do link.
## download.file(url, dest=basename(url))
## Arquivos de extensão ".txt".
## list.files(pattern="\\.txt$")
## Importa os dados a partir do arquivo local.
## da <- read.table(file="pimentel_racoes.txt", header=TRUE, sep="\t")
## da
## Importa os dados a partir do arquivo na web.
da <- read.table(file=url, header=TRUE, sep="\t")
## Mostra a estrutura do objeto.
str(da)
## 'data.frame': 20 obs. of 2 variables:
## $ racoes : Factor w/ 4 levels "A","B","C","D": 1 1 1 1 1 2 2 2 2 2 ...
## $ ganhopeso: int 35 19 31 15 30 40 35 46 41 33 ...
## Os dados no formato amplo (não empilhado).
unstack(da, form=ganhopeso~racoes)
## A B C D
## 1 35 40 39 27
## 2 19 35 27 12
## 3 31 46 20 13
## 4 15 41 29 28
## 5 30 33 45 30
## Os dados são de um experimento para verificar o ganho de peso em
## função de tipos de rações. Foram usandos 4 níveis do fator ração e 5
## repetições.
##-----------------------------------------------------------------------------
## Análise exploratória.
xyplot(ganhopeso~racoes, data=da, type=c("p","a"))
\[ \begin{aligned} Y|\text{fontes de variação} &\sim \text{Normal}(\mu_i, \sigma^2) \newline \mu_i &= \mu+\alpha_i \end{aligned} \]
em que \(\alpha_i\) é o efeito da ração \(i\) sobre a variável resposta \(Y\), \(\mu\) é a média de \(Y\) na ausência do efeito das rações, \(\mu_i\) é a média para as observações em cada nível de ração e \(\sigma^2\) é a variância das observações ao redor desse valor médio.
##-----------------------------------------------------------------------------
## Estimação.
m0 <- lm(ganhopeso~racoes, data=da)
## Estimativas dos efeitos.
coef(m0)
## (Intercept) racoesB racoesC racoesD
## 26 13 6 -4
## Estimativas e erros padrões.
summary(m0)
##
## Call:
## lm(formula = ganhopeso ~ racoes, data = da)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.00 -6.25 1.50 6.25 13.00
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 26.00 3.71 7.01 2.9e-06 ***
## racoesB 13.00 5.24 2.48 0.025 *
## racoesC 6.00 5.24 1.14 0.269
## racoesD -4.00 5.24 -0.76 0.457
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.29 on 16 degrees of freedom
## Multiple R-squared: 0.428, Adjusted R-squared: 0.321
## F-statistic: 3.99 on 3 and 16 DF, p-value: 0.0267
## Restrição paramétrica de zerer o primeiro nível.
contrasts(da$racoes)
## B C D
## A 0 0 0
## B 1 0 0
## C 0 1 0
## D 0 0 1
model.matrix(m0)
## (Intercept) racoesB racoesC racoesD
## 1 1 0 0 0
## 2 1 0 0 0
## 3 1 0 0 0
## 4 1 0 0 0
## 5 1 0 0 0
## 6 1 1 0 0
## 7 1 1 0 0
## 8 1 1 0 0
## 9 1 1 0 0
## 10 1 1 0 0
## 11 1 0 1 0
## 12 1 0 1 0
## 13 1 0 1 0
## 14 1 0 1 0
## 15 1 0 1 0
## 16 1 0 0 1
## 17 1 0 0 1
## 18 1 0 0 1
## 19 1 0 0 1
## 20 1 0 0 1
## attr(,"assign")
## [1] 0 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$racoes
## [1] "contr.treatment"
##-----------------------------------------------------------------------------
## Existe efeito de rações? Ou seja, todos os \alpha_i podem ser
## considerados iguais à zero?
## A hipótese H0: \alpha_i==0 para todo i pode ser avaliada de duas
## formas: 1) pelo valor de F do quadro de análise de variância, 2) pelo
## valor de F de um teste de Wald ou 3), pelo valor de chi-quadrado de
## um teste de razão de verossimilhanças.
## Quadro de análise de variância.
anova(m0)
## Analysis of Variance Table
##
## Response: ganhopeso
## Df Sum Sq Mean Sq F value Pr(>F)
## racoes 3 824 274.6 3.99 0.027 *
## Residuals 16 1100 68.8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##-----------------------------------------------------------------------------
## Análise de diagnóstico dos pressupostos via resíduos.
## As inferências acima só passam ser valiadas após confirmação de não
## haver afastamento dos pressupostos.
par(mfrow=c(2,2))
plot(m0); layout(1)
##-----------------------------------------------------------------------------
## Comparações múltiplas.
## Médias ajustadas (LSmeans).
p0 <- LSmeans(m0, effect="racoes", level=0.95)
## Criando a tabela com as estimativas.
p0 <- cbind(p0$grid, as.data.frame(p0$coef))
p0 <- equallevels(p0, da)
p0
## racoes estimate se df t.stat p.value lwr upr
## 1 A 26 3.708 16 7.012 2.935e-06 18.14 33.86
## 2 B 39 3.708 16 10.518 1.355e-08 31.14 46.86
## 3 C 32 3.708 16 8.630 2.047e-07 24.14 39.86
## 4 D 22 3.708 16 5.933 2.103e-05 14.14 29.86
## Gráfico de segmentos, média com IC.
## segplot(racoes~lwr+upr, centers=estimate,
## data=p0, draw=FALSE)
## Comparações múltiplas, contrastes de Tukey.
## Método de correção de p-valor: single-step.
g0 <- glht(m0, linfct=mcp(racoes="Tukey"))
summary(g0)
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lm(formula = ganhopeso ~ racoes, data = da)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## B - A == 0 13.00 5.24 2.48 0.102
## C - A == 0 6.00 5.24 1.14 0.669
## D - A == 0 -4.00 5.24 -0.76 0.870
## C - B == 0 -7.00 5.24 -1.33 0.555
## D - B == 0 -17.00 5.24 -3.24 0.024 *
## D - C == 0 -10.00 5.24 -1.91 0.264
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
## Resumo compacto com letras.
l0 <- cld(g0)
l0
## A B C D
## "ab" "b" "ab" "a"
## Formatando média seguida de letras.
p0$let <- sprintf("%.1f %s", p0$estimate, l0$mcletters$Letters)
p0$let
## [1] "26.0 ab" "39.0 b" "32.0 ab" "22.0 a"
##-----------------------------------------------------------------------------
## Representação gráfica dos resultados.
xlab <- "Ganho de peso (kg)"
ylab <- "Tipos de ração"
sub <- list(
"Médias seguidas de mesma letra indicam diferença nula à 5%.",
font=1, cex=0.8)
## Gráfico final.
segplot(racoes~lwr+upr, centers=estimate,
data=p0, draw=FALSE,
xlab=xlab, ylab=ylab, sub=sub, letras=p0$let,
panel=function(x, y, z, centers, letras, ...){
panel.segplot(x=x, y=y, z=z, centers=centers, ...)
panel.text(y=as.numeric(z), x=centers,
label=letras, pos=3)
})
##-----------------------------------------------------------------------------
## Com os níveis ordenados pela média.
p0 <- transform(p0, racoes=reorder(racoes, estimate))
levels(p0$racoes)
## [1] "D" "A" "C" "B"
p0 <- p0[order(p0$racoes),]
segplot(racoes~lwr+upr, centers=estimate,
data=p0, draw=FALSE,
xlab=xlab, ylab=ylab, sub=sub, letras=p0$let,
panel=function(x, y, z, centers, letras, ...){
panel.segplot(x=x, y=y, z=z, centers=centers, ...)
panel.text(y=as.numeric(z), x=centers,
label=letras, pos=3)
})
