Oficina de R

Semana Acadêmica - 2022

Instalação R e RStudio

• Deve-se instalar primeiro o R, cujo download pode ser feito em https://cloud.r-project.org/.

• Em seguida, pode-se instalar o RStudio, que pode ser definido como um editor integrado ao R, com várias funcionalidades. O download do RStudio pode ser feito em https://www.rstudio.com/products/rstudio/download/#download.

Introdução ao R

Definição de variáveis

• Crie um vetor chamado stock.prices com os seguintes pontos, referente a preços de ações: 23, 27, 23, 21, 34, 37.

stock.prices <- c(23, 27, 21, 34, 37)

Note que aqui temos um vetor numérico.

str(stock.prices)

 num [1:5] 23 27 21 34 37

• Crie um vetor chamado dias.semana com os dias da semana de segunda a sexta.

dias.semana <- c("segunda", "terça", "quarta", "quinta", "sexta")

Note que aqui temos um vetor de caracteres.

str(dias.semana)

 chr [1:5] "segunda" "terça" "quarta" "quinta" "sexta"

• Atribua nomes aos elementos do vetor stock.prices relacionado seus valores ao dia da semana.

names(stock.prices) <- dias.semana

Note que agora o vetor stock.prices tem um atributo a mais.

str(stock.prices)

 Named num [1:5] 23 27 21 34 37
 - attr(*, "names")= chr [1:5] "segunda" "terça" "quarta" "quinta" ...

attributes(stock.prices)

$names
[1] "segunda" "terça"   "quarta"  "quinta"  "sexta"  

Note que c, names, str e attributes são funções nativas do R. Há inúmeras funções, que conhecemos à medida que aprofundamos o estudo do software.

• Calcule a média do vetor stock.prices.

mean(stock.prices)

[1] 28.4

Podemos guardar essa média em um objeto, se for de interesse.

mean_sp <- mean(stock.prices)

Elementos lógicos: TRUE e FALSE

• Podemos testar condições e definir ações dependendo dessas condições.

• Por exemplo, podemos testar se 2 > 3.

2 > 3

[1] FALSE

• Crie um vetor chamado over.23 consistindo de valores lógicos (TRUE ou FALSE) que correspondam aos dias em que o preço das ações ultrapassa 23.

over.23 <- stock.prices > 23
print(over.23)

segunda   terça  quarta  quinta   sexta 
  FALSE    TRUE   FALSE    TRUE    TRUE 

• Podemos usar os valores lógicos para filtrar os valores alocados dentro de um objeto.

• Use o vetor over.23 para filtrar o vetor stock.prices, de modo a retornar apenas os elementos em que o preço foi superior a 23.

stock.prices[over.23]

 terça quinta  sexta 
    27     34     37 

• Elementos lógicos, TRUE ou FALSE, são entendidos, respectivamente, como 0 e 1 em operações numéricas.

2 * TRUE

[1] 2

2 * FALSE

[1] 0

• Suponha uma função f(x) tal que f(x) = 1, se x > 0 e f(x) = -1, se x \le 0. Podemos escrever isso da seguinte forma.

x <- 2 # exemplo
1 * (x > 0) -1 * (x <= 0)

[1] 1

• Adiante, veremos como escrever funções no R.

Funções no R - Introdução

• No R, as funções possuem argumentos. Valores default podem ser inseridos.

• Quando o nome da função é conhecido, podemos usar a função help para entender melhor a lista de argumentos de uma dada função.

• Por exemplo, vamos ver em detalhes a função matrix.

help(matrix)

• Podemos criar uma matriz preenchendo os valores para os argumentos. A função entenderá os valores dados aos seus argumentos pela ordem ou pelo nome do argumento.

mat_ex_1 <- matrix(1:12, nrow = 4, ncol = 3, byrow = TRUE)
print(mat_ex_1)

     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6
[3,]    7    8    9
[4,]   10   11   12

• De modo mais compacto

mat_ex_2 <- matrix(1:12, 4, 3, TRUE)
print(mat_ex_2)

     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6
[3,]    7    8    9
[4,]   10   11   12

Estrutura de dados

• Já vimos como escrever vetores no R.

Matrizes

• Uma matriz pode ser definida por concatenação de vetores.

• Crie um vetor denominado st.pr.sq que é igual ao quadrado de stock.prices. Em seguida crie uma matriz, denominada st_pr_mat, que concatena esses dois vetores.

st.pr.sq <- stock.prices^2
st_pr_mat <- cbind(stock.prices, st.pr.sq)
print(st_pr_mat)

        stock.prices st.pr.sq
segunda           23      529
terça             27      729
quarta            21      441
quinta            34     1156
sexta             37     1369

• Uma matriz aceita apenas um tipo de objeto. Por exemplo, vamos tentar concatenar os vetores stock.prices e dias.semana

st_ds_mat <- cbind(stock.prices, dias.semana)
print(st_ds_mat)

        stock.prices dias.semana
segunda "23"         "segunda"  
terça   "27"         "terça"    
quarta  "21"         "quarta"   
quinta  "34"         "quinta"   
sexta   "37"         "sexta"    

str(st_ds_mat)

 chr [1:5, 1:2] "23" "27" "21" "34" "37" "segunda" "terça" "quarta" ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:5] "segunda" "terça" "quarta" "quinta" ...
  ..$ : chr [1:2] "stock.prices" "dias.semana"

• Notamos que o numérico é transformado em carácter.

• A fim de concatenar vetores de mesmo tamanho, com estruturas diferentes, podemos usar o data frame, conforme veremos adiante.

• Podemos também usar a função matrix para criar matrizes, conforme vimos anteriormente.

mat_1 <- matrix(1:12, nrow = 4, ncol = 3, byrow = TRUE)
print(mat_1)

     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6
[3,]    7    8    9
[4,]   10   11   12

• Note que podemos escolher como os elementos são alocados na matriz: se por linhas ou por colunas.

mat_2 <- matrix(1:12, nrow = 4, ncol = 3, byrow = FALSE)
print(mat_2)

     [,1] [,2] [,3]
[1,]    1    5    9
[2,]    2    6   10
[3,]    3    7   11
[4,]    4    8   12

• Aqui usamos uma forma de criar sequências no R.

1:12

 [1]  1  2  3  4  5  6  7  8  9 10 11 12

seq(from = 1, to = 12, by = 1)

 [1]  1  2  3  4  5  6  7  8  9 10 11 12

• Podemos incluir nomes paras linhas e colunas de uma matriz.

colnames(mat_1) <- c("col.1", "col.2", "col.3")
rownames(mat_1) <- c("row.1", "row.2", "row.3", "row.4")
print(mat_1)

      col.1 col.2 col.3
row.1     1     2     3
row.2     4     5     6
row.3     7     8     9
row.4    10    11    12

Indexação de vetores e matrizes

• Podemos extrair um ou mais elementos de um vetor. Por exemplo, vamos pegar o segundo elemento do vetor stock.prices.

stock.prices[2]

terça 
   27 

• Vamos agora pegar o segundo e o quarto elementos do vetor stock.prices.

stock.prices[c(2, 4)]

 terça quinta 
    27     34 

• Caso o vetor seja nomeado, podemos usar os nomes para extrair valores. Por exemplo, vamos extrair os valores de terça e quinta do vetor stock.prices.

stock.prices[c("terça", "quinta")]

 terça quinta 
    27     34 

