## Objetivos:
## a) Apresentar os recursos básicos do R para análise de regressão
## b) Fundamentação em regressão linear simples
## c) Fundamentação em álgebra de matrizes
## Exemplos via algebra de matrizes
N <- c(10, 20, 30, 40, 50, 60, 70)
Pro <- c(1000, 2300, 2600, 3900, 5400, 5800, 6600)
(X <- cbind(1, N)) # Ajustando um modelo linear de 1 grau
#(X <- cbind(N)) # Ajustando um modelo linear de 1 grau forçando para a origem
#(X <- cbind(1, N, N^2)) # Ajustando um modelo linear de 2 grau
(Y <- cbind(Pro))
(XlX <- t(X) %*% X)
(XlX_inv <- solve(XlX))
(XlY <- t(X) %*% Y)
(b_est <- XlX_inv %*% XlY) # vetor das estimativas dos parâmetros do modelo ajustado: inv(X'X) * (X'Y)
(Y_est <- X %*% b_est) # vetor dos valores Y estimados pelo modelo ajustado
(res <- Y - Y_est) # vetor dos erros
plot(Pro ~ N,
xlim=c(0, 70),
ylim=c(0, 7000),
pch=19,
col='darkgreen')
lines(Y_est ~ N,
col='red')
points(Y_est ~ N,
col='red',
pch=19)
##===============================================================================
## ÊNFASE NAS FUNÇÕES DO R - USO ROTINEIRO
##===============================================================================
## Gerando dados de um experimento virtual
## x: variável independente (regressor/tratamento)
## y: variável dependente (fator resposta)
set.seed(10082006)
x <- sort(rep(1:5,
6)) # 5 tratamentos e 6 repetições
r <- rep(1:6,
5) # repetições
y <- 2 + 1 * x + rnorm(30,
mean=0,
sd=1) # definindo a variável de resposta, segundo um modelo linear
## pré-concebido (y = 2 + 1*x)
## com uma pertubação (erro experimental) intencional adicional
(dad <- data.frame(x,
r,
y))
## Visualizando os dados
with(dad,
plot(y ~ factor(x),
ylim=c(0, 10)))
## Observar que x visualmente influencia y
## ANOVA preliminar segundo delineamento inteiramente casualizado - DIC
av <- aov(y ~ factor(x),
data=dad)
summary(av)
## Pode concluir que o efeito dos tratamentos (níveis de x) é significativo.
## Como os tratamentos são quantitativos, a análise deve ser, necessariamente,
## conduzida via análise de regressão.
## Ou seja, não se pode realizar testes de comparação de médias múltiplas,
## nem análise de contrastes, para comparar os níveis de x.
## É necessário estudar se é possível estabelecer uma relação
## funcional de y em função de x. Em outras palavras, ajustar um modelo
## matemático que permita a predição de y a medida que x varia no intervalo
## considerado, além de testar a validade do modelo.
## Calculando as médias de tratamentos
dad_m <- aggregate(y ~ x,
data=dad,
mean)
dad_m
## Ajustando modelo linear e visualizando os resultados
with(dad_m,
plot(y ~ x,
ylab='m(y)',
ylim=c(0, 10),
pch=19,
col='blue'))
## Polinomio de 1 grau
reg <- lm(y ~ x,
data=dad_m)
## Polinomio de 2 grau
# reg <- lm(y ~ x + I(x^2),
# data=dad_m)
summary(reg)
ANOVAreg <- anova(reg)
ANOVAreg
## Conclusão: o r2 é elevado e 98,54% da variação na variável de resposta y,
## decorrente da variação na variável independente x, é explicada
## pelo modelo linear (y = 1.62941 + 1.09960x) ajustado.
## Em outras palavras, o modelo linear ajustado explica satisfatoriamente as
## as variações da variável dependente y em decorrência da variação na variável
## independente (ou regressora, x).
## Predizer valores y a partir de valores x
newdad <- data.frame(x=c(1.5, 2.5, 3.5, 4.5))
y_pred <- predict(reg,
newdata=newdad)
y_pred
## Visualizar médias e modelo linear ajustado
lines(fitted(reg) ~ dad_m$x,
col='red')
## Visualizar valores preditos
points(x=c(1.5, 2.5, 3.5, 4.5),
y=y_pred,
pch=19,
col='red')
Análise de Regressão
Análise didática de experimento em parcela subdividida
##=============================================================================== ## Título: Análise didática de experimento em parcela subdividida ## Versão: v4 ## Objetivos: ##=============================================================================== ## a) Mostrar diferentes alternativas de análise dos experimentos em parcelas ## subdivididas ## b) Testes de comparação de médias e análise de contrastes ##=============================================================================== # Leitura dos dados
# https://sites.google.com/site/agrobiostat/files
#=======================================
dados<-read.table('./dados/eps.txt',
header=T)
head(dados, 13)
str(dados)
is.factor(dados$bl)
is.factor(dados$adubo)
is.factor(dados$dose)
is.numeric(dados$y)
summary(dados)
#=======================================
#============== ANOVA ==================
#=======================================
av <-aov(y ~ bl + adubo * dose + Error(bl/adubo),
data=dados)
summary(av) # Notar que a interação entre adubo e dose
# é significativa
library(TukeyC)
cv(av)
#=======================================
#==== Ftab. GL adubo 3; GL resíduo 9 ===
#=======================================
(Ftab1 <- qf(5 / 100,
3,
9,
lower.tail=FALSE))
#=======================================
#Ftab. GL adubo:dose 6 ; GL resíduo 24
#=======================================
(Ftab2 <- qf(5 / 100,
6,
24,
lower.tail=FALSE))
#=======================================
#===== Ftab.GL dose 2; GL resíduo 24 ===
#=======================================
(Ftab3 <- qf(5 / 100,
2,
24,
lower.tail=FALSE))
#=======================================
#======== GRÁFICO ADUBO ================
#=======================================
curve(expr=df(x, 3, 9),
main="ADUBO - F calculado F(3, 9) gl",
xlab="F",
ylab="Densidade de probabilidades, f(F)",
xlim=c(0, 10),
n=1e3)
abline(v=Ftab1,
col="red")
abline(h=0,
lty=2)
xf <- seq(Ftab1,
10,
0.01)
ydf <- df(xf,
3,
9)
polygon(c(Ftab1, xf),
c(0, ydf),
col="red")
text(x=6,
y=.1,
paste("pf(>=",
round(Ftab1, 2),
") =",
round(pf(Ftab1,
3,
9,
lower.tail=F),
4)),
cex=0.8,
col="black")
#=======================================
#============ GRÁFICO DOSE =============
#=======================================
curve(expr=df(x, 2, 24),
main="DOSE F calculado F(2, 24) gl",
xlab="Valor F",
ylab="Densidade Probabilística F(f)",
xlim=c(0, 10),
n=1e3)
abline(v=Ftab3,
col="red")
abline(h=0,
lty=2)
xf <- seq(Ftab3,
10,
0.001)
ydf <- df(xf,
2,
24)
polygon(c(Ftab3, xf),
c(0, ydf),
col="red")
text(x=5,
y=.1,
paste("pf(>=",
round(Ftab3, 2),
") =",
round(pf(Ftab3,
2,
24,
lower.tail=F),
4)),
cex=0.8,
col="black")
#=======================================
#======== GRÁFICO ADUBO:DOSE ===========
#=======================================
curve(expr=df(x, 6, 24),
main="ADUBO:DOSE - F calculado F(6, 24) gl",
xlab="F",
ylab="Densidade Probabilística F(f)",
xlim=c(0, 10),
n=1e3)
abline(v=Ftab2,
col="red")
abline(h=0,
lty=2)
xf <- seq(Ftab2,
10,
0.001)
ydf <- df(xf,
6,
24)
polygon(c(Ftab2, xf),
c(0, ydf),
col="red")
text(x=4,
y=.1,
paste("pf(>=",
round(Ftab2, 2),
") =",
round(pf(Ftab2,
6,
24,
lower.tail=F),
4)),
cex=0.8,
col="black")
#=======================================
#============ INTERAÇÃO ==============
#=======================================
par(mfrow = c(2, 1))
with(dados,
interaction.plot(adubo,
dose,
y,
col = 1:4,
lwd = 2,
ylab = "Kg/ha"))
with(dados,
interaction.plot(dose,
adubo,
y,
col = 1:4,
lwd = 2,
ylab = "Kg/ha"))
#=======================================
#============== TUKEY ==================
#=======================================
# Testando efeito principal das parcelas (adubo)
# (finalidade didática pois a interação é significativa: necessário desdobrar!)
ta <- TukeyC(av,
which='adubo',
error='bl:adubo',
sig.level=.05)
summary(ta)
# Testando efeito principal da subparcela (dose)
# (finalidade didática pois a interação é significativa: necessário desdobrar!)
td <- TukeyC(av,
which='dose',
error='Within',
sig.level=0.05)
summary(td)
#=======================================================
# adubo/dose: estudando dose dentro dos níveis de adubo
#=======================================================
#===: nível 1
tad1 <- TukeyC.nest(av,
which='adubo:dose',
error='Within',
sig.level=0.05,
fl1=1) # 1: SALITRE C
summary(tad1)
#===: nível 2
tad2 <- TukeyC.nest(av,
which='adubo:dose',
error='Within',
sig.level=0.05,
fl1=2) # 2: SULFATO
summary(tad2)
#===: nível 3
tad3 <- TukeyC.nest(av,
which='adubo:dose',
error='Within',
sig.level=0.05,
fl1=3) # 3: URÉIA
summary(tad3)
#===: nível 4
tad4 <- TukeyC.nest(av,
which='adubo:dose',
error='Within',
sig.level=0.05,
fl1=4) # 4: CALNITRO
summary(tad4)
#=======================================
#============ CONTRASTES ===============
#=======================================
# É bem trabalhoso no caso dos experimentos em parcelas subdividadas
# devido a mais de uma estimativa do resíduo!
# Matriz de contrastes dos adubos
#A1 A2 A3 A4
mca <- rbind('A2 vs. (A1, A3, A4)' = c(-1, 3, -1, -1),
'A3 vs. (A1, A4)' = c(-1, 0, 2, -1),
'A1 vs. A4' = c( 1, 0, 0, -1))
# Matriz de contrastes das doses
#D1 D2 D3
mcd <- rbind('D3 vs. (D1, D2)' = c(-1, -1, 2),
'D1 vs. D2' = c( 1, -1, 0))
library(gmodels) # Necessário para a função: make.contrasts
## O modelo usado para a ANOVA
## não fornece todos os resultados corretos:
# av2 <-aov(y ~ bl + adubo * dose + Error(bl/adubo),
# data=dados,
# contrasts=list('adubo'=make.contrasts(mca),
# 'dose'=make.contrasts(mcd)))
##-----------------------------------------------------------------------------------
## Alternativa 1:
## Usar o modelo fatorial
## e corrigir manualmente as estimativas dos contrastes associados a subparcela
## (um pouco trabalhoso)
##-----------------------------------------------------------------------------------
## Necessário aninhar os efeitos!
av2 <-aov(y ~ bl + adubo + adubo/dose,
data=dados,
contrasts=list('adubo'=make.contrasts(mca),
'dose'=make.contrasts(mcd)))
effects(av2)[1:15]
sm <- summary(av2,
split=list('adubo'=list('A2 vs. (A1, A3, A4)' = 1,
'A3 vs. (A1, A4)' = 2,
'A1 vs. A4' = 3),
'adubo:dose'=list('A1/dose*' = c(1, 5),
'A1/d3 vs. (d1, d2)'= 1,
'A1/d1 vs. d2' = 5,
'A2/dose*' = c(2, 6),
'A2/d3 vs. (d1, d2)'= 2,
'A2/d1 vs. d2' = 6,
'A3/dose*' = c(3, 7),
'A3/d3 vs. (d1, d2)'= 3,
'A3/d1 vs. d2' = 7,
'A4/dose*' = c(4, 8),
'A4/d3 vs. (d1, d2)'= 4,
'A4/d1 vs. d2' = 8
)))
sm # Como o modelo informado não é correto,
# foi usado apenas para estimar os contrastes, observar que n erro b (subparcela):
# O qmd (141717) está errada, o correto é 74480!
# os graus de liberdade (33) estão errados, o correto é 24!
# Corrigindo os testes:
# Efeitos de dose em cada nível de adubo
round(pf(sm[[1]][c(7, 10, 13, 16), 3] / 74480, # 74480 erro associado a subparcela
2,
24, # Graus de liberdade associados a subparcela
lower.tail=FALSE),
3)
# Contrastes de dose em cada nível de adubo
round(pf(sm[[1]][c(8, 9, 11, 12, 14, 15, 17, 18), 3] / 74480,
1,
24,
lower.tail=FALSE),
3)
##-----------------------------------------------------------------------------------
## Alternativa 2:
## Com a participação do prof. Ivan Bezerra Allaman
## foi possível fazer a decomposição da SQD desse experimento corretamente!
##-----------------------------------------------------------------------------------
av3 <-aov(y ~ bl + adubo + adubo/dose + Error(bl/adubo),
data=dados,
contrasts=list('adubo'=make.contrasts(mca),
'dose'=make.contrasts(mcd)))
effects(av3) # Mensagem de erro -> método não aplicado a classe!
coef(av3) # Dá para pegar a mesma informação com a função coef!
summary(av3)
summary(av3,
split=list('adubo'=list('A2 vs. (A1, A3, A4)' = 1,
'A3 vs. (A1, A4)' = 2,
'A1 vs. A4' = 3),
'adubo:dose'=list('A1/dose*' = c(1, 5),
'A1/d3 vs. (d1, d2)'= 1,
'A1/d1 vs. d2' = 5,
'A2/dose*' = c(2, 6),
'A2/d3 vs. (d1, d2)'= 2,
'A2/d1 vs. d2' = 6,
'A3/dose*' = c(3, 7),
'A3/d3 vs. (d1, d2)'= 3,
'A3/d1 vs. d2' = 7,
'A4/dose*' = c(4, 8),
'A4/d3 vs. (d1, d2)'= 4,
'A4/d1 vs. d2' = 8
)))
Aplicação de contrastes
## Objetivos: ## Exemplificar, através do R, as análises e determinações numéricas ## referentes a análise de variância sob DIC e aplicação de ##contrastes.
===============================================================================
## Pacotes necessários:
## gplots (CRAN)
## gmodels (CRAN)
# install.packages(c('gplots',
# 'gmodels')) # Instala todos os pacotes necessários do CRAN.
##===============================
## Dados
##===============================
tra <- gl(4,
6,
labels=LETTERS[1:4])
r <- rep(1:6,
4)
pro <- c(58, 49, 51, 56, 50, 48,
60, 55, 66, 61, 54, 61,
59, 47, 44, 49, 62, 60,
45, 33, 34, 48, 42, 44)
dad <- data.frame(tra,
r,
pro)
rm(tra, r, pro) # Removendo objetos não mais necessários
##===============================
## Análise de variância
##===============================
is.data.frame(dad)
is.factor(dad$tra) # Necessário ser fator para fazer ANOVA
is.numeric(dad$pro) # Necessário ser numerico para fazer ANOVA
## Carrega pacote necessário
library(gplots)
## plotmeans (gplots)
plotmeans(pro ~ tra,
data=dad,
mean.labels=FALSE,
digits=3,
col='blue',
connect=FALSE,
ylab='Produção',
xlab='Tratamentos')
av <- aov(pro ~ tra,
data=dad)
summary(av)
##===============================
## Contrastes
##===============================
##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬
## Contrastes ortogonais
##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬
#A B C D
cmat <- rbind('(A,D) vs (B,C)' = c( 1, -1, -1, 1), # Define a matriz dos contrastes ortogonais
'A vs D' = c( 1, 0, 0, -1),
'B vs C' = c( 0, 1, -1, 0))
## Carrega pacote necessário
library(gmodels)
av1 <- aov(pro ~ tra,
data=dad,
contrasts=list(tra=make.contrasts(cmat))) # make.contrasts (gmodels): gera matriz dos contrastes
summary(av1, # ANOVA com a SQDtra e GLtra desdobrados em contrastes ortogonais
split=list(tra=list('(A,D) vs. (B,C)'=1,
'A vs. D ' =2,
'B vs. C' =3)))
##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬
## Cálculos alternativos
##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬
(tra <- gl(4,
6,
labels=LETTERS[1:4]))
#A B C D
fit.contrast(av,
tra,
c( 1, -1, -1, 1)) # fit.contrast (gmodels): testa constraste(s)
fit.contrast(av,
tra,
c( 1, 0, 0, -1))
fit.contrast(av,
tra,
c( 0, 1, -1, 0))
## Grupos de contrastes
fit.contrast(av,
tra,
rbind('(A,D) vs (B,C)' = c( 1, -1, -1, 1),
'A vs D' = c( 1, 0, 0, -1),
'B vs C' = c( 0, 1, -1, 0)))
fit.contrast(av,
tra,
cmat)
Teste de comparação de médias múltiplas
## Objetivos: ## Exemplificar, através do R, as análises e determinações numéricas ## referentes a análise de variância sob DIC e testes de comparação ## de médias usando pacotes alternativos
## Pacotes necessários:
## gplots (CRAN)
## multcomp (CRAN)
## TukeyC (CRAN)
## agricolae (CRAN)
# install.packages(c('gplots',
# 'multcomp',
# 'TukeyC',
# 'agricolae')) # Instala todos os pacotes necessários do CRAN.
##===============================
## Dados
##===============================
tra <- gl(4,
6,
labels=LETTERS[1:4])
r <- rep(1:6, 4)
pro <- c(58, 49, 51, 56, 50, 48,
60, 55, 66, 61, 54, 61,
59, 47, 44, 49, 62, 60,
45, 33, 34, 48, 42, 44)
dad <- data.frame(tra,
r,
pro)
rm(tra, r, pro) # Removendo objetos não mais necessários
##===============================
## Análise de variância
##===============================
is.data.frame(dad)
is.factor(dad$tra) # Necessário ser fator para fazer ANOVA
is.numeric(dad$pro) # Necessário ser numerico para fazer ANOVA
library(gplots)
## plotmeans (gplots)
plotmeans(pro ~ tra,
data=dad,
mean.labels=FALSE,
digits=3,
col='blue',
connect=FALSE,
ylab='Produção',
xlab='Tratamentos')
av <- aov(pro ~ tra,
data=dad)
summary(av)
#cv(av)
##===============================
## Teste de comparação de médias
##===============================
##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬
## Usando a função TukeyHSD do pacote stats
##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬
tkh <- TukeyHSD(av,
'tra')
tkh
plot(tkh)
##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬
## Usando o pacote multicomp
##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬
id.tra <- with(dad,
tapply(pro,
tra,
mean))
## Também pode ser como abaixo
# id.tra <- rep(6, 4); names(id.tra) <- LETTERS[1:4]
library(multcomp)
## contrMat (multcomp): Gera matriz de comparações
mc <- contrMat(id.tra,
type='Tukey')
mc
## glht (multcomp): Compara
tkm <- summary(glht(av,
linfct=mcp(tra=mc)))
tkm
plot(tkm)
##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬
## Usando o pacote agricolae
##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬
library(agricolae)
## HSD.test (agricolae)
with(dad,
HSD.test(pro, tra,
DFerror=df.residual(av),
MSerror=deviance(av)/df.residual(av),
alpha=0.05))
##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬
## Usando o pacote TukeyC
##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬
library(TukeyC)
## TukeyC (TukeyC)
tkc <- TukeyC(av)
summary(tkc)
summary(tkc, complete=FALSE)
plot(tkc)
## Informações adicionais
names(tkc)
tkc
tkc$Result
tkc$Diff_Prob
##====================================
## Alternativas e recursos adicionais
##====================================
## Recuperar algumas informações úteis contidas no objeto 'av'
(SQDtot <- anova(av)$Sum[1] + anova(av)$Sum[2])
(GLtot <- anova(av)$Df[1] + anova(av)$Df[2])
(SQDtra <- anova(av)$Sum[1])
(GLtra <- anova(av)$Df[1])
(SQDRes <- anova(av)$Sum[2])
(GLres <- anova(av)$Df[2])
## Cálculos alternativos de algumas medidas básicas
## Calcular médias dos tratamentos e armazenar na variável m
(m <- with(dad,
tapply(pro,
tra,
mean)))
sort(m,
d=TRUE)
## Calcular variâncias dos tratamentos e armazenar na variável vs
(vs <- with(dad,
tapply(pro,
tra,
var)))
sort(vs,
d=TRUE)
## Calcular a SQDtot em relação á média geral
(SQDtot <- with(dad,
sum((pro - mean(pro))^2)))
## Calcular os graus de liberdade total
(GLtot <- length(dad$pro) - 1)
## Coeficiente de variação
QMres <- sum(resid(av)^2) / av$df.residual
cv <- 100 * sqrt(QMres) / mean(dad$pro)
round(cv, 2)
Fundamentos da Análise de Variância
## Objetivos:
## a) Permitir a modelagem da ANOVA com número variado de grupos e ##repetições
## b) Exibir resultados numéricos e gráficos
## c) Comparar os resultados por 2 métodos diferentes de ##discriminação dos grupos
## Pacotes necessários:
## - plotrix
## - TukeyC
## - ScottKnott
##===============================================================================
## Dados para modelagem
alfa <- .08 # Probabilidade do erro tipo I na inferência da ANOVA
r <- 7 # Número de repetições de cada grupo (tratamento)
M <- c(60,
60,
# 60,
# 63,
65) # Médias dos grupos
SD <- 6 # Desvio padrão dos grupos
##===============================================================================
## Opção de cores
ocol <- rainbow(length(M))
## Funções
plotar_info_Ft <- function()
{
polygon(x=c(Ft,
seq(Ft,
20,
length=1e3)),
y=c(0,
df(seq(Ft,
20,
length=1e3),
length(M) - 1,
length(M) * (r-1))),
col='darkred')
segments(x0=Ft,
x1=Ft,
y0=0,
y1=1,
lty=3,
col='darkred')
text(x=Ft,
y=0.4,
pos=4,
offset=0,
bquote(alpha == .(alfa)),
cex=1.4,
col='darkred')
points(x=Ft,
y=-.06,
pch=17,
cex=1.5,
col='darkred')
}
plotar_info_Fc <- function()
{
polygon(x=c(Fc,
seq(Fc,
20,
length=1e3)),
y=c(0,
df(seq(Fc,
20,
length=1e3),
length(M) - 1,
length(M) * (r-1))),
col=gray(.5))
if (Fc <= 15)
{
segments(x0=Fc,
x1=Fc,
y0=0,
y1=.15,
lty=3)
}
if (Fc <= 15)
{
text(x=Fc,
y=0.12,
pos=4,
offset=0,
bquote(P == .(tab[1, 5])),
cex=1.4)
} else {
text(x=15,
y=0.12,
pos=4,
offset=0,
bquote(P == .(tab[1, 5])),
cex=1.4)
}
}
## Início da modelagem
## Gerando grupos com r repetições com distribuição normal em um data.frame
dad <- data.frame(tr=gl(length(M),
r,
labels=LETTERS[1:length(M)]),
y=rnorm(length(M) * r,
rep(M,
each=r),
sd=SD))
## Análise de variância
av <- aov(y ~ tr,
data=dad)
sav <- summary(av)
tab <- as.data.frame(round(sav[[1]][, 1:5],
2))
cv <- 100 * sqrt(tab[2, 3]) / mean(dad$y)
colnames(tab)=c('Gl',
'SQ',
'QM',
'Fcal',
'P(>Fcal)')
rownames(tab)=c('Grupos',
'Resíduo')
## Calculando médias e desvios dos grupos
medias <- with(dad,
tapply(y,
tr,
mean))
desvios <- with(dad,
tapply(y,
tr,
sd))
## O ordenamento irá facilitar a padronização de cores nos plots
ord <- order(medias)
medias <- sort(medias)
desvios <- desvios[ord]
## Valores F
Ft <- qf(alfa,
length(M) - 1,
length(M) * (r-1),
lower.tail=FALSE)
Fc <- tab[1,
4]
## Verificar dispositivo gráfico
if (is.null(dev.list()))
X11(w=8,
h=8)
## Layout do dispositivo gráfico
layout(matrix(1:6,
ncol=2,
byrow=T),
heights=c(3,
4,
4))
## Plot da tabela
plot(0.5,
type = 'n',
axes=F,
xlab='',
ylab='')
## Tabela da ANOVA no plot
suppressMessages(require(plotrix))
addtable2plot(.5,
.4,
tab,
display.rownames=TRUE,
bg=gray(.8),
cex=1.8)
mtext('ANOVA')
mtext(paste('cv = ',
round(cv,
2),
'%',
sep=''),
side=1)
## Distribuição F e decisão da inferência
curve(df(x,
length(M) - 1,
length(M) * (r-1)),
from=0,
to=20,
n=1e3,
xlab='F',
ylab='f(F)',
ylim=c(-.09,
1),
col=gray(.4))
text(x=Ft,
y=0.8,
lab=c('RAHo',
'RRHo'),
pos=c(2,
4),
offset=.5,
cex=1.4,
col=c('darkgreen',
'darkred'))
## Decisão da ANOVA
if(Fc < Ft){
points(Fc,
y=-.08,
pch=17,
cex=2,
col=gray(.5))
mtext(text = 'Não rejeita Ho',
cex=1.5,
col='darkgreen')
} else {
mtext('Rejeita Ho',
cex=1.5,
col='darkred')
if(Fc <= 15){
points(Fc,
y=-.08,
pch=17,
cex=2,
col=gray(.5))
} else {
text(15,
-.08,
pos=4,
offset=0,
expression(Delta > 15),
cex=1.3,
col=gray(.5))
}
}
## Plotando áreas sob a distribuição F
if(Fc >= Ft){
plotar_info_Ft()
plotar_info_Fc()
} else {
plotar_info_Fc()
plotar_info_Ft()
}
legend('topright',
pch=17,
col=c('darkred',
gray(.5)),
legend=c('Ftab',
'Fcal'),
cex=1.5)
means_dec <- with(dad,
reorder(tr,
-y,
mean))
boxplot(y ~ means_dec,
data=dad,
xlab='Grupos',
ylab='Y',
col=ocol)
## Plotando as densidades dos grupos
plot(1,
type='n',
ylab='Densidade',
xlab='Médias dos grupos',
xlim=c( .3 * min(medias),
1.7 * max(medias)),
ylim=c(-.005,
max(c(dnorm(medias[1:length(M)],
medias[1:length(M)],
desvios[1:length(M)])))))
for(i in 1:length(M)){
curve(dnorm(x,
medias[i],
desvios[i]),
.2 * min(medias[i]),
1.8 * max(medias[i]),
n=1e3,
add=TRUE,
col=rev(ocol)[i],
lw=1.5)
}
for(i in 1:length(M)){
segments(medias[i],
max(dnorm(medias[i],
medias[i],
desvios[i])),
medias[i],
0,
col=rev(ocol)[i],
lty=3)
}
for(i in 1:length(M)){
text(medias[i],
-.005,
substitute(mu[i],
list(i = LETTERS[ord][i])),
col=rev(ocol)[i])
}
legend('topright',
rev(LETTERS[ord]),
lty=rep(1,
length(M)),
lw=2,
col=ocol)
## Tukey
suppressMessages(require(TukeyC))
tk <- TukeyC(av,
sig.level=alfa)
plot(tk,
xlab='Grupos',
ylab='Y',
col=ocol,
title=paste('Tukey',
' (',
100*alfa,
'%',
')',
sep=''),
cex=1.2)
## ScottKnott
suppressMessages(require(ScottKnott))
sk <- SK(av,
sig.level=alfa)
plot(sk,
xlab='Grupos',
ylab='Y',
title=paste('Scott & Knott',
' (',
100*alfa,
'%',
')',
sep=''),
cex=1.2,
col=gray(seq(0,
.8,
length=length(unique(sk$groups)))))
# Coeficiente de Variação
# Arguments:
# av : an lm, aov or aovlist object
# round : rounds the values of cv to the specified number of decimal places (default 2)
cv <- function(av, round=2)
{
if(is.null(av))
stop('Please, check the parameter!')
if(inherits(av,
'lm')) {
qmee <- with(av,
sum(residuals^2) / df.residual)
m <- mean(av$fitted.values)
cv <- 100 * sqrt(qmee) / m
names(cv) <- 'cv'
return(round(cv,
round))
}
if(inherits(av,
'aovlist')) {
sm.av <- summary(av)
qmee <- sapply(sm.av,
function(i){
i[[1]]["Residuals", 3]
})
qmee <- qmee[!is.na(qmee)]
m <- av$'(Intercept)'$coefficients[[1]]
cv <- 100 * sqrt(qmee) / m
names(cv) <- paste('cv(',
letters[1:length(qmee)],
')',
sep='')
return(round(cv,
round))
}
}
Inferência Estatística - Distribuição F (Snedecor) - decisão
## Objetivos específicos # a) Gerar uma população normal; ## b) Amostrar repetidamente pares de amostras de tamanhos informados ## da população
## c) Estimar a variância a partir de cada amostra e calcular a razão ## entre as duas de cada par armazenando estes valores em um vetor.
## d) Usar este vetor para decidir, sob a distribuição F, se as ## razões anteriormente obtidas podem ser consideradas, ou não, como
## provenientes de uma mesma população, adotando um erro tipo I
## máximo para a inferência. ## Em outras palavras, o objetivo maior é ilustrar a flutuação ## normal do processo de amostragem e os possíveis erros de decisão
## na inferência influenciados por esta flutuação.
##===============================================================================
##-- Ini opções -----------------------------------------------------------------
## Tamanho das amostras
n1 <- 5 # Tamanho da amostra (n1)
n2 <- 50 # Tamanho da amostra (n2)
## Características da população
Mu <- 10
Sigma <- 2
## Número de pares de amostras
n <- 1e4 # Usar 1e2, 1e3, ..., 1e5
erro <- 5 # Define erro tipo I máximo a ser usado na inferência
# (informado em %)
##-- Fim opções -----------------------------------------------------------------
## Simulação
## Estimativas da variância a partir do tamanho n1
s2.n1 <- apply(matrix(rnorm(n * n1,
Mu,
Sigma),
ncol=n1),
1,
var)
## Estimativas da variância a partir do tamanho n2
s2.n2 <- apply(matrix(rnorm(n * n2,
Mu,
Sigma),
ncol=n2),
1,
var)
Fcal <- s2.n1 / s2.n2 # Valores da razão entre duas estimativas da variância
# de uma população normal observados
## Valor Ftab recebe valor limite da calda superior da distribuição Fi
## adotando erro tipo I máximo especificado = erro
Ftab <- qf(erro / 100,
n1-1,
n2-1,
lower.tail=FALSE)
## Inferência
dec <- Fcal < Ftab # Decisão: As amostras podem ser consideradas como
# provenientes de uma mesma população?
dec.cer <- length(dec[dec == TRUE])
dec.err <- n - dec.cer
# Plotar gráfico dos resultados observados
plot(Fcal,
ylab='Valores F calculados',
xlab='Número de pares de amostras',
font=2,
font.lab=2,
cex=0.5,
col='dark orange',
ylim=c(0, 10))
tab.inf <- parse(text=paste(paste('F[',
erro,
'~"%"]',sep=''),
'(',
n1 - 1,
', ',
n2 - 1,
')==',
round(Ftab,
2),
sep=''))
abline(h=Ftab)
text(n/2,
Ftab,
tab.inf,
font=4,
col='dark blue')
text(n/2,
Ftab/2,
paste('RAHo -> Decisões Corretas (%)',
100 * (dec.cer/n),
sep=' = '),
font=4,
col='dark green')
text(n/2,
Ftab + (max(Fcal) - Ftab)/2,
paste('RRHo -> Decisões Equivocadas (%)',
100 * (dec.err/n),
sep=' = '),
font=4,
col='red')
Inferência Estatística - Origem distribuição F (Snedecor)
## Objetivos específicos
## a) Amostrar repetidamente pares de amostras (tamanhos definidos pelo usuário)
## de uma população normal qualquer: Y ~ N(Mu, Sigma);
## b) Estimar a variância a partir de cada amostra e calcular a razão entre as
## duas estimativas;
## c) Plotar o histograma dos valores, que deve ser semelhante a função F
## (para os tamanhos de amostras defindos pelo usuário) se o número de
## repetições do processo for suficiente para se aproximar do limite;
## d) Sobrepor ao histograma a função densidade de probabilidade específica.
##===============================================================================
##-- Ini opções -----------------------------------------------------------------
## Tamanho das amostras
nN <- 4 # Tamanho da amostra (n) do Numerador (N)
nD <- 21 # Tamanho da amostra (n) do Denominador (D)
## Características da população
Mu <- 10
Sigma <- 2
## Normal padrão
# Mu <- 0
# Sigma <- 1
## Número de pares de amostras
n <- 1e4 # Usar 1e3, 1e4, 1e5
## Erro adotado na inferência
erro <- .05
#-- Fim opções ------------------------------------------------------------------
## Simulação
## Estimativas da variância do numerador
s2N <- apply(matrix(rnorm(n*nN,
Mu,
Sigma),
ncol=nN),
1,
var)
## Estimativas da variância do denominador
s2D <- apply(matrix(rnorm(n*nD,
Mu,
Sigma),
ncol=nD),
1,
var)
## Observando a distribuição da razão das estimativas da variância: observacional
require(fdth)
plot(fdt(s2N/s2D,
start=0,
end=10,
h=.05),
ty='d',
xlim=c(0, 10),
ylim=c(0, 1.2))
## Sobrepondo a curva de densidade de probabilidades: teórica
curve(df(x,
nN-1,
nD-1),
col='darkblue',
add=TRUE,
lwd=3)
# Em poucas palavras:
# A distribuição F (Snedecor) descreve como se distribui a razão
# entre duas estimativas da variâcia
# de uma distribuição normal ~N(mu, sigma) qualquer.
quantil <- qf(erro,
nN -1,
nD -1,
lower=FALSE)
# Linha decisão
segments(x0=quantil,
y0=0,
x1=quantil,
y1=1,
col='red',
lty=3)
# Texto das hipóteses
text(x=c(quantil - 1,
quantil + 1),
y=.8,
labels=c('RAH0',
'RRH0'),
col=c('darkgreen',
'red'))
# Texto do erro
text(x=quantil + 1,
y=.3,
labels=paste('Erro = ',
100*erro,
'%',
sep=''),
col='red')
# Valor do erro
text(x=quantil,
y=1.05,
label=round(quantil,
2))
Inferência Estatística - Origem distribuição F (Snedecor)
## Objetivos específicos
##===============================================================================
## a) Plotar várias funções de densidade de probabilidade da distribuição F.
## Sugere-se alterar a posição (numerador e denominador) da variável 'i'
## na função df(x,
## i,
## 20)
##===============================================================================
plot(NA,
xlab='F',
ylab='f(F)',
xlim=c(0, 5),
ylim=c(0, 1.2),
type='n')
# Número de curvas
nc <- seq(1,
5e1,
by=1)
# Define as cores
cores <- rainbow(length(nc))
for (i in nc)
{
curve(df(x,
i,
20),
col=cores[i], # Uma cor para cada curva
add=TRUE,
lwd=1)
}
Teorema central do limite
##== Ini opções =================================================================
## Pop
Mpop <- 10 # Média da pop
Vpop <- 4 # Variância da pop
## Amo
n1 <- 2 # Tamanho da amostra 1
n2 <- 10 # Tamanho da amostra 2
n3 <- 100 # Tamanho da amostra 3
nem <- 1e4 # Número de vezes que a média será estimada (usar 1e3, ..., 1e5)
##== Fim opções =================================================================
## Simulação
pop <- rnorm(nem,
mean=Mpop,
sd=sqrt(Vpop))
em1 <- apply(matrix(rnorm(nem*n1,
mean=Mpop,
sd=sqrt(Vpop)),
ncol=n1),
1,
mean)
em2 <- apply(matrix(rnorm(nem*n2,
mean=Mpop,
sd=sqrt(Vpop)),
ncol=n2),
1,
mean)
em3 <- apply(matrix(rnorm(nem*n3,
mean=Mpop,
sd=sqrt(Vpop)),
ncol=n3),
1,
mean)
## Plotar gráficos
x11(w=4, h=8)
op <- par(no.readonly=T)
par(mfcol=c(4,1),
oma=c(0,1,0,1),
mar=c(3,2,2,2))
## pop
mmY <- c(min(pop),
max(pop))
plot(pop,
main='pop',
cex=2,
pch='.',
ylab='',
xlab='',
ylim=mmY)
abline(Mpop,
0,
col='red')
## em1
plot(em1,
main=paste('n =',
n1),
cex=1,
col='red',
pch='.',
ylab='',
xlab='',
ylim=mmY)
abline(mean(em1),
0)
## em2
plot(em2,
main=paste('n =',
n2),
cex=1,
col='orange',
pch='.',
ylab='',
xlab='',
ylim=mmY)
abline(mean(em2), 0)
## em3
plot(em3,
main=paste('n =',
n3),
cex=1,
col='blue',
pch='.',
ylab='',
xlab='',
ylim=mmY)
abline(mean(em3),
0)
par(op, no.readonly=T)
## Imprimir resultados numéricos
cat('\n')
cat('pop:\n')
cat('M(pop) = ', round(mean(pop), 2)); cat('\n')
cat('V(pop) = ', round(var(pop), 2)); cat('\n\n')
cat('n1:\n')
cat('V(pop)/n1 = ', round((var(pop)/n1), 2)); cat('\n')
cat('s2(m1) = ', round(var(em1), 2)); cat('\n\n')
cat('n2:\n')
cat('V(pop)/n2 = ', round((var(pop)/n2), 2)); cat('\n')
cat('s2(m2) = ', round(var(em2), 2)); cat('\n\n')
cat('n3:\n')
cat('V(pop)/n3 = ', round((var(pop)/n3), 2)); cat('\n')
cat('s2(m2) = ', round(var(em3), 2)); cat('\n\n')
Estatísticas básicas no R
## Objetivo: Revisar tópicos de medidas estatísticas
##===================================================================
## 1- Apresentação do potencial do R
## 2- Familiarização com o R
##===================================================================
##-------------------------------------------------------------------
## Gerando amostras
##-------------------------------------------------------------------------------
A <- c(2, 1.2, 2.1, 1.6, 0.9, 2.2, 1.8) # Altura, m
B <- c(1.8, 1.7, 1.7, 1.7, 1.5, 1.7, 1.5) # Altura, m
##-------------------------------------------------------------------------------
## Visualização gráfica das amostras
##-------------------------------------------------------------------------------
par(mfrow=c(2,1),
bty='l')
## A
plot(A,
type='h',
ylim=c(0, 2.5),
xlab='A',
ylab='')
## Média de A
abline(h=mean(A),
col='red',
lw=2)
## B
plot(B,
type='h',
ylim=c(0, 2.5),
xlab='B',
ylab='Altura, m')
## Média de B
abline(h=mean(B),
col='red',
lw=2)
##-------------------------------------------------------------------
## Algumas medidas estatísticas importantes
##-------------------------------------------------------------------
## Média
mean(A)
mean(B)
## Variância
var(A)
var(B)
## Desvio padrão
sd(A)
sd(B)
## Coeficiente de variação
100 * sd(A) / mean(A)
100 * sd(B) / mean(B)
## Uma pequena e simples função
cv <- function(x) {
100 * sd(x) / mean(x)
}
cv(A)
cv(B)
## Verificando a influência da escala nas medidas
A100 <- 100 * A # (Altura, cm)
mean(A); mean(A100)
var(A); var(A100)
sd(A); sd(A100)
cv(A); cv(A100)
## Curiosidade
## Testando a igualdade
var(100 + A) == var(A)
## Testando a igualdade variando o número de casas decimais
for(ndec in 1:17)
print(round(var(100 + A), ndec) == round(var(A), ndec))
##-------------------------------------------------------------------------------
## Dividir SQD por n-1 no cálculo da estimativa da variância: por quê?
##-------------------------------------------------------------------------------
## Função para determinar o desvio quadrático médio - dqm
dqm <- function(x) {
sum((x - mean(x))^2) / length(x)
}
## Retirar 'na' amostras válidas de tamanho 'n' de uma população
## com distribuição normal Y ~ N(Mu, Sigma)
## (média 'Mu', desvio padrão 'Sigma')
## Parâmetros da população
Mu <- 10 # Média
Sigma <- 2 # Desvio Padrão
## Características da amostra
na <- 1e4 # Número de amostras (usar 1e3, 1e4, 1e5)
n <- 5 # Tamanho da amostra (usar 5, 10, 50, ...)
## Amostragem
amo <- matrix(rnorm(na * n,
Mu,
Sigma),
ncol=n)
## Parâmetros gráficos
par(mfrow=c(2,1),
pch='.',
cex=1)
## Determinando o dqm das amostras
plot(apply(amo,
1,
dqm),
ylim=c(0,30),
xlab='',
main='Desvio quadrático médio',
col='orange')
## Determinando a variância das amostras
plot(apply(amo,
1,
var),
ylim=c(0,30),
xlab='Amostra',
main='Variância',
col='darkgreen')
## Linha horizontal da média de dqm: E(dqm)
abline(h=mean(apply(amo,
1,
dqm)),
col='orange',
lw=2)
## Linha horizontal da média de var: E(var)
abline(h=mean(apply(amo,
1,
var)),
col='green',
lw=2)
##-------------------------------------------------------------------------------
## Medidas de associação
##-------------------------------------------------------------------------------
## Criando os vetores
Y1 <- 1:12
Y2 <- 10 * Y1
Y3 <- -10 * Y1
Y4 <- c(10, 24, 28, 40, 55, 62,
65, 80, 94, 95, 112, 116)
Y5 <- -1 * Y4
Y6 <- runif(n = 20, 1, 15)
Y7 <- Y6 + rnorm(20, m=2, sd=5)
Y8 <- -1 * Y6 + rnorm(20, m=2, sd=5)
Y9 <- c(0.03, 0.62, 0.07, 0.75, 0.88, 0.59,
0.93, 0.15, 0.45, 0.61, 0.33, 0.70)
Y10 <- c(0.78, 0.39, 0.40, 0.38, 0.68, 0.63,
0.66, 0.62, 0.19, 0.98, 0.75, 0.56)
## Definindo parâmetros gráficos
par(mfrow=c(2,2),
cex=1,
bty='l',
pch=19,
col='blue')
## Visualização gráfica das variáveis aleatórias em um diagrama de dispersão
## Y1 vs. Y2
plot(Y2 ~ Y1,
main='Perfeita e positiva',
sub=paste('cov(Y1, Y2)',
round(cov(Y1, Y2),
2),
sep='='))
## Y1 vs. Y3
plot(Y3 ~ Y1,
main='Perfeita e negativa',
sub=paste('cov(Y1, Y3)',
round(cov(Y1, Y3),
2),
sep='='))
## Y1 vs. Y4
plot(Y4 ~ Y1,
main='Forte e positiva',
sub=paste('cov(Y1, Y4)',
round(cov(Y1, Y4),
2),
sep='='))
## Y1 vs. Y5
plot(Y5 ~ Y1,
main='Forte e negativa',
sub=paste('cov(Y1, Y5)',
round(cov(Y1, Y5),
2),
sep='='))
## Y6 vs. Y7
plot(Y7 ~ Y6,
main='Fraca e positiva',
sub=paste('cov(Y6, Y7)',
round(cov(Y6, Y7),
2),
sep='='))
## Y6 vs. Y8
plot(Y8 ~ Y6,
main='Fraca e negativa',
sub=paste('cov(Y6, Y8)',
round(cov(Y6, Y8),
2),
sep='='))
## Y10 vs. Y9
plot(Y10 ~ Y9,
main='Nula',
sub=paste('cov(Y9, Y10)',
round(cov(Y9, Y10),
2),
sep='='))
## Covariância
cov(Y1, Y2) # Perfeita positiva
cov(Y1, Y3) # Perfeita negativa
cov(Y1, Y4) # Forte positiva
cov(Y1, Y5) # Forte negativa
cov(Y6, Y7) # Fraca e positiva
cov(Y7, Y8) # Fraca e negativa
cov(Y9, Y10) # Nula
## Correlação: observar as vantagens em relação a covariância!
cor(Y1, Y2) # Perfeita positiva
cor(Y1, Y3) # Perfeita negativa
cor(Y1, Y4) # Forte positiva
cor(Y1, Y5) # Forte negativa
cor(Y6, Y7) # Fraca e positiva
cor(Y7, Y8) # Fraca e negativa
cor(Y9, Y10) # Nula
## Verificando a influência da escala
## Covariância
cov(Y1, Y4)
cov(100 * Y1, 100 * Y4)
## Correlação
cor(Y1, Y4)
cor(100 * Y1, 100 * Y4)
Leitura e gravação de dados no R
##=============================================================================== ## Título: Leitura e gravação de dados no R ## Versão: v5 ## Objetivos: ##=============================================================================== ## a) Apresentar os recursos básicos para leitura e gravação de dados ##=============================================================================== ## Dados na WEB
## Download dos arquivos em https://sites.google.com/site/agrobiostat/files
## Local no computador
bm <- './dados/bussab_morettin.txt'
ms <- './dados/msfinal.csv'
pe <- './dados/peridon.txt'
se <- './dados/semente.csv'
tg <- './dados/tg.csv'
##===================================
## Leitura de dados remotos
##===================================
read.table(bm,
head=T,
dec=',') # observar que não foi definido o caracter usado para dado não diponível
read.table(bm,
head=T,
dec=',',
na.strings='.') # agora está OK
read.table(ms,
dec=',',
sep=';') # observar que faltou informar que a primeira linha do arquivo é o nome das variáveis
read.table(ms,
head=T,
dec=',',
sep=';') # agora está OK
read.table(se,
head=T,
dec=',',
sep=';') # um arquivo um pouco maior
read.table(tg,
head=T,
dec=',') # exemplo com um número maior de variáveis
## Observar que embora o R leia os dados remotamente, como a leitura não foi atribuída a nenhum objeto,
## Não é possível fazer nada com os dados.
## Para armazenar o objeto no espaço de trabalho (Workspace) para análises subsequente:
bm <- read.table(bm,
head=T,
dec=',',
na.strings='.') # agora está OK
ls()
str(bm)
summary(bm)
plot(bm)
##===================================
## Leitura e gravação de dados local
##===================================
## É necessário ajustar a localização dos arquivos para o computador do usuário
## Os arquivos texto estão disponíveis no LEC
read.table(bm,
head=T,
dec=',')
read.table(pe)
read.table(pe,
h=T)
dad <- read.table(pe,
h=T)
dad
## Gravação de dados
setwd('./tmp')
write.table(iris,
'iris.txt')
dad[1,1]=999
write.table(dad,
'peridon_alt.txt')
##===================================
## Salvar conteúdo de análises
##===================================
getwd() # Verificando onde está o dir. de trabalho (workdir)
dir()
sink('analise.txt') # O canal stdOUT é desviado para a conexão analise.txt
summary(iris)
library(fdth)
tb <- fdt(iris)
summary(tb)
sink() # O canal stdOUT retorna para o console do R
summary(iris)
Introdução ao R e programação
## Objetivos:
## - Apresentar os recursos gráficos básicos do R
## - Gerenciamento básico de pacotes
## - Documentação e ajuda
## - Operadores
## - Funções elementares
## - Estruturas de dados
## - Estruturas de controle de fluxo
## - Funções
## Alguns exemplos dos recursos gráficos básicos do R
demo(graphics,
echo=FALSE,
ask=TRUE) # Recursos gráficos genéricos
demo() # Lista os demos dos pacotes carregados na sessão
## Lista os demos de todos os pacotes instalados (carregados ou não)
demo(package=.packages(all.available=TRUE))
## Gerenciamento básico de pacotes
## Listar os pacotes em uso na sessão
search()
## Carregar um pacote já instalado no computador para uso
library(AER) # Applied Econometrics with R, Springer-Verlag, New York. ISBN 978-0-387-77316-2.
search()
## Instalar e carregar um pacote
install.packages('rgl')
library(rgl)
search()
## Dois demos do pacote rgl
demo(abundance,
echo=F,
ask=F) # Recursos gráficos 3D dinâmico (interagir com o mouse)
demo(bivar,
echo=F) # Recursos gráficos 3D dinâmico (interagir com o mouse)
## Atualizar todos os pacotes (deve ser feito periódicamente)
update.packages(ask=FALSE)
## Remover um pacote do caminho de busca da sessão
detach(package:rgl)
search()
## Remover (do computador) um pacote não mais necessário
remove.packages('rgl')
library(rgl) # Verificar mensagem de erro
##===================================================================
## Documentação e ajuda
##===================================================================
?round
?'for' # ou ?"for"
?'[[' # ou ?"[["
apropos('mea')
apropos('^mea') # Expressão regular
apropos('^mea',
ignore.case=F) # Expressão regular
help.search('mea')
help.start() # ou menu 'Help/Html help
vignette() # Documentos em pdf (dependente dos pacotes instalados)
vignette('grid') # Abre pdf relacionado ao pacote grid
# install.packages('sos')
library(sos)
findFn('biplot')
##===================================================================
## Operadores
##===================================================================
## Aritiméticos
##---------------------
## Operador Descrição
##---------------------
## + adição
## - subtração
## * multiplicação
## / divisião
## ^ or ** exponenciação
## x %% y verifica se a divisão é exata (módulo)
## x %/% y divisão inteira
1 + 2
4 - 2
2 * 3
5 / 2
2^3
2 ** 3
4 %% 2 # exata = 0
5 %% 2 # não exata = 1
5 %/% 2 # apenas a parte inteira da divisão
## Lógicos
##---------------------
## Operador Descrição
##---------------------
## < menor que
## <= menor ou igual que
## > maior que
## >= maior ou igual que
## == exatamente igual
## != não igual a
## !x não x
## x | y x ou y
## x & y x e y
2 < 1
2 <= 1
2 <= 2
1 > 2
1 != 1
1 != 2
!TRUE
(x <- 1:10)
x[(x < 3) | (x > 8)]
x[(x < 5) & (x > 5)]
x[(x <= 5) & (x >= 5)]
## Atribuição
##---------------------
## Operador Descrição
##---------------------
## <- recebe valor
## -> atribui valor
## = atribui valor
(x <- 1)
(2 -> x)
(x = 2)
## Extração e atribuição
##---------------------
## Operador Descrição
##---------------------
## [] usado em: vetor, matrix, array, frame e lista
## [[]] usado em: vetor, matrix, array, frame e lista
## $ usado em: frame, lista
## Extração
(x <- 1:10)
x[6]
x[c(1:2, 9:10)]
x[x >= 8]
(L <- list(x=c(4, 6, 8),
y=TRUE,
z=letters[1:6]))
L[1]
L[1] / 2 # erro proposital
L[[1]]
L[[1]] / 2
L$z[1:3]
L$z[c(1:2, 6)]
str(iris) # iris: conjunto de dados do pacote datasets
iris$Species
table(iris$Species)
## Atribuição
x[5] <- 999; x
L[1] <- NULL; L
## Definindo operadores próprios
## são funções com dois argumentos cujo nome começa e as extremidades em %
'%mop%' <- function(x, y) x * y
100 %mop% (1:10)
##===================================================================
## Algumas funções elementares
##===================================================================
## set.seed: Semente para poder reproduzir instruções aleatórias
## (runif nesse caso) sempre que necessário
set.seed(25)
(x <- round(runif(n=20,
min=0,
max=10),
digits=2))
length(x) # Número de elementos
min(x)
max(x)
sort(x)
sort(x,
dec=T)
median(x) # Mediana
mean(x) # Média
var(x) # Variância
sd(x) # Desvio padrão (standard deviation)
sqrt(var(x)) # Desvio padrão
sum(x) # Somatório
round(x)
round(x,
digits=1)
## Returns Tukey's five number summary
fivenum(x) # minimum, lower-hinge, median, upper-hinge, maximum
## Comparando resultados
fivenum(x) == quantile(x)
fivenum(x) == quantile(x,
type=5)
## Quantil arbitrário
quantile(x,
c(0, .33, .66, 1))
1:10
cumsum(1:10)
cumprod(1:10)
## Imprimir no console uma mensagem ou o valor de uma variável
print('Teste:')
x <- 10
x
print(x)
## Concatenação:
#cat('\nValor de x =',
# x)
#
#cat('\nValor de x =',
# x); cat('\n')
#
#cat('\n\tValor de x =',
# x); cat('\n')
##===================================================================
## Estruturas de dados: MUITO IMPORTANTE!!!
##===================================================================
##===============
## Vetores
##===============
## Algumas das diversas formas de criar:
c(1, 2, 3, 4, 5)
1:6
seq(from=1,
to=10,
by=1)
seq(1,
2,
length=10)
letters[1:5]
LETTERS[1:5]
c(a=1, b=5, c=10)
## Algumas formas de indexar:
(x <- seq(1, 10, by=1))
x[5:10]
x[c(5, 7:10)]
x[-(5:10)]
x > 5
x[x > 5]
x[x <= 6]
## Dar nomes aos componentes de um vetor:
names(x)
names(x) <- letters[1:length(x)]
x
x['b']
## Algumas operações básicas:
set.seed(3)
x <- round(runif(5,
0,
10),
d=1)
x
names(x) <- letters[1:length(x)]
x
x/2
x*2
x+10
sort(x)
sort(x,
dec=T)
rev(sort(x)) # idem anterior
set.seed(16)
(x <- sample(1:5,
10,
replace=T)) # Entre parênteses atribui e mostra
sort(x)
unique(x)
##===============
## Matrizes
##===============
(m <- matrix(c(1, 2, 3, 4),
nrow=2))
m[1,2]
## O produto matricial:
(x <- matrix(c(6, 7),
nrow=2))
m %*% x
## O determinante de uma matriz:
det(m)
## A transposta de uma matriz:
t(m)
## Uma matriz diagonal:
diag(c(1,2))
## A identidade da matriz:
diag(x=10,
nrow=2)
diag(nrow=2)
diag(rep(1, 3))
diag(2)
## Comandos cbind e o rbind para criar matrizes:
cbind(c(1, 2),
c(3, 4))
rbind(c(1, 3),
c(2, 4))
## O traço de uma matriz:
sum(diag(m))
## A inversa de uma matriz:
solve(m)
solve(m, x) # solução de sistemas lineares
solve(m) %*% x # ídem
##===============
## Arrays
##===============
## 2 linhas, 4 colunas, 3 dimensões
(ar <- array(letters[1:24],
c(2, 4, 3)))
class(ar)
ar[1,1,1] # ar[linha, coluna, dimensão] -> ar(x, y, z)
ar[1,1,2]
ar[1,2,3]
class(iris3) # iris3: array localizado no pacote datasets
iris3
##===============
## Fatores
##===============
set.seed(218)
(x <- factor(sample(c('a', 'b', 'c'),
5,
replace=T)))
class(x)
(l <- c('d', 'e', 'f'))
class(l)
levels(l) # observar a diferenças de fator!
set.seed(17)
(x <- factor(sample(l,
5,
rep=T)))
levels(x)
## Pode-se preferir uma tabela:
table(x)
## Se os valores estão de acordo com alguma razão, pode-se gerar níveis:
gl(n=1,
k=4)
gl(2, 4)
gl(2, 4,
labels=c(T, F))
gl(n=2,
k=1,
length=8)
gl(2, 1, 8,
labels=c(T, F))
## Pode fazer o produto cartesiano de dois fatores:
(x <- gl(2, 4))
(y <- gl(2,
1,
length=8))
interaction(x, y)
## O comando expand.grid é comparável (ele produz um frame), sendo muito útil
## para a geração de níveis de fatores para as matrizes provenientes de dados
## experimentais:
a <- c('a1', 'a2', 'a3')
b <- c('b1', 'b2')
c <- c('c1', 'c2')
dad <- expand.grid(a, b, c)
names(dad) <- c('A', 'B', 'C')
dad
##===============
## Frames
##===============
set.seed(17)
dF <- data.frame(x=rnorm(10,
m=10,
s=2),
y=sample(c(T, F),
10,
replace=T))
dF
## O comando str informa (retorna) a estrutura de um objeto e a parte dos dados
## que contém:
str(dF)
## A instrução 'summary' sumariza um objeto (aqui, um frame, mas vai bem com
## quase todos objetos):
summary(dF)
## Pode-se ter acesso aos dados das colunas de diversas maneiras:
dF
dF$x
dF[,1]
dF[['x']]
dim(dF)
names(dF)
row.names(dF)
## Pode-se mudar o nome das linhas ou das colunas:
dF
names(dF) <- LETTERS[1:ncol(dF)]
row.names(dF) <- letters[1:nrow(dF)]
dF
## Pode-se ter acesso direto as colunas de um frame usando o comando attach().
## Obs: Não deve esquecer-se de destacá-lo detach() quando terminar!
attach(dF)
table(A)
table(B)
detach(dF) # Nunca esqueça do detach após um attach!!!
## Opção melhor!
with(dF,
{
table(A)
table(B)
})
## Seleção de data.frames pelo valor das colunas
dF
subset(dF,
sub=(A > 10))
subset(dF,
sub=(B == TRUE))
subset(dF,
sub=(A > 10) & (B == TRUE))
## Ordenando um data.frame pelo valores de uma coluna
dF
(ord <- order(dF$A))
dF[ord, ]
dF2 <- dF[ord, ]
dF2
##===============
## Listas
##===============
## Gerando listas
l1 <- list()
l1[['foo']] <- 1
l1[['bar']] <- c('a', 'b', 'c')
str(l1)
l2 <- list(s1=iris,
s2=airquality,
s3=1:10,
s4=matrix(1:12,
nr=3,
nc=4))
str(l2)
## Indexando listas
l1
l1[['bar']] <- NULL
str(l1)
l1
## Exemplo, os parâmetros gráficos são armazenados em uma lista
str(par()) # a função par abre um novo dispositivo gráfico
## Indexando
l2$s3[l2$s3 >= 5 & l2$s3 < 9]
l2['s4']
l2[['s4']]
length(l2['s4']) # Observar a diferença com abaixo
length(l2[['s4']]) # Observar a diferença com acima
l2$s4
l2$s4[, 1:2]
l2$s4
l2$s4[1, 4]
l2$s4[, 4]
##===============================
## apply, tapply, lapply, sapply
##===============================
## O comando apply torna possível aplicar uma função (para o exemplo, a média)
## a cada coluna (ou linha) de um frame (ou de uma matriz):
(dF <- data.frame(x=1:10,
y=1:10,
z=1:10))
## Marginal nas linhas
apply(dF,
1,
mean)
apply(dF,
1,
sd)
## Marginal nas colunas
apply(dF,
2,
mean)
apply(dF,
2,
var)
## Em dimensões mais elevadas:
(m <- array(1:24,
dim=c(4, 3, 2)))
apply(m,
1,
mean)
mean(c( 1, 5, 9,
13, 17, 21)) # Observe!
apply(m,
c(1, 2),
mean)
mean(c(1, 13)) # Observe!
## A função tapply permite reagrupar as observações de acordo com o valor dos
## fatores e uma função (média, soma, etc..) para cada grupo obtido:
tapply(1:10,
gl(2, 5),
mean)
tapply(1:10,
gl(2, 5),
sum)
by(1:10,
gl(2, 5),
mean)
by(1:10,
gl(2, 5),
sum)
## A função sapply aplica a cada elemento de uma lista (ou de um vetor, etc..) e
## se possível retorna um vetor:
set.seed(2012)
x <- list(a=rnorm(10),
b=runif(100),
c=rgamma(50, 1))
sapply(x, sd)
## A função lapply faz a mesma coisa, mas retorna uma lista:
lapply(x, sd)
##===============================================================================
## Estruturas de controle
##===============================================================================
set.seed(15)
(x <- round(rnorm(10)))
(y <- ifelse(x > 0,
1,
-1))
(z <- ifelse(x > 0,
1,
ifelse(x < 0,
'< zero',
0)))
##===============
## Conexão:
##===============
x <- 'a'
y <- switch(x,
a='Bonjour',
b='Gutten Tag',
c='Hello',
d='Konnichi wa')
y
##===============
## Loop for:
##===============
for (i in 1:10)
print(i)
for (i in 2:9)
ifelse (i <= 5,
print(i),
print(paste(i,
'5',
sep=' > ')))
##===============
## Loop while:
##===============
a <- 0
while (a < 11) {
if (a >= 3)
print(a)
else
cat('não\n')
a <- a + 1 # Expressão avaliada..
}
##===============
## Loop repeat:
##===============
a <- 0
repeat {
a <- a + 1
if (a >= 3)
print(a)
else
cat('não\n')
if (a == 10) break
}
##===============================================================================
## Funções
##===============================================================================
f <- function(x) x/10 + 1
f(x=10)
f(10) # Chamada alternativa
f <- function(x) {
x/10 + 1
}
f(x=10)
f(10) # Chamada alternativa
## Pode atribuir valores aos argumentos:
f <- function(x, y=3) {
x/10 + 1 - y
}
f(10)
## Na chamada da função, pode-se usar o nome dos argumentos, passar novos valores
## para as variáveis, não sendo necessário que os mesmos sigam a ordem declarada
## na função (desde que os valores sejam acompanhados dos respectivos nomes):
f(y=1, x=10)
f <- function(x, y) {
x/10 + 1 - y
}
f(1, 10)
f(10, 1)
## No fim dos argumentos, pode haver três pontos, representando todos os
## argumentos não especificados:
f <- function(x,
y,
cex1=1,
cex2=3, ...) {
par(mfrow=c(2,1))
plot(x, cex=cex1, ...)
plot(y, cex=cex2, ...)
}
Introdução ao R
##Ajuda
help(mean)
help.start()
## Instalando pacotes ##
library()
help.start()
install.packages("vegan")
## Carregando um pacote já instalado: Demonstração com pacote MASS
x1 <- rnorm(n=15, mean=1, sd=3)
hist(x1)
truehist(x1)# erro: Pacote não carregado
search()
library(MASS)
truehist(x1)
## Listando e removendo objetos ##
## O que tem na area de trabalho?
ls()
## Apagar todos os objetos com
rm( list = ls() )
## Criando um vetor de 3 algarismos concatenados
A1 <- c(1,2,3)
A1 # digite o nome do objeto para exibir seu conteúdo
## Existe uma maneira + simples:
## coloque o comando de atribuição entre parenteses para exibir seu resultado
(A2 <- c(10,20,30)) #atribuição entre parenteses
(b <- c(A1,A2))
ls()
help(ls)
ls(pattern="A")
## Para apagar objetos com um certo padrão
rm(list=ls(pattern="A"))
## ou
rm(A1,A2)
## Classes de objetos ##
copa.70 <- "21/06/70"
copa.94 <- "17/07/94"
## diferença
copa.94 - copa.70 # não funciona
## Classes não adequadas:
class(copa.70)
class(copa.94)
## Mudando para a classe de data
copa.70 <- as.Date(copa.70,format="%d/%m/%y")
copa.94 <- as.Date(copa.94,format="%d/%m/%y")
class(copa.70)
class(copa.94)
copa.94 - copa.70
##Gerando dados
area <- c(303, 379, 961, 295, 332, 47, 122, 11, 53, 2749)
riqueza <- c(3, 10, 20, 7, 8, 4, 8, 3, 5, 23)
area
riqueza
summary(area)
summary(riqueza)
mean(x=area)
var(area)
sd(x=area)
mean(riqueza)
var(riqueza)
sd(riqueza)
plot(x=area, y=riqueza, xlab="Area (ha)", ylab="Número de Espécies")
modelo1 <- lm(riqueza~area)
plot(x=area, y=riqueza, xlab="Area (ha)", ylab="Número de Espécies")
abline(modelo1)
plot(x=area, y=riqueza, xlab="Log Area (ha)", ylab="Log Número de Espécies", log="xy")
modelo2 <- lm(log(riqueza,base=10)~log(area,base=10))
abline(modelo2)
#######################
## Funções Matemáticas
log( 10 )
# Logaritmo natural
log( 10, base = 10) # Log base 10
log10(10)
# Também log de base 10
log( 10, base = 3.4076) # base 3.4076
############################
#Constante de funções trigonométricas
sin(0.5*pi) # Seno
cos(2*pi) # Coseno
asin(1) # Arco seno (radianos)
asin(1) / pi * 180
#############################
1 - (1 + 10^(-15))
factorial(100) # Fatorial de 100
#############################
## arredondamentos
round( 4.3478 )
round( 4.3478 , digits=3)
round( 4.3478 , digits=2)
#######################################
#Valores Infinitos, Indefinidos e Inexistentes
-5/0
500000000000000000/Inf
2 * NA
#########################
#Vetores
a = c(3.4, pi, exp(-1))
a1= c(3.4,pi, "a")
a2= c(3.4,pi, a)
a1
a2
#Sequências
b = 1:8
b
seq(from=1, to=4)
seq(from=1, to=4, by=0.5)
seq(from=1, to=4, length=6)
## Sequencias com padrão
rep(5, times=3)
rep(1:5, 3)
rep(1:5,each=3)
####################
#Operações com vetores
a = seq(0,8,2)
a
b = c(1,15,18,3,6)
a+b
a^(1/b)
length(b)/length(a)
b
sort(b)
sort(b, decreasing=T)
#######################################
#P. gonoacantha é uma espécie arbórea com média de 40 cm e 15 cm de desvio-padrão do diâmetro à altura do peito (DAP), podendo atingir até 90 cm de DAP quando adulta
norm.ale=function(x){dnorm(x,mean=44, sd=15)}
curve(dnorm(x, mean=45, sd=15),from=5, to=90, xlab="DAP(cm)", ylab= "densidade probabilistica", col="blue")
## quantos tem DAP maior que 50
pnorm(q=50,mean=45,sd=15,lower.tail=FALSE)
## visualizar no gráfico
abline(v=50, col="red", lty=2)
# Qual o DAP que acomoda 95% da população?
qnorm(p=0.95,mean=45,sd=15)
# inserir no gráfico
abline(v=69.6728, col="green", lty=2)
help(mean)
help.start()
## Instalando pacotes ##
library()
help.start()
install.packages("vegan")
## Carregando um pacote já instalado: Demonstração com pacote MASS
x1 <- rnorm(n=15, mean=1, sd=3)
hist(x1)
truehist(x1)# erro: Pacote não carregado
search()
library(MASS)
truehist(x1)
## Listando e removendo objetos ##
## O que tem na area de trabalho?
ls()
## Apagar todos os objetos com
rm( list = ls() )
## Criando um vetor de 3 algarismos concatenados
A1 <- c(1,2,3)
A1 # digite o nome do objeto para exibir seu conteúdo
## Existe uma maneira + simples:
## coloque o comando de atribuição entre parenteses para exibir seu resultado
(A2 <- c(10,20,30)) #atribuição entre parenteses
(b <- c(A1,A2))
ls()
help(ls)
ls(pattern="A")
## Para apagar objetos com um certo padrão
rm(list=ls(pattern="A"))
## ou
rm(A1,A2)
## Classes de objetos ##
copa.70 <- "21/06/70"
copa.94 <- "17/07/94"
## diferença
copa.94 - copa.70 # não funciona
## Classes não adequadas:
class(copa.70)
class(copa.94)
## Mudando para a classe de data
copa.70 <- as.Date(copa.70,format="%d/%m/%y")
copa.94 <- as.Date(copa.94,format="%d/%m/%y")
class(copa.70)
class(copa.94)
copa.94 - copa.70
##Gerando dados
area <- c(303, 379, 961, 295, 332, 47, 122, 11, 53, 2749)
riqueza <- c(3, 10, 20, 7, 8, 4, 8, 3, 5, 23)
area
riqueza
summary(area)
summary(riqueza)
mean(x=area)
var(area)
sd(x=area)
mean(riqueza)
var(riqueza)
sd(riqueza)
plot(x=area, y=riqueza, xlab="Area (ha)", ylab="Número de Espécies")
modelo1 <- lm(riqueza~area)
plot(x=area, y=riqueza, xlab="Area (ha)", ylab="Número de Espécies")
abline(modelo1)
plot(x=area, y=riqueza, xlab="Log Area (ha)", ylab="Log Número de Espécies", log="xy")
modelo2 <- lm(log(riqueza,base=10)~log(area,base=10))
abline(modelo2)
#######################
## Funções Matemáticas
log( 10 )
# Logaritmo natural
log( 10, base = 10) # Log base 10
log10(10)
# Também log de base 10
log( 10, base = 3.4076) # base 3.4076
############################
#Constante de funções trigonométricas
sin(0.5*pi) # Seno
cos(2*pi) # Coseno
asin(1) # Arco seno (radianos)
asin(1) / pi * 180
#############################
1 - (1 + 10^(-15))
factorial(100) # Fatorial de 100
#############################
## arredondamentos
round( 4.3478 )
round( 4.3478 , digits=3)
round( 4.3478 , digits=2)
#######################################
#Valores Infinitos, Indefinidos e Inexistentes
-5/0
500000000000000000/Inf
2 * NA
#########################
#Vetores
a = c(3.4, pi, exp(-1))
a1= c(3.4,pi, "a")
a2= c(3.4,pi, a)
a1
a2
#Sequências
b = 1:8
b
seq(from=1, to=4)
seq(from=1, to=4, by=0.5)
seq(from=1, to=4, length=6)
## Sequencias com padrão
rep(5, times=3)
rep(1:5, 3)
rep(1:5,each=3)
####################
#Operações com vetores
a = seq(0,8,2)
a
b = c(1,15,18,3,6)
a+b
a^(1/b)
length(b)/length(a)
b
sort(b)
sort(b, decreasing=T)
#######################################
#P. gonoacantha é uma espécie arbórea com média de 40 cm e 15 cm de desvio-padrão do diâmetro à altura do peito (DAP), podendo atingir até 90 cm de DAP quando adulta
norm.ale=function(x){dnorm(x,mean=44, sd=15)}
curve(dnorm(x, mean=45, sd=15),from=5, to=90, xlab="DAP(cm)", ylab= "densidade probabilistica", col="blue")
## quantos tem DAP maior que 50
pnorm(q=50,mean=45,sd=15,lower.tail=FALSE)
## visualizar no gráfico
abline(v=50, col="red", lty=2)
# Qual o DAP que acomoda 95% da população?
qnorm(p=0.95,mean=45,sd=15)
# inserir no gráfico
abline(v=69.6728, col="green", lty=2)
Subscribe to:
Posts (Atom)