Uma primeira sessao com o R

Esta é uma primeira sessão com o R visando dar aos participantes uma idéia 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.

[1] -1.1584296 2.0265046 -0.9943093 1.2948114 0.9499178

[1] -1.158 2.027 -0.994 1.295 0.950

[1] 1.2299913 -0.9408620 0.9974967 -0.7250306 0.4638293

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ém 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ã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 idéias 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.

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 especidicadas 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. As argumentos e básicos são 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

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.

  > 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))

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.

[1] -1.1584296 -0.9943093 0.9499178 1.2948114 2.0265046

[1] -0.9943093 0.9499178 1.2948114

[1] -1.2 -1.0 0.9 1.3 2.0

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.

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.

[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  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

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.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

  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

Gráficos de resíduos são produzidos com plot(). Como a função produz 4 gráficos dividiremos a tela gráfica,

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))

   [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ã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(1, 20, c("linear simples", "ponderada", "loess", "verdadeira"),
  +     lty = c(1, 2, 3, 1), lwd = c(1, 1, 1, 2))

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 voce 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")

Lendo dados como um "data-frame"e inspecionando seu conteúdo. 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")
> 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
  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
  080    4  20   780
  081    5   1   890
  082    5   2   840
  083    5   3   780
  084    5   4   810
  085    5   5   760
  086    5   6   810
  087    5   7   790
  088    5   8   810
  089    5   9   820
  090    5  10   850
  091    5  11   870
  092    5  12   870
  093    5  13   810
  094    5  14   740
  095    5  15   810
  096    5  16   940
  097    5  17   950
  098    5  18   800
  099    5  19   810
  100    5  20   870

Devemos definir Expt e Run como fatores 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)

> 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.5  1.1053 0.363209
  Expt         4  94514 23628.5  4.3781 0.003071 **
  Residuals   76 410166  5396.9
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

   [1] "coefficients"  "residuals"     "effects"       "rank"          "fitted.values"
   [6] "assign"        "qr"            "df.residual"   "contrasts"     "xlevels"
  [11] "call"          "terms"         "model"

    (Intercept)          Run2          Run3          Run4          Run5          Run6
   9.506000e+02 -5.200000e+01 -2.800000e+01  6.000000e+00 -7.600000e+01 -1.040000e+02
           Run7          Run8          Run9         Run10         Run11         Run12
  -1.000000e+02 -4.000000e+01 -1.000000e+01 -3.800000e+01  4.000000e+00 -1.754382e-13
          Run13         Run14         Run15         Run16         Run17         Run18
  -3.600000e+01 -9.400000e+01 -6.000000e+01 -6.600000e+01 -6.000000e+00 -3.800000e+01
          Run19         Run20         Expt2         Expt3         Expt4         Expt5
  -5.000000e+01 -4.400000e+01 -5.300000e+01 -6.400000e+01 -8.850000e+01 -7.750000e+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 segunda a parametrização utilizada. Por default a parametrização é feita tomando o primeiro grupo como referência.

(Intercept) Expt2 Expt3 Expt4 Expt5
909.0 -53.0 -64.0 -88.5 -77.5

> mds <- tapply(Speed, Expt, mean)
> mds

1 2 3 4 5
909.0 856.0 845.0 820.5 831.5

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

> (mds[2] - mds[1])/2

> (2 * mds[3] - mds[1] - mds[2])/6

> (3 * mds[4] - mds[1] - mds[2] - mds[3])/12

> (4 * mds[5] - mds[1] - mds[2] - mds[3] - mds[4])/20

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

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

    [,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

    [,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

                  .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.287978e-17 -0.5345225  2.164914e-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 estivar claro experimente para cada uma destas examinar a matrix 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)

> detach()
> rm(fm, fm0)

No próximo exemplo x é um vetor de 50 valores igualmente espaçados no intervalo [-pi pi]. y idem. O objeto f é uma matrix quadrada com linhas e colunas indexadas por x e y respectivamente com os valores da função cos(y)∕(1 + x²).

  > 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 e definindo a região gráfica como quadrada e fazemos um mapa de contorno de f. Depois adicionamos mais linhas para melhor vizualização. fa é a “parte assimétrica” e t() é transposição. Ao final e 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)

  > oldpar <- par(no.readonly = TRUE)
  > par(pty = "s", mfrow = c(1, 2))
  > image(x, y, f)
  > image(x, y, fa)
  > par(oldpar)

[1] "f" "fa" "mds" "mm" "oldpar" "x" "y"

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((0+1i) * th)

Plotando complexos significa parte imaginária versus real Isto deve ser um 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. 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ã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

1 Uma primeira sessão com o R