• Para matrizes, a indexação acontece de forma semelhante, mas tanto nas linhas como nas colunas. Por exemplo, vamos pegar o elemento (3, 2) da matriz st_pr_mat.

st_pr_mat[3, 2]

[1] 441

• Caso a matriz seja nomeada, podemos usar os nomes para extrair valores. Por exemplo, vamos extrair os valores das linhas terça e quinta da coluna stock.prices da matriz st_pr_mat.

st_pr_mat[c("terça", "quinta"), "stock.prices"]

 terça quinta 
    27     34 

• Caso queiramos a segunda coluna da matriz st_pr_mat, basta fazer

st_pr_mat[, 2]

segunda   terça  quarta  quinta   sexta 
    529     729     441    1156    1369 

• Note que se o produto da extração de elementos de uma matriz for um vetor, o objeto resultante perde a estrutura de matriz.

print(st_pr_mat[, 2])

segunda   terça  quarta  quinta   sexta 
    529     729     441    1156    1369 

str(st_pr_mat[, 2])

 Named num [1:5] 529 729 441 1156 1369
 - attr(*, "names")= chr [1:5] "segunda" "terça" "quarta" "quinta" ...

Note que isso é diferente de

print(as.matrix(st_pr_mat[, 2]))

        [,1]
segunda  529
terça    729
quarta   441
quinta  1156
sexta   1369

str(as.matrix(st_pr_mat[, 2]))

 num [1:5, 1] 529 729 441 1156 1369
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:5] "segunda" "terça" "quarta" "quinta" ...
  ..$ : NULL

• Caso contrário, a estrutura é mantida.

str(st_pr_mat[c(2,4), c(1, 2)])

 num [1:2, 1:2] 27 34 729 1156
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:2] "terça" "quinta"
  ..$ : chr [1:2] "stock.prices" "st.pr.sq"

Função which

• Essa função retorna a posição tal que a condição seja verdadeira.

• Por exemplo, podemos encontrar a posição dos elementos de stock.prices superiores a 23.

which(stock.prices > 23)

 terça quinta  sexta 
     2      4      5 

• Ou ainda,

which(over.23)

 terça quinta  sexta 
     2      4      5 

• Note que os comandos a seguir retornam os mesmos valores.

stock.prices[over.23]

 terça quinta  sexta 
    27     34     37 

stock.prices[which(over.23)]

 terça quinta  sexta 
    27     34     37 

• A função which pode ser bem interessante para matrizes. Vamos considerar a matriz mat_1 como exemplo.

mat_1

      col.1 col.2 col.3
row.1     1     2     3
row.2     4     5     6
row.3     7     8     9
row.4    10    11    12

• Vamos pegar as posições dos elementos maiores do que 5 na matriz mat_1

pos_ge5 <- which(mat_1 > 5, arr.ind = TRUE)
print(pos_ge5)

      row col
row.3   3   1
row.4   4   1
row.3   3   2
row.4   4   2
row.2   2   3
row.3   3   3
row.4   4   3

• Com o objeto pos_ge5 podemos pegar os elementos maiores do que 5 da matriz mat_1.

Operações com vetores e matrizes e reciclagem no R

• Podemos fazer operações algébricas com vetores e matrizes.

• Vamos considerar os seguintes vetores e suas operações, elemento a elemento.

vec_1 <- 1:5
vec_2 <- seq(from = 3, by = 2, length = 5)

print(vec_1)

[1] 1 2 3 4 5

print(vec_2)

[1]  3  5  7  9 11

# adição de vetores
vec_1 + vec_2

[1]  4  7 10 13 16

# subtração de vetores
vec_1 - vec_2

[1] -2 -3 -4 -5 -6

# multiplicação de vetores
vec_1 * vec_2

[1]  3 10 21 36 55

# divisão de vetores
vec_1 / vec_2

[1] 0.3333333 0.4000000 0.4285714 0.4444444 0.4545455

# exponenciação
vec_1^0.3

[1] 1.000000 1.231144 1.390389 1.515717 1.620657

• Podemos calcular o produto interno entre dois vetores.

vec_1 %*% vec_2

     [,1]
[1,]  125

• O mesmo resultado poderia ser obtido através de

sum(vec_1 * vec_2)

[1] 125

• Note que no primeiro caso, o resultado tem estrutura matricial. No segundo caso, temos um escalar.

• Note também que o R entende que vec_1 %*% vec_2 é um produto de duas matrizes. vec_1 é entendido como uma matriz 5 \times 1, enquanto que vec_2 é entendido como uma matriz 1 \times 5.

• Se estivermos interessados no produto externo entre esses dois vetores, devemos transformá-los nas matrizes correspondentes.

mat_1_aux <- as.matrix(vec_1)
mat_2_aux <- as.matrix(vec_2)

print(mat_1_aux)

     [,1]
[1,]    1
[2,]    2
[3,]    3
[4,]    4
[5,]    5

print(mat_2_aux)

     [,1]
[1,]    3
[2,]    5
[3,]    7
[4,]    9
[5,]   11

mat_1_aux %*% t(mat_2_aux)

     [,1] [,2] [,3] [,4] [,5]
[1,]    3    5    7    9   11
[2,]    6   10   14   18   22
[3,]    9   15   21   27   33
[4,]   12   20   28   36   44
[5,]   15   25   35   45   55

• Note que o comando as.matrix aplicado a um vetor, retorna um vetor coluna. No exemplo, uma matriz 5 \times 1.

• A função t retorna a transposta de uma matriz.

• Podemos fazer operações semelhantes com matrizes. Vamos considerar, inicialmente, as operações, elemento a elemento.

print(mat_1)

      col.1 col.2 col.3
row.1     1     2     3
row.2     4     5     6
row.3     7     8     9
row.4    10    11    12

print(mat_2)

     [,1] [,2] [,3]
[1,]    1    5    9
[2,]    2    6   10
[3,]    3    7   11
[4,]    4    8   12

# adição de matrizes
mat_1 + mat_2

      col.1 col.2 col.3
row.1     2     7    12
row.2     6    11    16
row.3    10    15    20
row.4    14    19    24

# subtração de matrizes
mat_1 - mat_2

      col.1 col.2 col.3
row.1     0    -3    -6
row.2     2    -1    -4
row.3     4     1    -2
row.4     6     3     0

# multiplicação de matrizes
mat_1 * mat_2

      col.1 col.2 col.3
row.1     1    10    27
row.2     8    30    60
row.3    21    56    99
row.4    40    88   144

# divisão de matrizes
mat_1 / mat_2

         col.1     col.2     col.3
row.1 1.000000 0.4000000 0.3333333
row.2 2.000000 0.8333333 0.6000000
row.3 2.333333 1.1428571 0.8181818
row.4 2.500000 1.3750000 1.0000000

# exponenciação
mat_1^0.3

         col.1    col.2    col.3
row.1 1.000000 1.231144 1.390389
row.2 1.515717 1.620657 1.711770
row.3 1.792790 1.866066 1.933182
row.4 1.995262 2.053136 2.107436

• Podemos fazer o produto matricial.

mat_3 <- mat_1 %*% t(mat_2)
print(mat_3)

      [,1] [,2] [,3] [,4]
row.1   38   44   50   56
row.2   83   98  113  128
row.3  128  152  176  200
row.4  173  206  239  272

• Podemos obter a inversa de uma matriz.

mat_4 <- cbind(c(23, 2), c(5, 43))
print(mat_4)

     [,1] [,2]
[1,]   23    5
[2,]    2   43

