Baixe R embrapa e outras Notas de estudo em PDF para Engenharia Florestal, somente na Docsity!
Introdu¸c˜ao ao sistema estat´ıstico R
Mini-curso EMBRAPA
Paulo Justiniano Ribeiro Junior
Bras´ılia, 30/05 a 03/06 de 2005 (´ultima revis˜ao: 8 de outubro de 2008)
Estas notas est˜ao dispon´ıveis em formato html em http://www.leg.ufpr.br/~paulojus/embrapa/Rembrapa e tamb´em em no arquivo arquivo em formato PDF. Este curso foi montado visando uma introdu¸c˜ao ao sistema estat´ıstico R para profissionais da EMBRAPA. O objetivo ´e ilustrar aspectos b´asicos do sistema com ˆenfase na compreens˜ao de aspectos b´asicos da linguagem, a estrutura e a forma de operar o programa. O curso n˜ao tem o objetivo de discutir em detalhe nenhum m´etodo e/ou modelo estat´ıstico em particular. M´etodos estat´ısticos b´asicos s˜ao usados ao longo do texto simplesmente para ilustrar o uso da linguagem. Ser´a assumida apenas familiaridade com conceitos e m´etodos estat´ısticos b´asicos. N˜ao ser´a as- sumido nenhum conhecimento pr´evio do R. O curso foi preparado e ser´a ministrado em ambiente LINUX por´em n˜ao far´a uso de nenhum recurso espec´ıfico deste sistema operacional e participantes poder˜ao acompanhar usando outro sistema operacional, tal como Windows®. Vamos come¸car ”experimentando o R”, para ter uma id´eia de seus recursos e a forma de trabalhar com este programa. Para isto vamos rodar e estudar os comandos mostrados no texto e seus resultados para nos familiarizar com aspectos b´asicos do programa. Ao longo deste curso iremos ver com mais detalhes o uso do programa R. Siga os seguintes passos:
- inicie o R em seu computador;
- voce ver´a uma janela de comandos com o s´ımbolo >, este ´e o prompt do R indicando que o programa est´a pronto para receber comandos;
- a seguir digite (ou ”recorte e cole”) os comandos mostrados ao longo deste material.
No restante deste texto vamos seguir as seguintes conven¸c˜oes.
comandos do R s˜ao mostrados em fontes do tipo slanted verbatim como esta, e precedidas pelo s´ımbolo >,
sa´ıdas do R s˜ao sempre exibidas em fontes do tipo verbatim como esta,
linhas iniciadas pelo s´ımbolo # s˜ao coment´arios e s˜ao ignoradas pelo R.
1 Uma primeira sess˜ao com o R
Esta ´e uma primeira sess˜ao com o R visando dar aos participantes uma id´eia geral da aparˆencia e forma de opera¸c˜ao do programa. Os comandos abaixo ser˜ao reproduzidos e comentados durante o curso. Vamos come¸car gerando dois vetores x e y de coordenadas geradas a partir de n´umeros pseudo- aleat´orios e depois inspecionar os valores gerados.
x <- rnorm(5) x
[1] 1.8614407 -1.0874200 -0.5615027 -2.3187178 0.
print(x)
[1] 1.8614407 -1.0874200 -0.5615027 -2.3187178 0.
print(x, dig = 3)
[1] 1.861 -1.087 -0.562 -2.319 0.
y <- rnorm(x) y
[1] 0.1432350 0.5101738 -0.2760532 -0.2362307 1.
args(rnorm)
function (n, mean = 0, sd = 1) NULL
No exemplo acima primeiro geramos um vetor x com 5 elementos. Note que ao fazermos y <- rnorm(x) nao especificamos o tamanho da amostra explicitamente como anteriormente mas estamos definindo um vetor y que tem o mesmo tamanho de x, por isto y foi gerado com tamb´em 5 elementos. Note que se voce tentar reproduzir este exemplo deve obter valores simulados diferentes dos mostrados aqui. Ao digitar o nome do objeto x os elementos deste objetos s˜ao exibidos. O comando print(x) tamb´em exibe os elementos do objeto por´em ´e mais flex´ıvel pois oferece op¸c˜oes extras de visualiza¸c˜ao. O comando print(x, dig=3) exibe este particular objeto x com no m´ınimo 3 d´ıgitos significativos. Para controlar o n´umero de d´ıgitos globalmente, isto ´e, para impress˜ao de qualquer objeto, por exemplo com 4 d´ıgitos, usamos options(digits=4). Neste simples exemplo introduzimos v´arias id´eias e conceitos: objeto, atribui¸c˜ao de valores, veto- res, impress˜ao de objetos, fun¸c˜ao, argumentos de fun¸c˜oes, ”defaults”, gera¸c˜ao de n´umeros aleat´orios e controle de semente. Agora vamos colocar num gr´afico os pontos gerados usando o comando
plot(x, y)
Note que a janela gr´afica se abrir´a automaticamente e exibir´a o gr´afico. H´a muitas op¸c˜oes de controle e configura¸c˜ao da janela gr´afica que s˜ao especidicadas usando-se a fun¸c˜ao par(). Algumas destas op¸c˜oes ser˜ao vistas ao longo deste material. A fun¸c˜ao plot() oferece atrav´es de seus argumentos v´arias op¸c˜oes para visualiza¸c˜ao dos gr´aficos. As argumentos e b´asicos s˜ao mostrados a seguir.
args(plot.default)
function (x, y = NULL, type = "p", xlim = NULL, ylim = NULL, log = "", main = NULL, sub = NULL, xlab = NULL, ylab = NULL, ann = par("ann"), axes = TRUE, frame.plot = axes, panel.first = NULL, panel.last = NULL, asp = NA, ...) NULL
par(mfrow = c(4, 2), mar = c(2, 2, 0.3, 0.3), mgp = c(1.5, 0.6,
plot(x, y, type = "l") plot(x, y, type = "p") plot(x, y, type = "o") plot(x, y, type = "b") plot(x, y, type = "h") plot(x, y, type = "S") plot(x, y, type = "s") plot(x, y, type = "n") par(mfrow = c(1, 1))
−2 −1 0 1 2
x
y
−2 −1 0 1 2
x
y
−2 −1 0 1 2
x
y
−2 −1 0 1 2
x
y
−2 −1 0 1 2
x
y
−2 −1 0 1 2
x
y
−2 −1 0 1 2
x
y
−2 −1 0 1 2
x
y
[1] -3 -2 -1 0 1
trunc(x)
[1] -2 -1 0 0 1
Os objetos existentes na ´area de trabalho pode ser listados usando a fun¸c˜ao ls() e objetos podem ser removidos com a fun¸c˜ao rm(). Nos comandos a seguir estamos verificando os objetos existentes na ´area de trabalho e removendo objetos que julgamos n˜ao mais necess´arios.
ls()
[1] "x" "y"
rm(x, y)
A seguir vamos criar um vetor que chamaremos de x com uma sequˆencia de n´umeros de 1 a 20. Depois criamos um vetor w de pesos com os desvios padr˜oes de cada observa¸c˜ao. Na sequˆencia montamos um data-frame de 3 colunas com vari´aveis que chamamos de x, y e w. Inspecionando o conte´udo do objeto criado digitando o seu nome. A terminamos apagando objetos que n˜ao s˜ao mais necess´arios.
x <- 1: x
[1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
w <- 1 + sqrt(x)/ w
[1] 1.500000 1.707107 1.866025 2.000000 2.118034 2.224745 2.322876 2.414214 2. [10] 2.581139 2.658312 2.732051 2.802776 2.870829 2.936492 3.000000 3.061553 3. [19] 3.179449 3.
dummy <- data.frame(x = x, y = x + rnorm(x) * w, w = w) dummy
x y w 1 1 2.148754 1. 2 2 1.659649 1. 3 3 1.711935 1. 4 4 3.111563 2. 5 5 5.342233 2. 6 6 4.383622 2. 7 7 3.954104 2. 8 8 7.896386 2. 9 9 10.505363 2. 10 10 10.535822 2. 11 11 12.522613 2. 12 12 11.747249 2. 13 13 15.556417 2. 14 14 10.148046 2. 15 15 14.245631 2. 16 16 17.722934 3. 17 17 19.053369 3. 18 18 25.597813 3. 19 19 17.851351 3. 20 20 26.432684 3.
rm(x, w)
Nos comandos a seguir estamos ajustando uma regress˜ao linear simples de y em x e examinando os resultados. Na sequˆencia, uma vez que temos valores dos pesos, podemos fazer uma regress˜ao ponderada e comparar os resultados.
fm <- lm(y ~ x, data = dummy) summary(fm)
Call: lm(formula = y ~ x, data = dummy)
Residuals: Min 1Q Median 3Q Max -5.20702 -1.20003 -0.01178 0.98924 5.
Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -1.63969 1.16188 -1.411 0. x 1.21391 0.09699 12.516 2.56e-10 ***
par(mfrow = c(2, 2)) plot(fm)
0 5 10 15 20
−
−
−
0
2
4
6
Fitted values
Residuals
Residuals vs Fitted
18
14
20
−2 −1 0 1 2
−
−
0
1
2
Theoretical Quantiles
Standardized residuals
Normal Q−Q
18
14
20
0 5 10 15 20
Fitted values
Standardized residuals
Scale−Location
14 18 20
0.00 0.05 0.10 0.
−
−
0
1
2
Leverage
Standardized residuals Cook’s distance 0.
Residuals vs Leverage
18 20
19
lrf <- lowess(x, y) plot(x, y) lines(lrf, lty = 3) abline(coef(fm)) abline(coef(fm1), lty = 2) abline(0, 1, lwd = 2) legend(1, 20, c("linear simples", "ponderada", "loess", "verdadeira"),
- lty = c(1, 2, 3, 1), lwd = c(1, 1, 1, 2))
Ao final destas an´alises removemos o objeto dummy do caminho de procura.
detach()
Agora vamos fazer um gr´afico diagn´ostico padr˜ao para checar ajuste e pressupostos: o gr´afico de res´ıduos por valores preditos e gr´afico de escores normais para checar assimetria, curtose e outliers (n˜ao muito ´util aqui).
5 10 15 20
5
10
15
20
25
x
y
linear simples ponderada loess verdadeira
par(mfrow = c(1, 2)) plot(fitted(fm), resid(fm), xlab = "Fitted values", ylab = "Residuals",
- main = "Residuals vs Fitted")
qqnorm(resid(fm), main = "Residuals Rankit Plot")
E ao final retornamos ao gr´afico padr˜ao e ”limpamos”novamente o workspace, ou seja, apagando objetos.
par(mfrow = c(1, 1)) rm(fm, fm1, lrf, dummy)
Agora vamos inspecionar dados do experimento cl´assico de Michaelson e Morley para medir a velocidade da luz. Clique para ver o arquivo morley.tab de dados no formato texto. Se quiser voce pode ainda fazer o download deste arquivo para o seu micro. Pode-se visualizar um arquivo externo dentro do pr´oprio R utilizando file.show() e note que no comando abaixo assume-se que o arquivo est´a na ´area de trabalho do R, caso contr´ario deve ser precedido do caminho para o diret´orio adequado.
file.show("morley.tab") Lendo dados como um ”data-frame”e inspecionando seu conte´udo. H´a 5 experimentos (coluna Expt) e cada um com 20 “rodadas”(coluna Run) e sl ´e o valor medido da velocidade da luz numa escala apropriada
mm <- read.table("http://www.leg.ufpr.br/~paulojus/embrapa/morley.tab") mm
Devemos definir Expt e Run como fatores tornar o data-frame vis´ıvel na posi¸c˜ao 2 do caminho de procura.
mm$Expt <- factor(mm$Expt) mm$Run <- factor(mm$Run) attach(mm)
Podemos fazer um gr´afico para comparar visualmente os 5 experimentos
plot(Expt, Speed, main = "Speed of Light Data", xlab = "Experiment No.")
Depois analisamos como um experimento em blocos ao acaso com Run e Expt como fatores e inspecionamos os resultados.
fm <- aov(Speed ~ Run + Expt, data = mm) summary(fm)
Df Sum Sq Mean Sq F value Pr(>F) Run 19 113344 5965 1.1053 0. Expt 4 94514 23629 4.3781 0.003071 ** Residuals 76 410166 5397
Signif. codes: 0 '' 0.001 '' 0.01 '' 0.05 '.' 0.1 ' ' 1
names(fm)
[1] "coefficients" "residuals" "effects" "rank" "fitted.values" [6] "assign" "qr" "df.residual" "contrasts" "xlevels" [11] "call" "terms" "model"
fm$coef
fm0$coef
(Intercept) Expt2 Expt3 Expt4 Expt 909.0 -53.0 -64.0 -88.5 -77.
mds <- tapply(Speed, Expt, mean) mds
1 2 3 4 5 909.0 856.0 845.0 820.5 831.
mds[-1] - mds[1]
2 3 4 5 -53.0 -64.0 -88.5 -77.
E este comportamento ´e controlado por options(). Por exemplo, contrastes de Helmert s˜ao definidos como se segue.
options()$contrast
unordered ordered "contr.treatment" "contr.poly"
options(contrasts = c("contr.helmert", "contr.poly")) fm0 <- update(fm,. ~. - Run) fm0$coef
(Intercept) Expt1 Expt2 Expt3 Expt 852.400 -26.500 -12.500 -12.375 -5.
mean(Speed)
[1] 852.
(mds[2] - mds[1])/
2 -26.
(2 * mds[3] - mds[1] - mds[2])/
3 -12.
(3 * mds[4] - mds[1] - mds[2] - mds[3])/
4 -12.
(4 * mds[5] - mds[1] - mds[2] - mds[3] - mds[4])/
5 -5.
Enquanto que contrastes de cada tratamento contra a m´edia geral s˜ao obtidos da forma:
options(contrasts = c("contr.sum", "contr.poly")) fm0 <- update(fm,. ~. - Run) fm0$coef
(Intercept) Expt1 Expt2 Expt3 Expt 852.4 56.6 3.6 -7.4 -31.
mds - mean(Speed)
1 2 3 4 5 56.6 3.6 -7.4 -31.9 -20.
H´a algumas op¸c˜oes de contrastes implementadas no R e al´em disto o usu´ario pode implemen- tar contrastes de sua preferˆencia. Para entender melhor os resultados acima analise as sa´ıdas dos comandos abaixo.
contr.treatment(5)
2 3 4 5 1 0 0 0 0 2 1 0 0 0 3 0 1 0 0 4 0 0 1 0 5 0 0 0 1
contr.helmert(5)
[,1] [,2] [,3] [,4] 1 -1 -1 -1 - 2 1 -1 -1 - 3 0 2 -1 - 4 0 0 3 - 5 0 0 0 4
contr.sum(5)
[,1] [,2] [,3] [,4] 1 1 0 0 0 2 0 1 0 0 3 0 0 1 0 4 0 0 0 1 5 -1 -1 -1 -
contr.poly(5)
.L .Q .C ^ [1,] -6.324555e-01 0.5345225 -3.162278e-01 0. [2,] -3.162278e-01 -0.2672612 6.324555e-01 -0. [3,] -3.287978e-17 -0.5345225 1.595204e-16 0. [4,] 3.162278e-01 -0.2672612 -6.324555e-01 -0. [5,] 6.324555e-01 0.5345225 3.162278e-01 0.
Se ainda n˜ao estivar claro experimente para cada uma destas examinar a matrix do modelo com os comandos abaixo (sa´ıdas n˜ao s˜ao mostradas aqui).
options(contrasts = c("contr.treatment", "contr.poly")) model.matrix(Speed ~ Expt) options(contrasts = c("contr.helmert", "contr.poly")) model.matrix(Speed ~ Expt) options(contrasts = c("contr.sum", "contr.poly")) model.matrix(Speed ~ Expt)
Ao final desanexamos o objeto e limpamos novamente o workspace.
detach() rm(fm, fm0)
Vamos agora ver alguns gr´aficos gerados pelas fun¸c˜oes contour() e image(). No pr´oximo exemplo x ´e um vetor de 50 valores igualmente espa¸cados no intervalo [-pi pi]. y idem. O objeto f ´e uma matrix quadrada com linhas e colunas indexadas por x e y respectivamente com os valores da fun¸c˜ao cos(y)/(1 + x^2 ).
−3 −2 −1 0 1 2 3
−
−
−
0
1
2
3
x
y
−3 −2 −1 0 1 2 3
−
−
−
0
1
2
3
x
y
th <- seq(-pi, pi, len = 100) z <- exp((0+1i) * th)
Plotando complexos significa parte imagin´aria versus real Isto deve ser um c´ırculo: Suponha que desejamos amostrar pontos dentro do c´ırculo de raio unit´ario. uma forma simples de fazer isto ´e tomar n´umeros complexos com parte real e imagin´aria padr˜ao. E depois mapeamos qualquer externo ao c´ırculo no seu rec´ıproco:
par(pty = "s") plot(z, type = "l") w <- rnorm(100) + rnorm(100) * (0+1i) w <- ifelse(Mod(w) > 1, 1/w, w)
Desta forma todos os pontos est˜ao dentro do c´ırculo unit´ario, mas a distribui¸c˜ao n˜ao ´e uniforme. Um segundo m´etodo usa a distribui¸c˜ao uniforme. os pontos devem estar melhor distribu´ıdos sobre o c´ırculo
plot(w, xlim = c(-1, 1), ylim = c(-1, 1), pch = "+", xlab = "x",
lines(z) w <- sqrt(runif(100)) * exp(2 * pi * runif(100) * (0+1i)) plot(w, xlim = c(-1, 1), ylim = c(-1, 1), pch = "+", xlab = "x",
lines(z)
Apagamos novamente os objetos...
rm(th, w, z)
... e sa´ımos do R. q()
−1.0 −0.5 0.0 0.5 1.
−1.
−0.
Re(z)
Im(z)
++ +
++
−1.0 −0.5 0.0 0.5 1.
−1.
−0.
x
y
++
−1.0 −0.5 0.0 0.5 1.
−1.
−0.
x
y
2 Estat´ıstica computacional e o sistema R
Nesta se¸c˜ao iremos seguir a apresenta¸c˜ao dispon´ıvel no arquivo estcompR.pdf
4 Introdu¸c˜ao
O programa computational R ´e gratuito, de c´odigo aberto e livremente distribu´ıdo e proporciona um ambiente para an´alises estat´ısticas. Seguem algumas informa¸c˜oes b´asicas sobre este sistema.
4.1 O projeto R
O programa R ´e gratuito e de c´odigo aberto que propicia excelente ambiente para an´alises estat´ısti- cas e com recursos gr´aficos de alta qualidade. Detalhes sobre o projeto, colaboradores, documenta¸c˜ao e diversas outras informa¸c˜oes podem ser encontradas na p´agina oficial do projeto em: http://www.r-project.org. O programa pode ser copiado livremente pela internet. H´a alguns espelhos (mirrors) brasileiros da ´area de downloads do programa chamada de CRAN (Compreensive R Arquive Network), entre eles um situado no C3SL/UFPR que pode ser acessado em http://cran.br-r-project.org Ser´a feita uma apresenta¸c˜ao r´apida da p´agina do R durante o curso onde os principais recursos ser˜ao comentados assim como as id´eias principais que governam o projeto e suas dire¸c˜oes futuras.
4.2 Um tutorial sobre o R
Al´em dos materiais dispon´ıveis na p´agina do programa h´a tamb´em um Tutorial de Introdu¸c˜ao ao R dispon´ıvel em http://www.est.ufpr.br/Rtutorial. Sugerimos aos participantes deste curso que percorram todo o conte´udo deste tutorial e retornem a ele sempre que necess´ario no decorrer do curso.
4.3 Utilizando o R
Siga os seguintes passos.
- Inicie o R em seu computador. Para iniciar o Rno LINUX basta digitar R na linha de comando.
- Vocˆe ver´a o s´ımbolo > indicando onde vocˆe ir´a digitar comandos. Este ´e o prompt do R indicando que o programa est´a pronto para receber seus comandos.
- A seguir digite (ou ”recorte e cole”) os comandos mostrados neste material. No restante deste texto vamos seguir as seguintes conven¸c˜oes:
comandos do R s˜ao sempre mostrados em fontes do tipo typewriter como esta; linhas iniciadas pelo s´ımbolo # s˜ao coment´arios e s˜ao ignoradas pelo R.
4.4 Cart˜ao de referˆencia
Para operar o R ´e necess´ario conhecer e digitar comandos. Isto pode trazer alguma dificuldade no inicio at´e que o usu´ario se familiarize com os comandos mais comuns. Uma boa forma de aprender e memorizar os comandos b´asicos ´e utilizar um Cart˜ao de Referˆencia que ´e um documento que vocˆe pode imprimir e ter sempre com vocˆe e que cont´em os comandos mais frequentemente utilizados. Aqui v˜ao trˆes op¸c˜oes:
Cart˜ao de Referˆencia em formato HTML e traduzido para portuguˆes.
Cart˜ao de Referˆencia em formato PDF preparado por Jonathan Baron.
Cart˜ao de Referˆencia em formato PDF preparado por Tom Short.
4.5 Rcmdr - “The R commander” — “menus” para o R
Para operar o R, na forma usual, ´e necess´ario conhecer e digitar comandos. Alguns usu´arios acostumados com outros programas notar˜ao de in´ıcio a falta de ”menus”. Na medida que utilizam o programa, os usu´arios (ou boa parte deles) tendem a preferir o mecanismo de comandos pois ´e mais flex´ıvel e com mais recursos. Entretanto, alguns iniciantes ou usu´arios espor´adicos poder˜ao ainda preferir algum tipo de ”menu”. O pacote Rcmdr foi desenvolvido por John Fox visando atender a esta demanda. Para utilizar este pacote basta instal´a-lo e carregar com o comando require(Rcmdr) e o menu se abrir´a automa- ticamente. Aten¸cao: Note que o Rcmdr n˜ao provˆe acesso a toda funcionalidade do R mas simplesmente a alguns procedimentos estat´ısticos mais usuais. Maiores informa¸c˜oes sobre este pacote podem ser encontradas na p´agina do Rcmdr.
