Esta é uma primeira sessão com o R visando dar aos participantes uma ideia geral da aparência e forma de operação do programa. Os comandos abaixo motivam explicações sobre características básicas de linguagem e serão reproduzidos, comentados e discutidos com os participantes durante o curso.
Vamos começar gerando dois vetores x e y de coordenadas geradas a partir de números pseudo-aleatórios e depois inspecionar os valores gerados.
> x <- rnorm(5) > x
[1] -0.2068510 -1.2693623 0.5574482 2.0031405 0.3570333
> print(x)
[1] -0.2068510 -1.2693623 0.5574482 2.0031405 0.3570333
> print(x, dig=3)
[1] -0.207 -1.269 0.557 2.003 0.357
> y <- rnorm(x) > y
[1] 0.3027666 0.1811900 0.9989794 -0.9324467 -0.2098432
> 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) não 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ém 5 elementos. Note que se você tentar reproduzir este exemplo deve obter valores simulados diferentes dos mostrados aqui.
Ao digitar o nome do objeto x os elementos deste objetos são exibidos. O comando print(x) também exibe os elementos do objeto porém é mais flexível pois oferece opções extras de visualização. O comando print(x, dig=3) exibe este particular objeto x com no mínimo 3 dígitos significativos. Para controlar o número de dígitos globalmente, isto é, para impressão de qualquer objeto, por exemplo com 4 dígitos, usamos options(digits=4).
Neste simples exemplo introduzimos várias ideias e conceitos: objeto, atribuição de valores, vetores, impressão de objetos, função, argumentos de funções, defaults, geração de números aleatórios e controle de semente.
Agora vamos colocar num gráfico os pontos gerados usando o comando
> plot(x,y)
Note que a janela gráfica se abrirá automaticamente e exibirá o gráfico. Há muitas opções de controle e configuração da janela gráfica que são especificadas usando-se a função par(). Algumas destas opções serão vistas ao longo deste material.
A função plot() oferece através de seus argumentos várias opções para visualização dos gráficos. Os argumentos básicos são mostrados a seguir.
> args(plot.default)
function (x, y = NULL, type = "p", xlim = NULL, ylim = NULL, log = "", lim2 = FALSE, main = NULL, sub = NULL, xlab = NULL, ylab = NULL, ann = par("ann"), axes = TRUE, frame.plot = axes, panel.first = NULL, panel.last = NULL, asp = NA, xgap.axis = NA, ygap.axis = NA, ...) NULL
Para ilustração, no exemplo a seguir mostramos o uso do argumento type. Para facilitar esta ilustração vamos primeiro ordenar os valores de x e y na sequência crescente dos valores de x.
> x <- sort(x) > y <- y[order(x)]
Nos comandos abaixo iniciamos dividindo a janela gráfica em 8 partes e reduzindo as margens do gráfico. A seguir produzimos diversos gráficos com diferentes opções para o argumento type. Ao final retornamos a configuração original de apenas um gráfico na janela gráfica.
Um pouco mais sobre manipulação de vetores. Note que os colchetes [] são usados para selecionar elementos e há funções para arredondar valores.
> x
[1] -1.2693623 -0.2068510 0.3570333 0.5574482 2.0031405
> x[1]
[1] -1.269362
> x[3]
[1] 0.3570333
> x[2:4]
[1] -0.2068510 0.3570333 0.5574482
> round(x, dig=1)
[1] -1.3 -0.2 0.4 0.6 2.0
> ceiling(x)
[1] -1 0 1 1 3
> floor(x)
[1] -2 -1 0 0 2
> trunc(x)
[1] -1 0 0 0 2
Os objetos existentes na área de trabalho pode ser listados usando a função ls() e objetos podem ser removidos com a função rm(). Nos comandos a seguir estamos verificando os objetos existentes na área de trabalho e removendo objetos que julgamos não mais necessários.
> ls()
[1] "x" "y"
> rm(x, y)
A seguir vamos criar um vetor que chamaremos de x com uma sequência de números de 1 a 20. Depois criamos um vetor w de pesos com os desvios padrões de cada observação. Na sequência montamos um data.frame de 3 colunas com variáveis que chamamos de x, y e w. Inspecionando o conteúdo do objeto criado digitando o seu nome. A terminamos apagando objetos que não são mais necessá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
> 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 0.6493015 1.500000 2 2 4.0734352 1.707107 3 3 -0.3563404 1.866025 4 4 3.7202244 2.000000 5 5 6.0038068 2.118034 6 6 6.9167381 2.224745 7 7 5.1757240 2.322876 8 8 5.9294215 2.414214 9 9 5.1512831 2.500000 10 10 13.5005700 2.581139 11 11 12.9993392 2.658312 12 12 12.7587511 2.732051 13 13 13.9086405 2.802776 14 14 11.0851592 2.870829 15 15 14.4183039 2.936492 16 16 14.4569368 3.000000 17 17 14.3904028 3.061553 18 18 17.3722979 3.121320 19 19 20.9927450 3.179449 20 20 13.3542118 3.236068
> rm(x,w)
Nos comandos a seguir estamos ajustando uma regressão linear simples de y em x e examinando os resultados. Na sequência, uma vez que temos valores dos pesos, podemos fazer uma regressã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.0262 -1.5341 0.1448 1.6452 4.1258 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.36916 1.14083 0.324 0.75 x 0.90056 0.09524 9.456 2.1e-08 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 2.456 on 18 degrees of freedom Multiple R-squared: 0.8324, Adjusted R-squared: 0.8231 F-statistic: 89.42 on 1 and 18 DF, p-value: 2.097e-08
> fm1 <- lm(y ~ x, data=dummy, weight=1/w^2) > summary(fm1)
Call: lm(formula = y ~ x, data = dummy, weights = 1/w^2) Weighted Residuals: Min 1Q Median 3Q Max -1.74670 -0.62110 0.05046 0.60417 1.60510 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.13682 0.87043 0.157 0.877 x 0.92207 0.08875 10.389 4.95e-09 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.9619 on 18 degrees of freedom Multiple R-squared: 0.8571, Adjusted R-squared: 0.8491 F-statistic: 107.9 on 1 and 18 DF, p-value: 4.948e-09
Gráficos de resíduos são produzidos com plot(). Como a função produz 4 gráficos dividiremos a tela gráfica,
> par(mfrow=c(2,2)) > plot(fm)
Note que o comando acima par(mfrow=c(2,2)) dividiu a janela gráfica em 4 partes para acomodar os 4 gráficos. Para restaurar a configuração original usamos
> par(mfrow=c(1,1))
Tornando visíveis as colunas do data.frame.
> search()
[1] ".GlobalEnv" "package:knitr" "package:stats" "package:graphics" [5] "package:grDevices" "package:utils" "package:datasets" "package:methods" [9] "Autoloads" "package:base"
> attach(dummy) > search()
[1] ".GlobalEnv" "dummy" "package:knitr" "package:stats" [5] "package:graphics" "package:grDevices" "package:utils" "package:datasets" [9] "package:methods" "Autoloads" "package:base"
Fazendo uma regressão local não-paramétrica, e visualizando o resultado. Depois adicionamos a linha de regressão verdadeira (intercepto 0 e inclinação 1), a linha da regressão sem ponderação e a linha de regressão ponderada.
> 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("topleft", c("linear simples","ponderada","loess","verdadeira"), lty=c(1,2,3,1), lwd=c(1,1,1,2))
Ao final das análises removemos o objeto dummy do caminho de procura.
> detach()
Agora vamos fazer um gráfico diagnóstico padrão para checar ajuste e pressupostos: o gráfico de resíduos por valores preditos e gráfico de escores normais para checar assimetria, curtose e outliers (não muito ú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áfico padrão 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ássico de Michaelson e Morley para medir a velocidade da luz. Clique para ver o arquivo morley.tab de dados no formato texto. Se quiser você pode ainda fazer o download deste arquivo para o seu micro. Pode-se visualizar um arquivo externo dentro do próprio R utilizando file.show() e note que no comando abaixo assume-se que o arquivo está na área de trabalho do R, caso contrário deve ser precedido do caminho para o diretório adequado.
> file.show("morley.tab")
Vamos importar os dados como um data.frame e inspecionar seu conteúdo listando as seis
primeiras linhas do objeto criado com o comando head(). Há 5 experimentos (coluna Expt) e cada
um com 20 “rodadas”(coluna Run) e sl é o valor medido da velocidade da luz numa escala apropriada
> mm <- read.table("http://www.leg.ufpr.br/~paulojus/dados/morley.tab") > head(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
Devemos definir Expt e Run como fatores e tornar o data.frame visível na posição 2 do caminho de procura.
> mm$Expt <- factor(mm$Expt) > mm$Run <- factor(mm$Run) > attach(mm)
Podemos fazer um gráfico 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.105 0.36321 Expt 4 94514 23629 4.378 0.00307 ** 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
(Intercept) Run2 Run3 Run4 Run5 Run6 9.50600e+02 -5.20000e+01 -2.80000e+01 6.00000e+00 -7.60000e+01 -1.04000e+02 Run7 Run8 Run9 Run10 Run11 Run12 -1.00000e+02 -4.00000e+01 -1.00000e+01 -3.80000e+01 4.00000e+00 -8.70878e-14 Run13 Run14 Run15 Run16 Run17 Run18 -3.60000e+01 -9.40000e+01 -6.00000e+01 -6.60000e+01 -6.00000e+00 -3.80000e+01 Run19 Run20 Expt2 Expt3 Expt4 Expt5 -5.00000e+01 -4.40000e+01 -5.30000e+01 -6.40000e+01 -8.85000e+01 -7.75000e+01
Podemos redefinir o modelo, por exemplo ajustando um sub-modelo sem o fator “runs” e comparar os dois modelos lineares via uma análise de variância.
> fm0 <- update(fm, . ~ . - Run) > anova(fm0, fm)
Analysis of Variance Table Model 1: Speed ~ Expt Model 2: Speed ~ Run + Expt Res.Df RSS Df Sum of Sq F Pr(>F) 1 95 523510 2 76 410166 19 113344 1.1053 0.3632
É importante saber interpretar os coeficientes segundo a parametrização utilizada. Por default a parametrização é feita tomando o primeiro grupo como referência.
> fm0$coef
(Intercept) Expt2 Expt3 Expt4 Expt5 909.0 -53.0 -64.0 -88.5 -77.5
> (mds <- tapply(Speed, Expt, mean))
1 2 3 4 5 909.0 856.0 845.0 820.5 831.5
> mds[-1] - mds[1]
2 3 4 5 -53.0 -64.0 -88.5 -77.5
E este comportamento é controlado por options(). Por exemplo, contrastes de Helmert são 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 Expt4 852.400 -26.500 -12.500 -12.375 -5.225
> mean(Speed)
[1] 852.4
> (mds[2] - mds[1])/2
2 -26.5
> (2*mds[3] - mds[1] - mds[2])/6
3 -12.5
> (3*mds[4] - mds[1] - mds[2] - mds[3])/12
4 -12.375
> (4*mds[5] - mds[1] - mds[2] - mds[3] - mds[4])/20
5 -5.225
Enquanto que contrastes de cada tratamento contra a média geral são obtidos da forma:
> options(contrasts=c("contr.sum", "contr.poly")) > fm0 <- update(fm, . ~ . - Run) > fm0$coef
(Intercept) Expt1 Expt2 Expt3 Expt4 852.4 56.6 3.6 -7.4 -31.9
> mds - mean(Speed)
1 2 3 4 5 56.6 3.6 -7.4 -31.9 -20.9
Há algumas opções de contrastes implementadas no R e além disto o usuário pode implementar contrastes de sua preferência. 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 -1 2 1 -1 -1 -1 3 0 2 -1 -1 4 0 0 3 -1 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 -1
> contr.poly(5)
.L .Q .C ^4 [1,] -6.324555e-01 0.5345225 -3.162278e-01 0.1195229 [2,] -3.162278e-01 -0.2672612 6.324555e-01 -0.4780914 [3,] -3.510833e-17 -0.5345225 1.755417e-16 0.7171372 [4,] 3.162278e-01 -0.2672612 -6.324555e-01 -0.4780914 [5,] 6.324555e-01 0.5345225 3.162278e-01 0.1195229
Se ainda não estiver claro experimente para cada uma destas examinar a matriz do modelo com os comandos abaixo (saídas não são 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áficos gerados pelas funções contour() e image().
No próximo exemplo x é um vetor de 50 valores igualmente espaçados no intervalo \([-\pi , \pi ]\). O objeto y é
uma cópia de x. O objeto f é uma matriz quadrada com linhas e colunas indexadas por x e y
respectivamente com os valores da função \(cos(y)/(1 + x^2)\).
> x <- seq(-pi, pi, len=50) > y <- x > f <- outer(x, y, function(x, y) cos(y)/(1 + x^2))
Agora gravamos parâmetros gráficos definindo a região gráfica como quadrada e fazemos um mapa de contorno de f. Depois adicionamos mais linhas para melhor visualização. fa é a “parte assimétrica” e t() é transposição. Ao final restauramos os parâmetros gráficos iniciais.
> oldpar <- par(no.readonly = TRUE) > par(pty="s", mfrow=c(1,2)) > contour(x, y, f) > contour(x, y, f, nlevels=15, add=TRUE) > fa <- (f-t(f))/2 > contour(x, y, fa, nlevels=15) > par(oldpar)
Fazendo um gráfico de imagem
> oldpar <- par(no.readonly = TRUE) > par(pty="s", mfrow=c(1,2)) > image(x, y, f) > image(x, y, fa) > par(oldpar)
E apagando objetos novamente antes de prosseguir.
> objects()
[1] "f" "fa" "mds" "mm" "oldpar" "x" "y"
> rm(x, y, f, fa)
Para encerrar esta sessão vejamos mais algumas funcionalidades do R. O R pode fazer operação com
complexos, note que 1i denota o número complexo i.
> th <- seq(-pi, pi, len=100) > z <- exp(1i*th)
Plotando a sequência de complexos com a parte imaginária versus a real obtém-se círculo. Suponha que desejamos amostrar pontos dentro do círculo de raio unitário. Uma forma simples de fazer isto é tomar números complexos com parte real e imaginária padrão. Depois mapeamos qualquer externo ao círculo no seu recíproco:
> par(pty="s") > plot(z, type="l") > w <- rnorm(100) + rnorm(100)*1i > w <- ifelse(Mod(w) > 1, 1/w, w)
Desta forma todos os pontos estão dentro do círculo unitário, mas a distribuição não é uniforme. Um segundo método usa a distribuição 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", ylab="y") > lines(z) > > w <- sqrt(runif(100))*exp(2*pi*runif(100)*1i) > plot(w, xlim=c(-1,1), ylim=c(-1,1), pch="+", xlab="x", ylab="y") > lines(z)
Apagamos novamente os objetos …
> rm(th, w, z)
…e saímos do R.
q()