solve(mat_4)

             [,1]         [,2]
[1,]  0.043922370 -0.005107252
[2,] -0.002042901  0.023493361

• Note que mat_4^(-1) fornece o inverso de cada elemento e não a inversa da matriz.

mat_4^(-1)

           [,1]       [,2]
[1,] 0.04347826 0.20000000
[2,] 0.50000000 0.02325581

Reciclagem no R

• Um ponto bastante importante no R é a reciclagem de elementos.

• Por exemplo, se desejamos multiplicar um vetor por um escalar, podemos fazer:

print(vec_1)

[1] 1 2 3 4 5

a <- -2
a * vec_1

[1]  -2  -4  -6  -8 -10

• Se tivermos dois vetores de tamanhos diferentes e fizermos a multiplicação entre eles, teremos

vec_3 <- c(-1, -2, -3)

vec_1 * vec_3

Warning in vec_1 * vec_3: longer object length is not a multiple of shorter
object length

[1]  -1  -4  -9  -4 -10

• Note que a multiplicação é feita, mas recebemos um warning. Note também que o quarto (quinto) elemento de vec_1 é multiplicado pelo primeiro (segundo) elemento de vec_3.

• No caso de matrizes, no produto, elemento a elemento, de um vetor com uma matriz haverá a reciclagem do vetor por coluna e não por linha.

• Por exemplo, vamos analisar a situação abaixo.

vec_4 <- 1:4
mat_1

      col.1 col.2 col.3
row.1     1     2     3
row.2     4     5     6
row.3     7     8     9
row.4    10    11    12

vec_4 * mat_1

      col.1 col.2 col.3
row.1     1     2     3
row.2     8    10    12
row.3    21    24    27
row.4    40    44    48

• Agora, vamos considerar a multiplicação de um vetor que não seja múltiplo da matriz.

vec_1

[1] 1 2 3 4 5

vec_1 * mat_1

Warning in vec_1 * mat_1: longer object length is not a multiple of shorter
object length

      col.1 col.2 col.3
row.1     1    10    12
row.2     8     5    30
row.3    21    16     9
row.4    40    33    24

• Note que a operação é realizada, mas obtemos um warning.

• A reciclagem entre matrizes não é permitida.

as.matrix(vec_4)
mat_1
as.matrix(vec_4) * mat_1

[1] "Error in as.matrix(vec_4) * mat_1 : non-conformable arrays"

• Se quiséssemos escrever o produto elemento a elemento na forma matricial, precisaríamos fazer.

mat_5 <- matrix(
  data = vec_4,
  nrow = nrow(mat_1),
  ncol = ncol(mat_1),
  byrow = FALSE
)

mat_5

     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    2    2    2
[3,]    3    3    3
[4,]    4    4    4

mat_1

      col.1 col.2 col.3
row.1     1     2     3
row.2     4     5     6
row.3     7     8     9
row.4    10    11    12

mat_5 * mat_1

      col.1 col.2 col.3
row.1     1     2     3
row.2     8    10    12
row.3    21    24    27
row.4    40    44    48

Data frames

• Uma estrutura bastante útil para banco de dados é o data frame.

• Um data frame permite a combinação de variáveis com estruturas diferentes.

• Vamos considerar como exemplo o data frame mtcars, presente no banco de dados nativo do R.

• Para carregá-lo, devemos fazer

data(mtcars)

help(mtcars)

• Vamos criar um objeto igual a esse banco de dados e explorar algumas de suas características.

df1 <- mtcars

str(df1)

'data.frame':   32 obs. of  11 variables:
 $ mpg : num  21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
 $ cyl : num  6 6 4 6 8 6 8 4 4 6 ...
 $ disp: num  160 160 108 258 360 ...
 $ hp  : num  110 110 93 110 175 105 245 62 95 123 ...
 $ drat: num  3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
 $ wt  : num  2.62 2.88 2.32 3.21 3.44 ...
 $ qsec: num  16.5 17 18.6 19.4 17 ...
 $ vs  : num  0 0 1 1 0 1 0 1 1 1 ...
 $ am  : num  1 1 1 0 0 0 0 0 0 0 ...
 $ gear: num  4 4 4 3 3 3 3 4 4 4 ...
 $ carb: num  4 4 1 1 2 1 4 2 2 4 ...

Indexação de data frames

• Podemos acessar cada variável usando o operador \$.

mpg <- df1$mpg
head(mpg)

[1] 21.0 21.0 22.8 21.4 18.7 18.1

• Podemos indexar pelo nome ou pela posição, exatamente como fazemos com matrizes.

• Por exemplo, construa um data.frame que reúne apenas as variáveis am, gear e carb.

df2 <- df1[, c("am", "gear", "carb")]
head(df2)

                  am gear carb
Mazda RX4          1    4    4
Mazda RX4 Wag      1    4    4
Datsun 710         1    4    1
Hornet 4 Drive     0    3    1
Hornet Sportabout  0    3    2
Valiant            0    3    1

• Construa um data.frame que reúne todos os carros de 6 cilindros presentes em df4. (O número de cilindros está na coluna cyl).

df3 <- df1[df1$cyl == 6, ]
print(df3)

                mpg cyl  disp  hp drat    wt  qsec vs am gear carb
Mazda RX4      21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag  21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
Hornet 4 Drive 21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
Valiant        18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
Merc 280       19.2   6 167.6 123 3.92 3.440 18.30  1  0    4    4
Merc 280C      17.8   6 167.6 123 3.92 3.440 18.90  1  0    4    4
Ferrari Dino   19.7   6 145.0 175 3.62 2.770 15.50  0  1    5    6

• Note que aqui usamos um operador lógico, df1$cyl == 6.

• Modifique o objeto df1 a fim de incluir uma coluna denominada performance, que é dada pela divisão de hp por wt.

df4 <- df1
df4$performance <- df1$hp / df1$wt

head(df4)

                   mpg cyl disp  hp drat    wt  qsec vs am gear carb
Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1
                  performance
Mazda RX4            41.98473
Mazda RX4 Wag        38.26087
Datsun 710           40.08621
Hornet 4 Drive       34.21462
Hornet Sportabout    50.87209
Valiant              30.34682

Atribuição de característica

• Exemplo mtcars. Podemos classificar o carro quanto ao consumo de combustível.

• crie um objeto denominado fuel.con que assume os seguintes valores:

  • "baixo" se a variável mpg é menor ou igual a 15;
  • "medio" se a variável mpg é maior do que 15 e menor ou igual a 23;
  • "alto" se a variável mpg é maior do que 23.

data(mtcars)

fuel.con <- rep(NA, nrow(mtcars))

fuel.con[mtcars$mpg <= 15] <- "baixo"
fuel.con[mtcars$mpg > 15 & mtcars$mpg <= 23] <- "medio"
fuel.con[mtcars$mpg > 23] <- "alto"

• Agora, vamos definir um novo objeto, denominado, mymtcars, que concatena o data.frame mtcars e a variável fuel.con.

mymtcars <- data.frame(mtcars, fuel.con)
head(mymtcars)

                   mpg cyl disp  hp drat    wt  qsec vs am gear carb fuel.con
Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4    medio
Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4    medio
Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1    medio
Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1    medio
Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2    medio
Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1    medio

• Podemos considerar a variável fuel.con como fator (factor), o que é interessante se quisermos considerá-la como variável categórica.

fuel.con <- as.factor(fuel.con)
head(fuel.con)

[1] medio medio medio medio medio medio
Levels: alto baixo medio

