Introduç˜ao ao Ambiente Estat´ıstico R

290 Pages·2011·4.43 MB·Portuguese

Introdu¸c˜ao ao Ambiente Estat´ıstico R
Paulo Justiniano Ribeiro Junior
´
Ultima atualiza¸c˜ao: 29 de maio de 2011
Estas notas1 foram inicialmente escritas para um curso de introdu¸c˜ao ao sistema estat´ıstico
R ministrado para profissionais da EMBRAPA em Bras´ılia, 30/05 a 03/06 de 2005. Desde
sua versa˜o inicial o material tem sido constantemente modificado com a expansa˜o, corre¸ca˜o e
inclusa˜o de to´picos.
O objetivo ´e ilustrar aspectos b´asicos do sistema com ˆenfase na compreensa˜o de aspectos
ba´sicosdalinguagem,aestruturaeaformadeoperaroprograma. Nenhumm´etodoe/oumodelo
estat´ıstico em particular´e discutido em detalhes seja em seus fundamentos ou alternativas para
an´alises. Os m´etodos estat´ısticos sa˜o usados ao longo do texto simplesmente para ilustrar
aspectos do uso da linguagem.
Na maior parte do texto assume-se apenas familiaridade com conceitos e m´etodos ba´sicos
de estat´ıstica. Alguns to´picos especializados sa˜o usados em algumas Sesso˜es e, na˜o sendo de
interesse de leitor, podem ser deixados de lado sem preju´ızo ao acompanhamento das demais
partes do texto. Na˜o sera´ assumido nenhum conhecimento pr´evio do R. O curso foi preparado e
ministrado em ambiente LINUX por´em na˜o faz uso de nenhum recurso espec´ıfico deste sistema
operacional. O material pode ser acompanhado utilizando o R instalado em outros sistemas
operacionais, tal como Windowsfi ou Macintosh.
O texto come¸ca com uma Se¸ca˜o que tem como objetivo ”experimentar o R”, o que permite
ter uma id´eia de seus recursos e a forma de trabalhar com este programa. Sugere-se reproduzir
e estudar os comandos indicados bem como inspecionar e interpretar os resultados produzidos
por tais comandos, o que vai permitir uma familiaridade com aspectos b´asicos do uso do pro-
grama. Espera-se que ao final desta Se¸c˜ao o leitor se sinta a vontade para iniciar o programa
e experimentar o seu uso em outros contextos de an´alises. Ao longo do material mais detalhes
o uso do programa R sera˜o apresentados, na maior parte das vezes motivados por exemplos de
an´alise de dados.
Para utilizar o R siga os seguintes passos:
1. inicie o R em seu computador;
2. voce vera´ uma janela de comandos com o s´ımbolo >; que ´e chamado de prompt do R,
indicando que o programa esta´ pronto para receber comandos;
3. a seguir digite (ou ”recorte e cole”) os comandos mostrados ao longo deste material ou
seus pro´prios comandos.
No restante deste texto vamos seguir as seguintes conven¸co˜es.
(cid:136) comandos do R sa˜o mostrados em fontes do tipo slanted verbatim como esta, e pre-
cedidas pelo s´ımbolo >,
1Estas notas est˜ao dispon´ıveis em formato html em http://www.leg.ufpr.br/~paulojus/embrapa/
Rembrapa e tamb´em em arquivo no formato PDF.
1
2
(cid:136) sa´ıdas do R sa˜o sempre exibidas em fontes do tipo verbatim como esta,
(cid:136) linhas iniciadas pelo s´ımbolo # sa˜o coment´arios e sa˜o ignoradas pelo R.
˜
3 INTRODUC¸AO AO R
1 Uma primeira sess˜ao com o R
Esta ´e uma primeira sessa˜o com o R visando dar aos participantes uma id´eia geral da
aparˆencia e forma de opera¸ca˜o do programa. Os comandos abaixo motivam explica¸co˜es sobre
caracter´ısticas ba´sicas de linguagem e ser˜ao reproduzidos, comentados e discutidos com os
participantes durante o curso.
Vamos comec¸ar gerando dois vetores x e y de coordenadas geradas a partir de nu´meros
pseudo-aleat´orios e depois inspecionar os valores gerados.
> x <- rnorm(5)
> x
[1] -1.1584296 2.0265046 -0.9943093 1.2948114 0.9499178
> print(x)
[1] -1.1584296 2.0265046 -0.9943093 1.2948114 0.9499178
> print(x, dig = 3)
[1] -1.158 2.027 -0.994 1.295 0.950
> y <- rnorm(x)
> y
[1] 1.2299913 -0.9408620 0.9974967 -0.7250306 0.4638293
> 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.
Aodigitaronomedoobjetoxoselementosdesteobjetoss˜aoexibidos. Ocomandoprint(x)
tamb´em exibe os elementos do objeto por´em ´e mais flex´ıvel pois oferece op¸c˜oes extras de
visualizac¸˜ao. O comando print(x, dig=3) exibe este particular objeto x com no m´ınimo 3
d´ıgitos significativos. Para controlar o nu´mero de d´ıgitos globalmente, isto ´e, para impressa˜o
de qualquer objeto, por exemplo com 4 d´ıgitos, usamos options(digits=4).
Neste simples exemplo introduzimos va´rias id´eias e conceitos: objeto, atribui¸c˜ao de valores,
vetores, impress˜ao de objetos, fun¸c˜ao, argumentos de fun¸c˜oes, ”defaults”, gera¸c˜ao de nu´meros
aleat´orios e controle de semente.
Agora vamos colocar num gra´fico os pontos gerados usando o comando
> plot(x, y)
Note que a janela gra´fica se abrira´ automaticamente e exibira´ o gra´fico. H´a muitas op¸co˜es
de controle e configurac¸˜ao da janela gra´fica que sa˜o especidicadas usando-se a fun¸ca˜o par().
Algumas destas op¸co˜es sera˜o vistas ao longo deste material.
A fun¸ca˜o plot() oferece atrav´es de seus argumentos va´rias op¸c˜oes para visualizac¸˜ao dos
gra´ficos. As argumentos e ba´sicos 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
˜
1 UMA PRIMEIRA SESSAO COM O R 4
l
1.0 l
5
0. l
y
0
0.
5
0.
−
l
0 l
1.
−
−1.0 −0.5 0.0 0.5 1.0 1.5 2.0
x
Para ilustra¸ca˜o, no exemplo a seguir mostramos o uso do argumento type. Para facilitar esta
ilustrac¸˜ao vamos primeiro ordenar os valores de x e y na sequˆencia crescente dos valores de x.
> x <- sort(x)
> y <- y[order(x)]
Nos comandos abaixo iniciamos dividindo a janela gra´fica em 8 partes e reduzindo as margens
do gra´fico. A seguir produzimos diversos gr´aficos com diferentes op¸c˜oes para o argumento type.
Ao final retornamos a configurac¸˜ao original de apenas um gra´fico na janela gra´fica.
Um pouco mais sobre manipula¸c˜ao de vetores. Note que os colchetes [] sa˜o usados para
selecionar elementos e ha´ fun¸co˜es para arredondar valores.
> x
[1] -1.1584296 -0.9943093 0.9499178 1.2948114 2.0265046
> x[1]
[1] -1.15843
> x[3]
[1] 0.9499178
> x[2:4]
[1] -0.9943093 0.9499178 1.2948114
> round(x, dig = 1)
[1] -1.2 -1.0 0.9 1.3 2.0
> ceiling(x)
[1] -1 0 1 2 3
> floor(x)
[1] -2 -1 0 1 2
˜
5 INTRODUC¸AO AO R
> par(mfrow = c(4, 2), mar = c(2, 2, 0.3, 0.3), mgp = c(1.5, 0.6,
+ 0))
> 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))
1.0 1.0 l l
l
y0.0 y0.0
1.0 1.0 l l
− −
−1.0 −0.5 0.0 0.5 1.0 1.5 2.0 −1.0 −0.5 0.0 0.5 1.0 1.5 2.0
x x
1.0 l l 1.0 l l
l l
y0.0 y0.0
1.0 l l 1.0 l l
− −
−1.0 −0.5 0.0 0.5 1.0 1.5 2.0 −1.0 −0.5 0.0 0.5 1.0 1.5 2.0
x x
0 0
1. 1.
y0.0 y0.0
0 0
1. 1.
− −
−1.0 −0.5 0.0 0.5 1.0 1.5 2.0 −1.0 −0.5 0.0 0.5 1.0 1.5 2.0
x x
0 0
1. 1.
y0.0 y0.0
0 0
1. 1.
− −
−1.0 −0.5 0.0 0.5 1.0 1.5 2.0 −1.0 −0.5 0.0 0.5 1.0 1.5 2.0
x x
> trunc(x)
[1] -1 0 0 1 2
Os objetos existentes na a´rea de trabalho pode ser listados usando a fun¸ca˜o ls() e objetos
podem ser removidos com a fun¸ca˜o rm(). Nos comandos a seguir estamos verificando os objetos
existentes na a´rea de trabalho e removendo objetos que julgamos na˜o 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 nu´meros de 1
a 20. Depois criamos um vetor w de pesos com os desvios padro˜es de cada observac¸˜ao. Na
sequˆencia montamos um data-frame de 3 colunas com vari´aveis que chamamos de x, y e w.
Inspecionando o conteu´do do objeto criado digitando o seu nome. A terminamos apagando
objetos que n˜ao s˜ao mais necessa´rios.
> x <- 1:20
> x
[1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
˜
1 UMA PRIMEIRA SESSAO COM O R 6
> w <- 1 + sqrt(x)/2
> w
[1] 1.500000 1.707107 1.866025 2.000000 2.118034 2.224745 2.322876 2.414214 2.500000
[10] 2.581139 2.658312 2.732051 2.802776 2.870829 2.936492 3.000000 3.061553 3.121320
[19] 3.179449 3.236068
> dummy <- data.frame(x = x, y = x + rnorm(x) * w, w = w)
> dummy
x y w
1 1 1.5134167 1.500000
2 2 0.5108889 1.707107
3 3 3.8646610 1.866025
4 4 2.9698978 2.000000
5 5 3.6280840 2.118034
6 6 5.1372141 2.224745
7 7 4.0745573 2.322876
8 8 8.3395924 2.414214
9 9 10.0991717 2.500000
10 10 8.5394195 2.581139
11 11 14.9127871 2.658312
12 12 12.8874653 2.732051
13 13 8.5623611 2.802776
14 14 13.0209883 2.870829
15 15 13.8633465 2.936492
16 16 16.5102130 3.000000
17 17 12.2013277 3.061553
18 18 23.6649292 3.121320
19 19 17.5954538 3.179449
20 20 12.7695068 3.236068
> rm(x, w)
Nos comandos a seguir estamos ajustando uma regressa˜o linear simples de y em x e exami-
nando os resultados. Na sequˆencia, uma vez que temos valores dos pesos, podemos fazer uma
regressa˜o 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.6783 -1.1512 0.0336 1.1914 7.0518
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.1014 1.3423 0.076 0.941
x 0.9173 0.1121 8.186 1.76e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
˜
7 INTRODUC¸AO AO R
Residual standard error: 2.89 on 18 degrees of freedom
Multiple R-squared: 0.7883, Adjusted R-squared: 0.7765
F-statistic: 67.02 on 1 and 18 DF, p-value: 1.763e-07
> fm1 <- lm(y ~ x, data = dummy, weight = 1/w^2)
> summary(fm1)
Call:
lm(formula = y ~ x, data = dummy, weights = 1/w^2)
Residuals:
Min 1Q Median 3Q Max
-1.80255 -0.53011 -0.02346 0.56315 2.22021
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.07449 0.91136 -0.082 0.936
x 0.93386 0.09293 10.049 8.28e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.007 on 18 degrees of freedom
Multiple R-squared: 0.8487, Adjusted R-squared: 0.8403
F-statistic: 101 on 1 and 18 DF, p-value: 8.28e-09
Gra´ficos de res´ıduos s˜ao produzidos com plot(). Como a fun¸ca˜o produz 4 gra´ficos dividi-
remos a tela gra´fica,
Note que o comando acima par(mfrow=c(2,2)) dividiu a janela gr´afica em 4 partes para
acomodar os 4 gra´ficos. Para restaurar a configura¸ca˜o original usamos
> par(mfrow = c(1, 1))
Tornando vis´ıveis as colunas do data-frame.
> search()
[1] ".GlobalEnv" "package:tools" "package:stats" "package:graphics"
[5] "package:grDevices" "package:utils" "package:datasets" "package:methods"
[9] "Autoloads" "package:base"
> attach(dummy)
> search()
[1] ".GlobalEnv" "dummy" "package:tools" "package:stats"
[5] "package:graphics" "package:grDevices" "package:utils" "package:datasets"
[9] "package:methods" "Autoloads" "package:base"
Fazendo uma regressa˜o local na˜o-param´etrica, e visualizando o resultado. Depois adi-
cionamos a linha de regressa˜o verdadeira (intercepto 0 e inclina¸ca˜o 1), a linha da regress˜ao
sem pondera¸ca˜o e a linha de regress˜ao ponderada.
> lrf <- lowess(x, y)
> plot(x, y)
> lines(lrf, lty = 3)
> abline(coef(fm))
> abline(coef(fm1), lty = 2)
˜
1 UMA PRIMEIRA SESSAO COM O R 8
> par(mfrow = c(2, 2))
> plot(fm)
Residuals vs Fitted Normal Q−Q
3
18l 18l
s
al
5 l11 u 2
d 11l
duals l l ll l l ed resi 1 lllll
Resi 0 l lll l ll l ardiz 0 lllllllll
l l l and −1 l ll
5 St
− 20l 2
− l20
5 10 15 −2 −1 0 1 2
Fitted values Theoretical Quantiles
Scale−Location Residuals vs Leverage
duals1.5 l11 18l20l duals 23 l l18 01.5
esi0 l l esi 1
andardized r0.51. llllllllll l l andardized r −10 lllllll ll ll ll17l ll l
St l l St −2 Cook's distance 20l 0.5
0 l
0.
5 10 15 0.00 0.05 0.10 0.15
Fitted values Leverage
> 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()
Agoravamosfazerumgra´ficodiagn´osticopadra˜oparachecarajusteepressupostos: ogra´fico
de res´ıduos por valores preditos e gr´afico de escores normais para checar assimetria, curtose e
outliers (na˜o muito u´til aqui).
> 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 padra˜o e ”limpamos”novamente o workspace, ou seja, apa-
gando objetos.
> par(mfrow = c(1, 1))
> rm(fm, fm1, lrf, dummy)
˜
9 INTRODUC¸AO AO R
l
0
2
linear simples
ponderada
l
loess
l
verdadeira
5 l
1
l
l l l
y l
0 l
1
l l l
5 l
l l l
l
l
l
0
5 10 15 20
x
Residuals vs Fitted Residuals Rankit Plot
l l
6 6
l l
4 4
s
e
als 2 l l l antil 2 l l l
du l l l Qu lll
Resi 0 l l l l l l l mple 0 lllllll
l a l
2 S 2
− l − l
l l l l
4 4
− −
6 l 6 l
− −
5 10 15 −2 −1 0 1 2
Fitted values Theoretical Quantiles
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
vocepodeaindafazerodownload destearquivoparaoseumicro. Pode-sevisualizarumarquivo
externo dentro do pr´oprio R utilizando file.show() e note que no comando abaixo assume-se
que o arquivo est´a na a´rea de trabalho do R, caso contr´ario deve ser precedido do caminho para
o direto´rio adequado.
˜
1 UMA PRIMEIRA SESSAO COM O R 10
> file.show("morley.tab")
Lendo dados como um ”data-frame”e inspecionando seu conteu´do. Ha´ 5 experimentos (col-
una 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/dados/morley.tab")
> mm
Expt Run Speed
001 1 1 850
002 1 2 740
003 1 3 900
004 1 4 1070
005 1 5 930
006 1 6 850
007 1 7 950
008 1 8 980
009 1 9 980
010 1 10 880
011 1 11 1000
012 1 12 980
013 1 13 930
014 1 14 650
015 1 15 760
016 1 16 810
017 1 17 1000
018 1 18 1000
019 1 19 960
020 1 20 960
021 2 1 960
022 2 2 940
023 2 3 960
024 2 4 940
025 2 5 880
026 2 6 800
027 2 7 850
028 2 8 880
029 2 9 900
030 2 10 840
031 2 11 830
032 2 12 790
033 2 13 810
034 2 14 880
035 2 15 880
036 2 16 830
037 2 17 800
038 2 18 790
039 2 19 760
040 2 20 800
041 3 1 880
042 3 2 880
043 3 3 880