Introduç˜ao ao sistema estat´ıstico R Mini-curso EMBRAPA

287 Pages·2009·2.48 MB·Portuguese

R
Introduc¸˜ao ao sistema estat´ıstico
Mini-curso EMBRAPA
Paulo Justiniano Ribeiro Junior
Bras´ılia, 30/05 a 03/06 de 2005 (u´ltima revis˜ao: 29 de janeiro de 2009)
Estasnotasest˜aodispon´ıveisemformatohtmlemhttp://www.leg.ufpr.br/~paulojus/embrapa/Rembrapa
e tamb´em em no arquivo arquivo em formato PDF.
Este curso foi montado visando uma introduc¸˜ao ao sistema estat´ıstico R para profissionais da
EMBRAPA.O objetivo´e ilustraraspectosb´asicosdo sistemacomˆenfase na compreens˜aode 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 comec¸ar ”experimentando o R”, para ter uma id´eia de seus recursos e a forma de trabalhar
comesteprograma. Paraistovamosrodareestudaroscomandosmostradosnotextoeseusresultados
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:
1. inicie o R em seu computador;
2. 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;
3. a seguir digite (ou ”recorte e cole”) os comandos mostrados ao longo deste material.
No restante deste texto vamos seguir as seguintes convenc¸˜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
Introduc¸˜ao ao R 2
R
1 Uma primeira sess˜ao com o
Esta ´e uma primeira sess˜ao com o R visando dar aos participantes uma id´eia geral da aparˆencia
e forma de operac¸˜ao do programa. Os comandos abaixo ser˜ao reproduzidos e comentados 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.8614407 -1.0874200 -0.5615027 -2.3187178 0.3776864
> print(x)
[1] 1.8614407 -1.0874200 -0.5615027 -2.3187178 0.3776864
> print(x, dig = 3)
[1] 1.861 -1.087 -0.562 -2.319 0.378
> y <- rnorm(x)
> y
[1] 0.1432350 0.5101738 -0.2760532 -0.2362307 1.1996061
> 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.
Notequesevocetentarreproduziresteexemplodeveobtervaloressimuladosdiferentesdosmostrados
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´ıvelpois oferece opc¸˜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 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, atribuic¸˜ao de valores, veto-
res, impressa˜o de objetos, func¸˜ao, argumentos de func¸˜oes, ”defaults”, gera¸c˜ao de nu´meros 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 exibira´ o gr´afico. H´a muitas opc¸˜oes de controle
e configurac¸˜ao da janela gr´afica que s˜ao especidicadas usando-se a func¸˜ao par(). Algumas destas
opc¸˜oes ser˜ao vistas ao longo deste material.
A func¸˜ao plot() oferece atrav´es de seus argumentos v´arias opc¸˜oes para visualizac¸˜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
Introduc¸˜ao ao R 3
0
1.
5
y 0.
0
0.
−2 −1 0 1 2
x
Para ilustrac¸˜ao, no exemplo a seguir mostramos o uso do argumento type. Para facilitar esta ilus-
trac¸˜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 gr´afica em 8 partes e reduzindo as margens do
gr´afico. A seguir produzimos diversos gr´aficos com diferentes opc¸˜oes para o argumento type. Ao
final retornamos a configurac¸˜ao original de apenas um gr´afico na janela gr´afica.
Umpoucomaissobremanipulac¸˜aodevetores. Notequeoscolchetes[]s˜aousadosparaselecionar
elementos e h´a func¸˜oes para arredondar valores.
> x
[1] -2.3187178 -1.0874200 -0.5615027 0.3776864 1.8614407
> x[1]
[1] -2.318718
> x[3]
[1] -0.5615027
> x[2:4]
[1] -1.0874200 -0.5615027 0.3776864
> round(x, dig = 1)
[1] -2.3 -1.1 -0.6 0.4 1.9
> ceiling(x)
[1] -2 -1 0 1 2
> floor(x)
Introduc¸˜ao ao R 4
> 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))
0 0
1. 1.
y y
0 0
0. 0.
−2 −1 0 1 2 −2 −1 0 1 2
x x
0 0
1. 1.
y y
0 0
0. 0.
−2 −1 0 1 2 −2 −1 0 1 2
x x
0 0
1. 1.
y y
0 0
0. 0.
−2 −1 0 1 2 −2 −1 0 1 2
x x
0 0
1. 1.
y y
0 0
0. 0.
−2 −1 0 1 2 −2 −1 0 1 2
x x
[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 func¸˜ao ls() e objetos podem
ser removidos com a func¸˜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 nu´meros de 1 a 20.
Depois criamos um vetor w de pesos com os desvios padr˜oes 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
necess´arios.
> x <- 1:20
> x
Introduc¸˜ao ao R 5
[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)/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 2.148754 1.500000
2 2 1.659649 1.707107
3 3 1.711935 1.866025
4 4 3.111563 2.000000
5 5 5.342233 2.118034
6 6 4.383622 2.224745
7 7 3.954104 2.322876
8 8 7.896386 2.414214
9 9 10.505363 2.500000
10 10 10.535822 2.581139
11 11 12.522613 2.658312
12 12 11.747249 2.732051
13 13 15.556417 2.802776
14 14 10.148046 2.870829
15 15 14.245631 2.936492
16 16 17.722934 3.000000
17 17 19.053369 3.061553
18 18 25.597813 3.121320
19 19 17.851351 3.179449
20 20 26.432684 3.236068
> 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.38711
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.63969 1.16188 -1.411 0.175
x 1.21391 0.09699 12.516 2.56e-10 ***
---
Introduc¸˜ao ao R 6
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
' ' ' ' ' ' ' ' ' '
Residual standard error: 2.501 on 18 degrees of freedom
Multiple R-Squared: 0.8969, Adjusted R-squared: 0.8912
F-statistic: 156.6 on 1 and 18 DF, p-value: 2.556e-10
> 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.74545 -0.50251 0.03886 0.33719 1.87258
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.92001 0.82522 -1.115 0.280
x 1.14849 0.08414 13.649 6.18e-11 ***
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
' ' ' ' ' ' ' ' ' '
Residual standard error: 0.9119 on 18 degrees of freedom
Multiple R-Squared: 0.9119, Adjusted R-squared: 0.907
F-statistic: 186.3 on 1 and 18 DF, p-value: 6.185e-11
Gr´aficos de res´ıduos s˜ao produzidos com plot(). Como a func¸˜ao produz 4 gr´aficos dividiremos a
tela gr´afica,
Note que o comando acima par(mfrow=c(2,2)) dividiu a janela gr´afica em 4 partes para acomo-
dar os 4 gr´aficos. Para restaurar a configurac¸˜ao 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 regress˜ao local n˜ao-param´etrica, e visualizando o resultado. Depois adicionamos a
linha de regress˜ao verdadeira (intercepto 0 e inclinac¸˜ao 1), a linha da regress˜ao sem ponderac¸˜ao e a
linha de regress˜ao ponderada.
Introduc¸˜ao ao R 7
> par(mfrow = c(2, 2))
> plot(fm)
Residuals vs Fitted Normal Q−Q
6
18 18
s
4 20 ual 2 20
d
si
s 2 e 1
al r
u d
d 0 e
si z 0
e di
R 2 r
− a
d 1
n −
4 a
− t
S
14 2
6 − 14
−
0 5 10 15 20 −2 −1 0 1 2
Fitted values Theoretical Quantiles
Scale−Location Residuals vs Leverage
5 18
. 14
1 18
uals 20 uals 2 20 0.5
d d
d resi1.0 d resi 1
ze ze 0
di di
dar0.5 dar −1
n n
Sta Sta 2 19
− Cook’s distance 0.5
0
.
0
0 5 10 15 20 0.00 0.05 0.10 0.15
Fitted values Leverage
> 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 u´til aqui).
Introduc¸˜ao ao R 8
5
2
0
2
linear simples
ponderada
loess
verdadeira
5
1
y
0
1
5
5 10 15 20
x
> 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 conteu´do. 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
Introduc¸˜ao ao R 9
Residuals vs Fitted Residuals Rankit Plot
4 4
s
2 e 2
s ntil
al a
u u
d 0 Q 0
si e
Re pl
m
2 a 2
− S −
4 4
− −
0 5 10 15 20 −2 −1 0 1 2
Fitted values Theoretical Quantiles
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
Introduc¸˜ao ao R 10
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
044 3 4 860
045 3 5 720
046 3 6 720
047 3 7 620
048 3 8 860
049 3 9 970
050 3 10 950
051 3 11 880
052 3 12 910
053 3 13 850
054 3 14 870
055 3 15 840
056 3 16 840
057 3 17 850
058 3 18 840
059 3 19 840
060 3 20 840
061 4 1 890
062 4 2 810
063 4 3 810
064 4 4 820
065 4 5 800
066 4 6 770
067 4 7 760
068 4 8 740
069 4 9 750
070 4 10 760
071 4 11 910
072 4 12 920
073 4 13 890
074 4 14 860
075 4 15 880
076 4 16 720
077 4 17 840
078 4 18 850
079 4 19 850