mymtcars <- data.frame(mtcars, fuel.con)
head(mymtcars)

                   mpg cyl disp  hp drat    wt  qsec vs am gear carb fuel.con
Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4    medio
Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4    medio
Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1    medio
Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1    medio
Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2    medio
Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1    medio

Importação e exportação de dados

• Podemos importar dados em vários formatos no R: .csv, .xlsx, .dta, .ods, etc.

• Vamos ver como importar os dados no formato .csv. A importação dos outros formatos ocorre de modo semelhante, mas com funções específicas.

• Com exemplo, vamos considerar os seguintes dados: dados_exemplo_1.csv.

• Fazendo o download desse arquivo, devemos definir o diretório de trabalho e executar o seguinte comando:

dados_ex1 <- read.table(file = "dados_exemplo_1.csv", header = TRUE)

head(dados_ex1)

     lwage exper educ expersq
1 1.131402     2   11       4
2 1.175573    22   12     484
3 1.098612     2   11       4
4 1.791759    44    8    1936
5 1.667707     7   12      49
6 2.169054     9   16      81

str(dados_ex1)

'data.frame':   526 obs. of  4 variables:
 $ lwage  : num  1.13 1.18 1.1 1.79 1.67 ...
 $ exper  : int  2 22 2 44 7 9 15 5 26 22 ...
 $ educ   : int  11 12 11 8 12 16 18 12 12 17 ...
 $ expersq: int  4 484 4 1936 49 81 225 25 676 484 ...

• Se quisermos exportar um banco de dados no R, podemos fazer:

educsq <- dados_ex1$educ^2
dados_ex2 <- data.frame(dados_ex1, educsq)

write.table(
  x = dados_ex2,
  file = "dados_exemplo_2.csv",
  row.names = FALSE
)

Exemplo: Modelos de regressão

• Um modelo de regressão linear pode ser descrito como: y_i = \vec{x}_i^{\intercal} \beta + \epsilon_i, \ \qquad i = 1,\dots, n, onde (\vec{x}_i^{\intercal} ,y_i), \ i = 1, \dots, n é uma sequência iid.

Além disso,

• y_i é a variável dependente;

• \epsilon_i é o termo de erro, com média zero, não observável;

• x_i^{\intercal} = (1, x_{i1}, x_{i2}, \dots, x_{iK}) é o vetor linha 1 \times (K+1) de regressores (uma constante e K covariáveis) e

• \beta = (\beta_0, \beta_1, \dots, \beta_K)^{\intercal} é o vetor coluna (K+1) \times 1 de coeficientes.

• O estimador de Mínimos Quadrados Ordinários (MQO) para esse problema é dado por: \begin{aligned} \vec{\hat \beta} &= (X^{\intercal} X)^{-1} X^{\intercal} y; \\ \hat \sigma^2 &= \frac{1}{n - K - 1} \sum_{i = 1}^{n} \hat \epsilon_i^2 \end{aligned} aqui X é uma matriz n \times K tal que \vec{x}_i^{\intercal} é a i-ésima linha de X e y é um vetor coluna com os elementos y_1, y_2, \dots, y_n.

• De modo mais claro,

X = \begin{bmatrix} 1 & x_{11} & x_{12} & \dots & x_{1K}\\ 1 & x_{21} & x_{22} & \dots & x_{2K}\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n1} & x_{n2} & \dots & x_{nK}\\ \end{bmatrix} é a matriz dos regressores contendo as n observações de cada covariável.

Regressão simples - formato somatório

• No caso de regressão com apenas uma variável explicativa, regressão simples, temos: y_i = \beta_0 + \beta_1 x_i + \epsilon_i, \ \qquad i = 1,\dots, n,

• Nesse caso, os estimadores de MQO para \beta_0 e \beta_1 são: \begin{align} \hat \beta_0 &= \bar{y} - \hat \beta_1 \bar{x} \\ \hat \beta_1 &= \frac{ \sum_{i = 1}^{n} (x_i - \bar{x}) y_i }{ \sum_{i = 1}^{n} (x_i - \bar{x})^2 } \end{align} onde \bar{x} é a média da variável x.

• Considere o banco de dados mtcars, disponível no R.

• Estime por MQO a seguinte equação: mpg_i = \beta_0 + \beta_1 cyl_i + \epsilon_i, \ i = 1, \dots, n.

data(mtcars)

x <- mtcars$cyl
y <- mtcars$mpg

xb <- mean(x)
yb <- mean(y)

b1_hat <- sum((x - xb) * y) / sum((x - xb)^2)
print(b1_hat)

[1] -2.87579

b0_hat <- yb - b1_hat * xb
print(b0_hat)

[1] 37.88458

Regressão simples - formato matricial

• Agora, vamos considerar o formato matricial do estimador MQO. \vec{\hat \beta} = (X^{\intercal} X)^{-1} X^{\intercal} y.

• Vamos repetir o exemplo anterior, mas agora considerando o formato matricial.

x_mat <- cbind(cte = 1, x)
y_vec <- as.matrix(y)

beta_hat <- solve(t(x_mat) %*% x_mat) %*% t(x_mat) %*% y_vec

print(beta_hat)

        [,1]
cte 37.88458
x   -2.87579

Regressão simples - função pronta

• Agora, a ideia e utilizar a função lm, disponível na distribuição básica do R.

reg_1 <- lm(mpg ~ 1 + cyl, data = mtcars)
reg_1

Call:
lm(formula = mpg ~ 1 + cyl, data = mtcars)

Coefficients:
(Intercept)          cyl  
     37.885       -2.876  

Listas

• Uma estrutura de dados muito útil no R é a lista.

• Uma lista é capaz de concatenar objetos de estruturas variadas e comprimentos diferentes.

• Vejamos um exemplo.

l1 <- list(
  vetor = vec_1,
  matriz = mat_1,
  character = dias.semana
)

print(l1)

$vetor
[1] 1 2 3 4 5

$matriz
      col.1 col.2 col.3
row.1     1     2     3
row.2     4     5     6
row.3     7     8     9
row.4    10    11    12

$character
[1] "segunda" "terça"   "quarta"  "quinta"  "sexta"  

str(l1)

List of 3
 $ vetor    : int [1:5] 1 2 3 4 5
 $ matriz   : int [1:4, 1:3] 1 4 7 10 2 5 8 11 3 6 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:4] "row.1" "row.2" "row.3" "row.4"
  .. ..$ : chr [1:3] "col.1" "col.2" "col.3"
 $ character: chr [1:5] "segunda" "terça" "quarta" "quinta" ...

Indexação em listas

• A indexação nas listas apresenta algumas diferenças. Um elemento i da lista pode ser obtido por [[i]]. Por exemplo,

l1[[2]]

      col.1 col.2 col.3
row.1     1     2     3
row.2     4     5     6
row.3     7     8     9
row.4    10    11    12

• Caso a lista seja nomeada, podemos extrair seu elemento usando $ e seu nome.

l1$matriz

      col.1 col.2 col.3
row.1     1     2     3
row.2     4     5     6
row.3     7     8     9
row.4    10    11    12

• Note que se usarmos [] teremos uma lista menor.

l1[1]

$vetor
[1] 1 2 3 4 5

str(l1[1])

List of 1
 $ vetor: int [1:5] 1 2 3 4 5

l1[1]$vetor

[1] 1 2 3 4 5

• Note que o objeto reg_1, oriundo da função lm é, de fato, uma lista.

str(reg_1)

List of 12
 $ coefficients : Named num [1:2] 37.88 -2.88
  ..- attr(*, "names")= chr [1:2] "(Intercept)" "cyl"
 $ residuals    : Named num [1:32] 0.37 0.37 -3.58 0.77 3.82 ...
  ..- attr(*, "names")= chr [1:32] "Mazda RX4" "Mazda RX4 Wag" "Datsun 710" "Hornet 4 Drive" ...
 $ effects      : Named num [1:32] -113.65 -28.6 -3.7 0.71 3.82 ...
  ..- attr(*, "names")= chr [1:32] "(Intercept)" "cyl" "" "" ...
 $ rank         : int 2
 $ fitted.values: Named num [1:32] 20.6 20.6 26.4 20.6 14.9 ...
  ..- attr(*, "names")= chr [1:32] "Mazda RX4" "Mazda RX4 Wag" "Datsun 710" "Hornet 4 Drive" ...
 $ assign       : int [1:2] 0 1
 $ qr           :List of 5
  ..$ qr   : num [1:32, 1:2] -5.657 0.177 0.177 0.177 0.177 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:32] "Mazda RX4" "Mazda RX4 Wag" "Datsun 710" "Hornet 4 Drive" ...
  .. .. ..$ : chr [1:2] "(Intercept)" "cyl"
  .. ..- attr(*, "assign")= int [1:2] 0 1
  ..$ qraux: num [1:2] 1.18 1.02
  ..$ pivot: int [1:2] 1 2
  ..$ tol  : num 1e-07
  ..$ rank : int 2
  ..- attr(*, "class")= chr "qr"
 $ df.residual  : int 30
 $ xlevels      : Named list()
 $ call         : language lm(formula = mpg ~ 1 + cyl, data = mtcars)
 $ terms        :Classes 'terms', 'formula'  language mpg ~ 1 + cyl
  .. ..- attr(*, "variables")= language list(mpg, cyl)
  .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:2] "mpg" "cyl"
  .. .. .. ..$ : chr "cyl"
  .. ..- attr(*, "term.labels")= chr "cyl"
  .. ..- attr(*, "order")= int 1
  .. ..- attr(*, "intercept")= int 1
  .. ..- attr(*, "response")= int 1
  .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. ..- attr(*, "predvars")= language list(mpg, cyl)
  .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
  .. .. ..- attr(*, "names")= chr [1:2] "mpg" "cyl"
 $ model        :'data.frame':  32 obs. of  2 variables:
  ..$ mpg: num [1:32] 21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
  ..$ cyl: num [1:32] 6 6 4 6 8 6 8 4 4 6 ...
  ..- attr(*, "terms")=Classes 'terms', 'formula'  language mpg ~ 1 + cyl
  .. .. ..- attr(*, "variables")= language list(mpg, cyl)
  .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : chr [1:2] "mpg" "cyl"
  .. .. .. .. ..$ : chr "cyl"
  .. .. ..- attr(*, "term.labels")= chr "cyl"
  .. .. ..- attr(*, "order")= int 1
  .. .. ..- attr(*, "intercept")= int 1
  .. .. ..- attr(*, "response")= int 1
  .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. .. ..- attr(*, "predvars")= language list(mpg, cyl)
  .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
  .. .. .. ..- attr(*, "names")= chr [1:2] "mpg" "cyl"
 - attr(*, "class")= chr "lm"

• Desse modo, podemos, por exemplo, recuperar os valores estimados para os coeficientes.

reg_1$coef

(Intercept)         cyl 
   37.88458    -2.87579 

• Note que o início do nome é suficiente para se extrair um elemento da lista, desde que não haja ambiguidade.

Criando nossa própria função

• Para criar funções no R, basta usar a função function.

• A função function recebe os argumentos desejados e retorna o que for definido pelo usuário.

• Vamos, por exemplo, criar uma função que tem como argumentos:

  • uma matriz de regressores,
  • um vetor com a variável dependente. E retorna os valores estimados para os coeficientes de uma regressão linear estimada por MQO.

f_mqo <- function(
  x, # matriz dos regressores
  y # vetor referente à variável dependente
) {
  
  beta_hat <- solve(t(x) %*% x) %*% t(x) %*% as.matrix(y)

  return(beta_hat)
}

• Vamos também incluir em nossa função os valores preditos \vec{\hat{y}} = X \vec{\hat{\beta}} e os resíduos \vec{\hat{\epsilon}} = \vec{y} - \vec{\hat{\beta}}.

• É interessante a função retornar uma lista, pois, nessa caso, podemos misturar objetos de diferentes estruturas.

f_mqo <- function(
  x, # matriz dos regressores
  y # vetor referente à variável dependente
) {
  
  beta_hat <- solve(t(x) %*% x) %*% t(x) %*% as.matrix(y)
  y_hat <- x %*% beta_hat
  e_hat <- y - y_hat

  return(
    list(
      coeficientes = beta_hat,
      y_predito = as.vector(y_hat),
      residuos = as.vector(e_hat)
    )
  )
}

• Vamos testar nossa função para o exemplo anterior.

data(mtcars)

x <- mtcars$cyl
y <- mtcars$mpg

x_mat <- cbind(cte = 1, x)
y_vec <- as.matrix(y)

reg_2 <- f_mqo(x = x_mat, y = y_vec)

str(reg_2)

List of 3
 $ coeficientes: num [1:2, 1] 37.88 -2.88
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:2] "cte" "x"
  .. ..$ : NULL
 $ y_predito   : num [1:32] 20.6 20.6 26.4 20.6 14.9 ...
 $ residuos    : num [1:32] 0.37 0.37 -3.58 0.77 3.82 ...

• Comparando com a função lm, temos

str(reg_1[1:5])

List of 5
 $ coefficients : Named num [1:2] 37.88 -2.88
  ..- attr(*, "names")= chr [1:2] "(Intercept)" "cyl"
 $ residuals    : Named num [1:32] 0.37 0.37 -3.58 0.77 3.82 ...
  ..- attr(*, "names")= chr [1:32] "Mazda RX4" "Mazda RX4 Wag" "Datsun 710" "Hornet 4 Drive" ...
 $ effects      : Named num [1:32] -113.65 -28.6 -3.7 0.71 3.82 ...
  ..- attr(*, "names")= chr [1:32] "(Intercept)" "cyl" "" "" ...
 $ rank         : int 2
 $ fitted.values: Named num [1:32] 20.6 20.6 26.4 20.6 14.9 ...
  ..- attr(*, "names")= chr [1:32] "Mazda RX4" "Mazda RX4 Wag" "Datsun 710" "Hornet 4 Drive" ...

• As funções criadas procuram dentro de si os objetos, caso eles não sejam encontrados, a procura é feita no workspace. Vejamos um exemplo.

f1 <- function(x) {
  x + y
}

y <- 2

f1(3)

[1] 5

f2 <- function(x, y) {
  x + y
}

f2(3, 4)

[1] 7

Exemplo: Integração Numérica - Regra Trapezoidal

A computação numérica de uma integral I = \int_{a}^{b} h(x) dx

pode ser feita pela seguinte regra trapezoidal: \hat{I} = \frac{1}{2} \sum_{i=1}^{n-1}(x_{i+1} - x_{i})(h(x_i) + h(x_{i+1}));

onde os x_i^{'}s constituem uma partição ordenada do intervalo [a,b].

• Calcule a seguinte integral usando a função integrate do R. I = \int_{-1}^{1} (3x^2 + 2x) dx

### Numerical Integration 

f1 <- function(x) 3 * x^2 + 2 * x

integrate(f1, -1, 1) # a

2 with absolute error < 2.5e-14

• Construa uma função que implemente a regra trapezoidal.

• Calcule a integral acima usando a função criada no item anterior.

f.trap <- function(xg, h) {  # b
  x1 <- sort(xg)
  ng <- length(xg)
  res <- 0.5 * (x1[-1] - x1[-ng]) * (h(x1[-1]) + h(x1[-ng]))

  return(sum(res))
}

int.ex1 <- f.trap(
  xg = seq(from = -1, to = 1, length = 10000),
  h = f1
)

print(int.ex1)

[1] 2

Loops

Loop for

• Nessa estrutura, a sequência na qual o índice do loop toma valores está fixa. Por exemplo, se quisermos fazer uma soma acumulada sem usar a função cumsum.

x <- 1
for(i in 2:10){
  x[i] <- x[i-1] + i
}
print(x)

 [1]  1  3  6 10 15 21 28 36 45 55

print(cumsum(1:10))

 [1]  1  3  6 10 15 21 28 36 45 55

Exemplo: estimador não-paramétrico de densidade (kernel)

• O estimador de densidade não-paramétrico kernel (KDE) pode ser definido da seguinte forma: para uma dada amostra x_1,x_2,\cdots,x_n, a densidade estimada num ponto de interesse x_0 é definida por: \hat f(x_0) = \frac{1}{nh}\sum_{i=1}^{n} k\left(\frac{x_i - x_0}{h} \right), onde k é uma função densidade simétrica ao redor de zero e h > 0 uma bandwidth, que controla o grau de suavidade do estimador.

• Considerando k(x) = 0.75 (1 - x^2) I_{[-1, 1]}(x), implemente uma função que tem como argumentos x, x_0 e h e retorna \hat f(x_0).

# k Epanechnikov #

k.epa <- function(x) {
  0.75 * (1 - x^2) * (x > -1) * (x < 1)
}

# Densidade Kernel #

f.kde.0 = function(x, x0, h, k = k.epa) {
  # x: vetor com a amostra
  # x0: valor real (constante)
  # h: janela (bandwidth)
  # k: densidade kernel (default: dnorm)
  # retorna: f.hat.x0
  x1 <- (x - x0) / h
  n <- length(x)
  f.hat.x0 <- sum(k(x1)) / (n * h)
  return(f.hat.x0)
}

• Agora considere uma amostra de tamanho n = 300 da distribuição N(0, 1). Estime \hat f(x_0), x_0 = 1 usando o estimador KDE. Considere h = 2.7 \hat \sigma n^{-1/5}.

set.seed(1234)

n <- 300
x <- rnorm(n, 0, 1)
h <- 2.7 * sd(x) * n^(-1 / 5)
x0 <- 1
f.hat <- f.kde.0(x, x0, h)

print(f.hat)

[1] 0.2222478

print(dnorm(x0))

[1] 0.2419707

• Pode ser interessante estimar a densidade em um grid de pontos. Isso pode ser feito com um loop for.

• Implemente uma função que tem como argumentos x, x_0 e h e retorna \hat f(x_0), onde x_0 agora é um vetor de pontos de interesse.

# Densidade Kernel #

f.kde = function(x, x0, h, k = k.epa) {
  # x: vetor com a amostra
  # x0: vetor de pontos de interesse
  # h: janela (bandwidth)
  # k: densidade kernel (default: k.epa)
  # retorna: f.hat (vetor de tamanho length(x0))
  ng <- length(x0)
  f.hat <- rep(NA, ng)
  for (i in 1:ng) {
    f.hat[i] <- f.kde.0(x = x, x0 = x0[i], h = h, k = k)
  }
  return(f.hat)
}

• Considerando a amostra gerada no exemplo anterior, estime \hat f(x_0), usando o estimador KDE. Agora x_0 é uma sequência igualmente espaçada de comprimento 100 entre -3 e 3. Considere h = 2.7 \hat \sigma n^{-1/5}.

set.seed(1234)

n <- 300
x <- rnorm(300, 0, 1)
h <- 2.7 * sd(x) * n^(-1 / 5)

x0 <- seq(from = -3, to = 3, length = 100)
f.hat <- f.kde(x, x0, h)

dnorm(x0[20:25])

[1] 0.07226707 0.08068571 0.08975477 0.09947714 0.10984842 0.12085626

f.hat[20:25]

[1] 0.06536703 0.07568038 0.08639090 0.09841405 0.11091243 0.12290779

Loop while

• O loop while no R é usado para repetir um bloco de instruções até que a expressão especificada seja FALSE. O loop while começa com a expressão lógica e, se a expressão for TRUE, as instruções dentro do loop while serão executadas.

cc <- 1 # valor para tornar a condição verdadeira.

while (cc > 0) { # condição fundamental do loop while
  print("Hello, world!")
  cc <- cc - .1 # tornará a condição FALSE, impedindo nova execução do loop.
}

[1] "Hello, world!"
[1] "Hello, world!"
[1] "Hello, world!"
[1] "Hello, world!"
[1] "Hello, world!"
[1] "Hello, world!"
[1] "Hello, world!"
[1] "Hello, world!"
[1] "Hello, world!"
[1] "Hello, world!"
[1] "Hello, world!"

Exemplo: Método de Newton

• Vamos iniciar nosso estudo com o método de Newton-Raphson para obter a raiz de uma equação. Ou seja, queremos encontrar x tal que f(x)=0.

• Seja f uma função diferenciável. A expansão de Taylor nos diz que uma boa aproximação para f(x_1) ao redor de x_0 é dada por: f(x_1) = f(x_0) + f'(x_0)(x_1-x_0) + resto

• Suponha que f(x_1)=0 e que o resto seja negligenciável, então: x_1 = x_0 - \frac{f(x_0)}{f'(x_0)}

• Podemos usar essa ideia para resolver uma equação do tipo f(x) = \alpha. Nesse caso, x_1 = x_0 - \frac{f(x_0)-\alpha}{f'(x_0)}.

• Seja X \sim N(0,1). Dado um valor \alpha \in (0,1), o percentil de nível \alpha é o valor que soluciona a equação F(x_{\alpha}) = \alpha. Encontre o percentil de nível 0.025 para a distribuição N(0,1) usando o método de Newton acima. Note que F(x_{\alpha}) - \alpha = 0.

set.seed(1234)
cc <- 1 # condição inicial. Precisa ser verdadeira.
conta <- 0 # contador
a <- 0.025 # percentil de interesse

x0 <- 0  # valor inicial
eps <- 0.0001 # nível de tolerância escolhido
while (cc > eps) {
  x1 <- x0 - (pnorm(x0) - a) / dnorm(x0)
  cc <- (x1 - x0)^2 # métrica de convergência
  x0 <- x1
  conta <- conta + 1
}

# valor estimado
print(x1)

[1] -1.959947

# valor verdadeiro
print(qnorm(a))  # valor função pronta do R.

[1] -1.959964

Loop repeat

Seja X o número de ensaios independentes de Bernoulli até que o primeiro sucesso ocorra. Suponha que a probabilidade de sucesso seja p=0.3. Qual é o número esperado de repetições independentes até que o primeiro sucesso ocorra?

set.seed(1234) # semente
nr <- 5000 # numero de repetições para se ter uma boa estimativa da média
y <- rep(NA, nr)
for (j in 1:nr) {
  i <- 1
  x <- NA # guarda as realizações da v.a. X
  repeat{
    x[i] <- sample(c("s", "f"), size = 1, prob = c(0.3, 0.7))
    if (x[i] == "s") break # condição e comando para interromper o repeat
    i <- i + 1
  }
  y[j] <- length(x)
}
mean(y) # valor verdadeiro 1/0.3

[1] 3.3238

Um pouco mais sobre comandos Lógicos

• Em diversas situações desejamos executar parte do código apenas se alguma condição for atendida.

if

O operador if executa uma expressão, caso uma condição de comprimento um seja verdadeira. Por exemplo, se temos um dia de sol podemos programar ir à praia.

set.seed(1234)

dia <- sample(c("chuva", "sol"), size = 1)

if (dia == "sol") {
  print("Hoje é dia de praia!")
}

[1] "Hoje é dia de praia!"

print(dia)

[1] "sol"

else

Podemos usar if para criar várias condições em sequência, mas se desejamos apenas uma situação e o caso contrário, podemos usar o comando else. No exemplo anterior,

set.seed(123)

dia <- sample(c("chuva", "sol"), size = 1)

if (dia == "sol") {
  print("Hoje é dia de praia!")
} else {
  print("Hoje é dia de ficar em casa.")
}

[1] "Hoje é dia de ficar em casa."

print(dia)

[1] "chuva"

ifelse

Uma forma ainda mais compacta, útil em algumas situações, é a função ifelse.

set.seed(123)

dia <- sample(c("chuva", "sol"), size = 1)

ds <- "Hoje é dia de praia!"
dc <- "Hoje é dia de ficar em casa."

ifelse(dia == "sol", yes = ds, no = dc)

[1] "Hoje é dia de ficar em casa."

Exemplo

Um mineiro está perdido dentro de uma mina e se depara com três portas:

• A primeira porta conduz o mineiro à liberdade após três horas de viagem.

• A segunda porta conduz o mineiro de volta ao ponto inicial (em frente às três portas) após cinco horas de viagem.

• A terceira conduz o mineiro de volta ao ponto inicial após sete horas de viagem.

• Suponha que o mineiro tenha problema de memória, de tal modo que todas as vezes escolhe qualquer uma das três portas com igual probabilidade.

• Sendo assim, qual é o número esperado de horas até que o mineiro consiga a liberdade?

• Construa um script para estimar o referido valor esperado. A ideia aqui é usar os comandos if e repeat, além de um loop de for para aproximar o valor esperado.

nr <- 500

y <- rep(NA, nr) # resultado do experimento em cada volta

for (i in seq_len(nr)) {
  tempo_saida <- 0
  repeat {
    porta <- sample(1:3, size = 1)
    if (porta == 1) {
      tempo_saida <- tempo_saida + 3
      break
    }
    if (porta == 2) {
      tempo_saida <- tempo_saida + 5
    }
    if (porta == 3) {
      tempo_saida <- tempo_saida + 7
    }
  }
  y[i] <- tempo_saida
}

resp_estimada <- mean(y)
print(resp_estimada)

[1] 15.146

• A resposta verdadeira é 15 horas.

Gráficos

• No que diz respeito à utilização de gráficos no R, a função mais comum é a plot. Com vários argumentos, é uma função bastante flexível.

• Voltando ao nosso exemplo de densidade, temos:

## Exemplo ##

set.seed(1234)
n <- 500
x <- rchisq(n, 10)

h1 <- 2.7 * sd(x) * n ^ (-1 / 5)
f.ex2 <- f.kde(x = x, x0 = x, h = h1) 
# avaliando nos pontos da amostra

plot(
  x = sort(x),
  y = f.ex2[order(x)],
  type = "l",
  ylim = c(0, 0.12)
)
lines(sort(x), dchisq(sort(x), 10), col = 2)

Analisando o comando

plot(sort(x), f.ex2[order(x)], type = "l")

é importante notar que os valores de x foram ordenados e os valores de f.ex2 foram ligados aos correspondentes valores de x.

• Desse modo, quando escolhemos type = "l" os pontos são ligados na sequência correta. Caso contrário teríamos:

plot(x, f.ex2, type = "l")

Outros argumentos interessantes da função plot

plot(
  x = sort(x), # ordenada
  y = f.ex2[order(x)], # abscissa
  type = "l", # tipo de gráfico. ver `?plot` para mais opções
  col = "black", # cor do que está sendo plotado. Podem ser usados números
  main = "Gráfico densidade Kernel", # título para o gráfico
  sub = "exemplo", # subtítulo para o gráfico
  xlab = "x", # nome do eixo x
  ylab = "densidade", # nome do eixo y
  ylim = c(0, 0.15), # intervalo mostrado do eixo y
  xlim = c(0, 35), # intervalo mostrado do eixo x
)

Função lines

• O comando

lines(sort(x), dchisq(sort(x), 10), col = 2)

nos permite adicionar um gráfico por cima do grafico previamente adicionado. Ou seja, a função lines tem nível mais alto de que a função plot.

• Algumas funções apresentam. Por exemplo, as funções hist que calcula o histograma e a função curve, útil para gráfico de funções:

plot(
  x = sort(x),
  y = f.ex2[order(x)],
  type = "l",
  main = "Estimadores de densidade",
  sub = expression(paste("Densidade verdadeira" - chi^2)),
  xlab = "x",
  ylab = "densidade",
  ylim = c(0, 0.12),
  xlim = c(min(x), max(x))
)

curve(
  dchisq(x, 10), # função que se deseja plotar
  from = 0, # extremo esquerdo do intervalo de interesse dentro do suporte da função
  to = 35, # extremo direito
  add = TRUE, # TRUE significa que o gráfico será plotado por cima da janela gráfica ativa
  col = 2 # cor
)

hist(
  x = x,
  freq = FALSE, # FALSE para densidade. TRUE para frequência absoluta
  add = TRUE,
  col = NULL, # cor do preenchimento. NULL deixa sem cor
  border = 1, # cor da borda do histograma
)

Exemplo

Considere o seguinte fato:

Se X \sim \chi_n^2, Z \sim N(0,1) e X e Z são independentes, então, Y = \frac{Z}{\sqrt{X/n}} tem distribuição com n graus de liberdade. Denota-se Y \sim t_{n}.

1. Usando a função rnorm para simular valores de uma distribuição N(0,1), simule 500 observações da distribuição t_5.

2. Faça um histograma de densidade para a amostra gerada e plot por cima a densidade verdadeira.

set.seed(1234)

nr <- 1000
n <- 5

x <- rchisq(nr, df = 5)
z <- rnorm(nr)

# item 1

y <- z / sqrt(x / n)

# item 2

par(mfrow = c(1, 2))

hist(y, freq = FALSE, ylim = c(0, 0.5))
curve(dt(x, df = n), add = TRUE, col = 2)

# alternativamente

hist(y, freq = FALSE, ylim = c(0, 0.5))
lines(sort(y), dt(sort(y), df = n), col = 2)

Função qplot do pacote ggplot2

• O pacote ggplot2 trouxe novidades em termos de estética e, mais ainda, em termos de sintaxe para gráficos no R.

• Para quem já está acostumado com o comando plot e se interessa pelo ggplot2, uma boa forma de iniciar a transição é através da função qplot.

• A função qplot se assemelha em termos de sintaxe à função plot.

library(ggplot2)

set.seed(1234)
n <- 200
x <- rchisq(n, 10)

f_hat <- f.kde(x = x, x0 = x, h = 2.7 * sd(x) * n^(-1 / 5))

qplot(
  x = sort(x),
  y = f_hat[order(x)],
  data = NULL,
  geom = "line", 
  size = I(1), 
  color = I("black"),
  xlab = "", 
  ylab = "",
  main = "Estimadores de densidade"
) +
geom_line(aes(y = dchisq(sort(x), 10)), color = "red")

Adicionando layers

• Contudo, a função qplot não é genérica como a função plot, não sendo possível integrá-la com comandos como lines, por exemplo. Para isso precisamos adicionar um novo layer, já na sintaxe do pacote ggplot2.

library(ggplot2)

aux_hist <- hist(
  x = x,
  plot = FALSE
)$breaks 

qplot(
  x = x,
  y = ..density..,
  geom = "histogram", 
  breaks = aux_hist,
  color = I("black"),
  fill = I("grey"),
  xlab = "", 
  ylab = "density",
  main = "Estimadores de densidade"
) +
geom_line(aes(y = dchisq(x, 10)), color = "red") +
geom_line(aes(y = f_hat), color = "blue")

Estrutura básica da função ggplot

• A ideia aqui é fornecer uma introdução à estrutura básica para utilização da função ggplot em situações de maior interesse em nossa disciplina.

• Para quem tiver interesse em algo mais detalhado, um blog interessante no Github, de um ex-aluno, pode ser acessado pelo link

• Há um grande número de possibilidades e aplicações das funções contidas no pacote ggplot2. A ideia aqui não é cobrir essas possibilidades, mas sim fornecer uma estrutura introdutória.

• A ideia básica do pacote ggplot2 é trabalhar com camadas (layers). As camadas fundamentais são:

  • DATA: é necessário definir o data.frame com os dados a serem utilizados;

  • AESTHETICS: pode ser entendida como propriedades visuais dos gráficos. É bastante flexível e pode ser modificada à medida que novas camadas sejam incorporadas. Representada pela função aes.

  • GEOMETRIES: é a camada utilizada para representar os pontos dos dados. Usando as propriedades da função aes dentro de cada particular função geom_* permite a representação de variáveis. Cada função geom_* retorna uma camada.

Voltando ao exemplo anterior

-Vamos reproduzir os gráficos vistos acima. Voltando ao nosso exemplo de densidade:

## Exemplo ##

set.seed(1234)
n <- 200
x <- rchisq(n, 10)

f_hat <- f.kde(x = x, x0 = x, h = 2.7 * sd(x) * n^(-1 / 5))

df_ex2 <- data.frame(
  "x" = x,
  "densidade.estimada" = f_hat,
  "densidade.verdadeira" = dchisq(x, df = 10)
)

g_ex2 <- ggplot(
  data = df_ex2,
  aes(x = x, y = densidade.estimada)
) +
geom_line() +
geom_line(
  aes(x = x, y = densidade.verdadeira),
  color = "red"
)

print(g_ex2)

• É interessante notar que na primeira utilização da função geom_line, os elementos definidos na função aes foram herdados da primeira camada. Entretanto, a função aes pode ser utilizada novamente, a fim de modificar as características do gráfico, como é feito na segunda utilização da função geom_line.

Incorporando o histograma

• Vamos agora refazer o gráfico gerado através da função qplot.

g_ex3 <- ggplot(
  data = df_ex2
) +
geom_histogram(
  aes(x = x, y = ..density..),
  breaks = aux_hist,
  color = "black",
  fill = "grey",
  alpha = 0.5 # opacidade
) +
geom_line(
  aes(x = x, y = densidade.estimada),
  color = "blue"
) +
geom_line(
  aes(x = x, y = densidade.verdadeira),
  color = "red"
)

print(g_ex3)

Alguns elementos gráficos

• Agora podemos adicionar alguns elementos ao gráfico:

g_ex4 <- g_ex3 +
labs(
  title = "Estimadores de densidade",
  ylab = "density"
)

print(g_ex4)

Família apply¹

apply():

• Geralmente é aplicado em uma matriz/dataframe de forma a executar uma mesma função em todas as linhas ou colunas daquele objeto.

################################
######       APPLY         #####
################################

# matriz com 20 colunas e 10 linhas
x <- matrix(rnorm(200), ncol = 20)

# média na linha (MARGIN = 1)
media_linha <- apply(x, MARGIN = 1, mean)   
# aplicar a função mean nas linhas

# média na coluna (MARGIN = 2)
media_coluna <- apply(x, MARGIN = 2, mean)  
# aplicar a função mean nas colunas

lapply():

• É uma função que é aplicada em cada elemento de um vetor ou cada nó de uma lista.

• O output é uma lista obrigatoriamente. Recebe como argumentos um vetor/lista e uma função.

################################
######       LAPPLY        #####
################################

# Uma lista qualquer com 7 vetores de temperaturas em cada dia
temp <- list(
  c(3, 7,  9,  6, -1),
  c(6,  9, 12, 13,  5),
  c(4,  8,  3, -1, -3),
  c(1,  4,  7,  2, -2),
  c(5, 7, 9, 4, 2),
  c(-3,  5,  8,  9,  4),
  c(3, 6, 9, 4, 1)
  )

# Temperatura mínima em cada dia. Retorna uma lista.
lapply(temp, min)

[[1]]
[1] -1

[[2]]
[1] 5

[[3]]
[1] -3

[[4]]
[1] -2

[[5]]
[1] 2

[[6]]
[1] -3

[[7]]
[1] 1

# Exemplo com um vetor

lapply(1:10, sqrt)

[[1]]
[1] 1

[[2]]
[1] 1.414214

[[3]]
[1] 1.732051

[[4]]
[1] 2

[[5]]
[1] 2.236068

[[6]]
[1] 2.44949

[[7]]
[1] 2.645751

[[8]]
[1] 2.828427

[[9]]
[1] 3

[[10]]
[1] 3.162278

sapply():

• Similar ao lapply, porém a saída geralmente é simplificada, sendo apenas um vetor.

• Caso sua saída seja mais de um elemento, a saída deixa de ser um vetor e passa a ser uma matriz.

• Recebe como argumentos um vetor/lista e uma função.

################################
######       SAPPLY        #####
################################

# Temperatura mínima em cada dia. Retorna um vetor.
sapply(temp, min)

[1] -1  5 -3 -2  2 -3  1

# Exemplo com um vetor

sapply(1:10, sqrt)

 [1] 1.000000 1.414214 1.732051 2.000000 2.236068 2.449490 2.645751 2.828427
 [9] 3.000000 3.162278

Exemplo

• Sabe-se que se Z_1, \cdots, Z_n é uma sequência iid com distribuição N(0,1), então X = \sum_{i=1}^n Z_i^2 tem distribuição de Qui-quadrado com n graus de liberdade. Denota-se X \sim \chi_n^2.

• Usando a função rnorm para simular valores de uma distribuição N(0,1), simule 500 observações das distribuições \chi_5^2.

• Use as funções apply e replicate.

nr <- 500
n <- 5

z <- replicate(n = nr, expr = rnorm(n = n, 0, 1))

x <- apply(z^2, 2, sum)

hist(x, freq = FALSE)
curve(dchisq(x, df = n), add = TRUE, col = 2)

Notas de rodapé

1. Agradecimento ao prof. Hudson Chaves por ter sedido o material referente a essa parte.