aula de hoje:

I. introdução ao R

II. estatística descritiva e gráficos

III. exercícios desta semana

apêndice 1: outras introduções ao R

apêndice 2: figuras com o ggplot

I. Introdução ao R

Luis Argerich, professor de Data Science na Universidad de Buenos Aires (UBA):

R is very different than other programming languages, you don’t learn R, you just use it.

This is very important and is something that most programmers will find puzzling. You don’t read a book about R control structures, functions and so on and start doing things to learn the language.

You have a problem, you find how to solve it using R and you learn whatever part of the language you need for it. Repeat this until you can do a lot of things with R. You will find the very curious case of people using R for years and being able to do a lot of things with R and they might not claim they know the language.

para instalar o R no computador de sua casa: visite http://cran.r-project.org/

1. R como linguagem

em R todos os comandos são funções que agem sobre objetos de vários tipos e produzem outros objetos

o usuário pode escrever suas próprias funções

* comandos básicos: pode-se rodar R de várias formas:

1 criando um ambiente para rodar R em um terminal

2 usando rstudio (este documento foi feito com rstudio)

3 rodando um script em um terminal

* comando para criar um ambiente para rodar R em um terminal:

R

* comando para sair do ambiente R:

q()

* comando para entrar no rstudio: digite no terminal

rstudio &

* para rodar um script em R da linha de comando do linux:

Rscript nome_do_script

nome_do_script: arquivo em ascii com comandos em R

* para rodar R com um jupyter notebook:

veja https://hub.gke2.mybinder.org/user/binder-examples-r-qnw5j0x3/notebooks/index.ipynb

* ajuda sobre uma função (no console):

?nome_da_função

exemplo: ajuda sobre a função plot: ?plot

2. R como calculadora:

Inicie uma seção em R digitando R no terminal e dando enter; agora você pode começar a entrar comandos:

2+2

## [1] 4

exp(-2)

## [1] 0.1353353

# pi
pi

## [1] 3.141593

# log na base 10
log10(pi)

## [1] 0.4971499

# logaritmo neperiano
log(pi)

## [1] 1.14473

# log na base 2
log(pi, 2)

## [1] 1.651496

# arredondamento de uma resposta em 3 casas decimais:
round(log(pi),3)

## [1] 1.145

faça ?log para ver as opções da função log

sin(2*pi/3)

## [1] 0.8660254

sqrt(19)

## [1] 4.358899

# expoentes
13^(1/5)

## [1] 1.670278

# atribuindo um valor a uma variável: <- ou =
x <- 10
x = 10
X <- x+x     # x e' diferente de X; nem todos os nomes sao possiveis: T e F são usados para as variáveis lógicas TRUE e FALSE
2*x+X^(1/3)

## [1] 22.71442

y <- x/exp(-x/X)
y

## [1] 16.48721

R: operadores

3. sequências e repetições

# sequências
1:5

## [1] 1 2 3 4 5

# colocando uma sequência em um vetor
y <- 1:10

#sequência de 1 a 77 com incrementos de 5 unidades
s <- seq(1,77,5)
s

##  [1]  1  6 11 16 21 26 31 36 41 46 51 56 61 66 71 76

# repetições:
rep(2,10)

##  [1] 2 2 2 2 2 2 2 2 2 2

4. vetores

# concatenando as componentes com o comando c  (concatenate)
vetor1 <- c(1.,-2.,-7.)
vetor1 <- c(vetor1,3.)
vetor1

## [1]  1 -2 -7  3

# acessando o segundo elemento do vetor
# os índices dos vetores em R começam em 1 (em python eles começam em 0)
vetor1[2]

## [1] -2

# acessando os elementos 2 a 4
vetor1[2:4]

## [1] -2 -7  3

# acessando os elementos > 0
vetor1[vetor1 > 0]

## [1] 1 3

# operadores elemento a elemento: +, -, *, /, ^, etc
vetor2 <- -10+c(12,4,8,11)
vetor2

## [1]  2 -6 -2  1

# R é uma linguagem "vetorizada":
# operações elemento a elemento de cada vetor
vetor1+vetor2

## [1]  3 -8 -9  4

v3 <- vetor1/vetor2^2
v3

## [1]  0.25000000 -0.05555556 -1.75000000  3.00000000

vetor1*vetor2

## [1]  2 12 14  3

# produto escalar: usa-se %*%
vetor1%*%vetor2

##      [,1]
## [1,]   31

# funções podem ser aplicadas a todos os elementos de um vetor de uma única vez:
# criamos um vetor
x<-seq(0,1,0.1)
# função aplicada ao vetor
y <- x*exp(-x)
y

##  [1] 0.00000000 0.09048374 0.16374615 0.22224547 0.26812802 0.30326533
##  [7] 0.32928698 0.34760971 0.35946317 0.36591269 0.36787944

# máximo e mínimo de um vetor
max(vetor1)

## [1] 3

min(vetor1)

## [1] -7

# dimensão do vetor
length(vetor1)

## [1] 4

# soma dos elementos
sum(vetor1)

## [1] -5

# soma cumulativa
cumsum(vetor1)

## [1]  1 -1 -8 -5

#produto dos elementos
prod(vetor1)

## [1] 42

# produto acumulado
cumprod(vetor1)

## [1]  1 -2 14 42

# ordenação dos elementos de um vetor
# vetor1
vetor1

## [1]  1 -2 -7  3

# ordenação dos elementos
sort(vetor1)

## [1] -7 -2  1  3

# índice da ordenação dos elementos
order(vetor1)

## [1] 3 2 1 4

5. matrizes

# vetor com 4 elementos
c(1:3, 8)

## [1] 1 2 3 8

# criando matrizes com matrix()
# Matriz 2x2
matrix(c(1:3, 8), 2)

##      [,1] [,2]
## [1,]    1    3
## [2,]    2    8

# matriz 3x2
A <- matrix(c(1, 1, 1, 2, 4, 20),nrow=3,ncol=2) 
A

##      [,1] [,2]
## [1,]    1    2
## [2,]    1    4
## [3,]    1   20

# acessando os elementos de uma matriz (linha,coluna)
# para acessar o elemento (2,1):
A[2,1]

## [1] 1

# para acessar toda a coluna 1:
A[,1]

## [1] 1 1 1

# para acessar toda a linha 3:
A[3,]

## [1]  1 20

# para acessar os elementos das linhas 2 a 3 e da coluna 2:
A[2:3, 2]

## [1]  4 20

# para acessar todos os elementos maiores que 3
A[A > 3]

## [1]  4 20

# dimensão de uma matriz
dim(A)

## [1] 3 2

operadores matriciais:

A <- matrix(c(1, 1, 2, 4),2) #matriz 2x2
A

##      [,1] [,2]
## [1,]    1    2
## [2,]    1    4

B <- matrix(c(2, 1, 3, 1),2) #matriz 2x2
B

##      [,1] [,2]
## [1,]    2    3
## [2,]    1    1

# multiplicacao de matrizes  %*%
A%*%B

##      [,1] [,2]
## [1,]    4    5
## [2,]    6    7

# matriz transposta
t(A)

##      [,1] [,2]
## [1,]    1    1
## [2,]    2    4

# matriz inversa
solve(B)

##      [,1] [,2]
## [1,]   -1    3
## [2,]    1   -2

# determinante de uma matriz
det(A)

## [1] 2

# teste de sanidade:
# multiplicação de uma matriz por sua inversa (com %*%)
B%*%solve(B)

##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1

arrays multi-dimensionais

podemos criar arrays multidimensionais

por exemplo, a partir de um vetor:

ar=array(c(1,2,3,4,5,6),dim=c(2,2,3))
ar

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

# o vetor preenche o array dimensão por dimensão, com o primeiro índice variando mais rápido que o segundo e este mais rápido que o terceiro

6. funções

uma das vantagens do R é que ele tem um grande número de funções pré-definidas

exemplo 1: simulação de 10 números aleatórios uniformemente distribuidos entre 0 e 7 com a função runif()

reprodutibilidade (importante!):

para garantir que os resultados desta simulação possam ser reproduzidos, deve-se definir uma “semente”:

# usamos a função set.seed()
set.seed(42)

# a função runif() permite gerar números com distribuição uniforme
# simulando 10 números com distribuição unifirme entre 0 e 7:
x = runif(10,min=0,max=7)  
x

##  [1] 6.4036423 6.5595279 2.0029767 5.8131334 4.4922186 3.6336716 5.1561182
##  [8] 0.9426662 4.5989460 4.9354535

# se você gerar agora uma outra sequência, ela será diferente da anterior: verifique!

# simule duas sequências, estabelecendo antes de cada simulação a mesma seed

exemplo 2: regressão linear $y = a+bx$

vamos simular dados e ajustá-los com uma reta

# vamos começar simulando dados de uma reta y = a + b x
a0 = 1
b0 = 0.5

# vamos simular 10 valores de x e y como
npts = 10
x = runif(npts,min=0,max=7)
y = a0 + b0*x + rnorm(npts,mean=0,sd=0.5)
# onde adicionamos ruído gaussiano com média 0 e desvio padrão 0.5 (sd = standard deviation) 

# para fazer a regressão linear (= ajuste da reta), usamos a função lm() (de linear model):
ajuste = lm(y~x)
# resultado
ajuste

## 
## Call:
## lm(formula = y ~ x)
## 
## Coefficients:
## (Intercept)            x  
##      0.4104       0.6230

# coeficientes ajustados:
a1=ajuste$coef[1]
b1=ajuste$coef[2]

# visualização do resultado:
plot(x,y)
abline(a=a0,b=b0)
abline(a=a1,b=b1,col='red')
legend('topleft',legend=c('modelo','ajuste'),col=c('black','red'),  text.col = c('black','red'))

definindo funções

sintaxe: nome = function(parâmetros da função){cálculo da função}

vamos considerar alguns exemplos:

duplica = function(u){2*u}
duplica(2)

## [1] 4

# distribuição de Cauchy: ela tem 2 parâmetros, a e b
cauchy = function(x,a,b){
b/(pi*(b^2+(x-a)^2))
}

# exemplo:
x = seq(-10,13,0.1)
plot(x,cauchy(x,3,1),type='l',ylab='f(x)',col='salmon',lwd=4)
abline(v=3)

7. loops

loops em R são lentos, e sempre que possível deve-se procurar trabalhar com vetorização, isto é, com expressões que envolvam apenas operações com vetores e matrizes

loops em R: for, while, repeat

# para ilustrar vamos considerar dois vetores, que escreveremos usando *vetorização*:
x = c(1,2,3,4,5)
# dimensão do vetor x:
n = length(x)
n

## [1] 5

# vetor y usando vetorialização
y = 2*x
y

## [1]  2  4  6  8 10

# agora vetor y, usando loop com for()
y = NULL # para definir um objeto y
# loop:
for(i in 1:n){
y[i]=2*x[i]
}
y

## [1]  2  4  6  8 10

# exemplo de um loop com while()
x = 1
while (x < 6) {
x = 2*x
print(x)
}

## [1] 2
## [1] 4
## [1] 8

# exemplo de um loop usando repeat() 
x = 1
repeat {
x = 2*x
print(x)
if (x == 8){break}
}

## [1] 2
## [1] 4
## [1] 8

# sai-se do loop com 'break'

8. combinação de variáveis em tabelas

já vimos como concatenar variáveis com o comando c(), concatenate. Por exemplo, podemos criar três vetores u, v e z de dimensão 6 como

u = c(3,-4,1,0,1,-1)
v = c(1,-2,5,-1,0,0)
z = c(0,2,2,-2,1,1)

# tamanhos dos vetores:
c(length(u),length(v),length(z))

## [1] 6 6 6

Se quisermos criar uma tabela (matriz) tendo como colunas u, v e z usamos cbind()
(column bind):

tab1 = cbind(u,v,z)
tab1

##       u  v  z
## [1,]  3  1  0
## [2,] -4 -2  2
## [3,]  1  5  2
## [4,]  0 -1 -2
## [5,]  1  0  1
## [6,] -1  0  1

# dimensão da tabela/matriz
dim(tab1)

## [1] 6 3

# nomes das colunas
colnames(tab1)

## [1] "u" "v" "z"

# vamos atribuir novos nomes às colunas: note que os nomes são variáveis nominais
colnames(tab1) = c('var_u','var_v', 'var_z')
tab1

##      var_u var_v var_z
## [1,]     3     1     0
## [2,]    -4    -2     2
## [3,]     1     5     2
## [4,]     0    -1    -2
## [5,]     1     0     1
## [6,]    -1     0     1

Podemos também criar tabelas juntando linhas. Por exemplo, juntando por linha os vetores associados a cada coluna de tab1 (u, v e z):

tab2 = rbind(u,v,z) # row bind
tab2

##   [,1] [,2] [,3] [,4] [,5] [,6]
## u    3   -4    1    0    1   -1
## v    1   -2    5   -1    0    0
## z    0    2    2   -2    1    1

# atribuindo nomes às linhas:
rownames(tab2) = c('var_u','var_v','var_z')
tab2

##       [,1] [,2] [,3] [,4] [,5] [,6]
## var_u    3   -4    1    0    1   -1
## var_v    1   -2    5   -1    0    0
## var_z    0    2    2   -2    1    1

pode-se aplicar funções para variáveis em tabelas com apply():

# vamos ilustrar calculando médias de linhas e colunas de uma tabela
x = 1:5
y = 2^x
# vamos criar uma tabela com os x e y acima
tab = cbind(x,y)  # column bind
tab

##      x  y
## [1,] 1  2
## [2,] 2  4
## [3,] 3  8
## [4,] 4 16
## [5,] 5 32

# média das colunas
apply(tab,2,mean)

##    x    y 
##  3.0 12.4

# média das linhas
apply(tab,1,mean)

## [1]  1.5  3.0  5.5 10.0 18.5

# pode-se usar outras funções
apply(tab,2,duplica)

##       x  y
## [1,]  2  4
## [2,]  4  8
## [3,]  6 16
## [4,]  8 32
## [5,] 10 64

9. a função list()

A função list() pode combinar variáveis de tipos e dimensões diferentes. Exemplo: suponha que v1 seja um vetor de dimensão 3, v2 um conjunto de 4 caracteres, v3 uma variável de dimensão 1 e v4 uma matriz 2x2:

v1 <- c(1, 2, 3)
v2 <- c("a", "b", "c", "d")
v3 <- 3
v4 <- matrix(seq(1,4,1),nrow = 2, ncol = 2)

# criando a lista:
Y <- list(x1 = v1, x2 = v2, x3 = v3, x4 = v4)
Y

## $x1
## [1] 1 2 3
## 
## $x2
## [1] "a" "b" "c" "d"
## 
## $x3
## [1] 3
## 
## $x4
##      [,1] [,2]
## [1,]    1    3
## [2,]    2    4

# para acessar um elemento da lista:
Y$x2

## [1] "a" "b" "c" "d"

# para saber o que tem na lista
names(Y)

## [1] "x1" "x2" "x3" "x4"

Uma coisa que torna a função list() importante é que a saída de muitas funções em R é uma lista.

10. A função merge()

Às vezes deseja-se combinar variáveis de duas tabelas diferentes a partir de uma variável/identificador comum às duas tabelas. Para isso podemos usar a função merge():

# vamos criar uma tabela com 50 linhas (objetos) e 5 colunas (parâmetros)
dados = matrix(runif(250),ncol=5,nrow=50)
colnames(dados)=c('a','b','c','d','e')

# vamos criar duas tabelas a partir da tabela 'dados'
# t1: linhas 1 a 30, colunas 1 a 3
t1 <- dados[1:30,1:3]

# t2: linhas 20 a 50 e colunas 1,4,5
t2 <- dados[20:50,c(1,4,5)]

# e agora criamos uma nova tabela com os dados em comum das duas tabelas usando a coluna comum 1 (variável 'a')
t3 <- merge(t1,t2,by = 'a')
t3

##            a           b         c          d         e
## 1  0.1894739 0.537376695 0.1795557 0.83875503 0.1352305
## 2  0.2405447 0.828942131 0.1861022 0.89271859 0.8867536
## 3  0.2712866 0.001380844 0.4518865 0.57463733 0.8692578
## 4  0.5144129 0.002272966 0.6801642 0.25968998 0.6189515
## 5  0.5664884 0.452731573 0.1969945 0.06094975 0.9457392
## 6  0.6756073 0.608937453 0.9348230 0.54201594 0.7148537
## 7  0.6932048 0.612133090 0.1161747 0.54742608 0.9250459
## 8  0.7595443 0.751522563 0.6017662 0.33641913 0.3110496
## 9  0.8281585 0.355665954 0.3170534 0.35335038 0.2050496
## 10 0.8496897 0.535789994 0.5352366 0.45131085 0.5000251
## 11 0.9828172 0.836801559 0.5504941 0.64987584 0.1233006

11. métodos numéricos

Em R existem rotinas e bibliotecas prontas para inúmeros métodos numéricos que usamos em nosso dia a dia como cientistas: vamos ilustrar com a integração numérica de uma função:

f = function(x) 1/sqrt(x)
integrate(f, 0, 1)

## 2 with absolute error < 5.8e-15

# note que a saida da função integrate() é uma lista
# e, para acessar apenas o resultado da integral, pode-se fazer
g = integrate(f, 0, 1)
g$value

## [1] 2

# outro exemplo:
fgauss = function(x) exp(-x^2/2)
I = integrate(fgauss, -Inf, Inf)
I

## 2.506628 with absolute error < 0.00023

# o valor exato dessa integral é sqrt(2*pi)
I$value/sqrt(2*pi)

## [1] 1

12. pacotes e bibliotecas

O R tem uma infinidade de pacotes e bibliotecas especializadas. Para instalá-las use o comando install.packages(‘nome da biblioteca’).

Para usar uma biblioteca no R, antes de rodar o comando desejado, use o comando library(nome da biblioteca). Veremos alguns exemplos mais abaixo.

13. leitura de um arquivo de dados:

Vamos considerar aqui dados sobre a emissão na linha H$\alpha$ e outros parâmetros de uma amostra de galáxias do CALIFA estudada por Novais e Sodré (2019, MNRAS, 482, 2717).

Há muitas formas de se importar arquivos no R, dependendo do formato dos dados e da presença ou não de um “cabeçalho” (header), por exemplo:

dados <- read.table(file="novais.txt", header=TRUE)

# tamanho da tabela: número de linhas e colunas
dim(dados)

## [1] 86 13

# primeiras linhas
head(dados)

##    Galaxy     Mr  umr logt  logZ logMstar Hubble_type Morphological_class
## 1  IC0776 -18.69 1.94 8.34 -0.27     9.59          Sd               Slate
## 2  IC1256 -20.81 2.21 9.08 -0.29    10.72          Sb              Searly
## 3  IC1683 -20.75 2.54 9.31 -0.24    10.76          Sb              Searly
## 4  IC4566 -21.51 2.88 9.30 -0.26    11.01          Sb              Searly
## 5 NGC0001 -21.11 2.24 8.89 -0.35    10.82         Sbc               Slate
## 6 NGC0036 -21.86 2.48 9.30 -0.34    11.22          Sb              Searly
##   Halpha_profile   cr ceHalpha ccHalpha  eps
## 1             CL 1.95     0.55     0.20 0.28
## 2             CL 2.19     0.57     0.22 0.20
## 3             CL 2.59     0.30     0.77 0.28
## 4             EX 2.24     0.63     0.08 0.26
## 5             CE 3.04     0.52     0.40 0.25
## 6             EX 2.49     0.72     0.05 0.25

# nomes das colunas
colnames(dados)

##  [1] "Galaxy"              "Mr"                  "umr"                
##  [4] "logt"                "logZ"                "logMstar"           
##  [7] "Hubble_type"         "Morphological_class" "Halpha_profile"     
## [10] "cr"                  "ceHalpha"            "ccHalpha"           
## [13] "eps"

# significado das colunas:
# galaxy name, absolute magnitude $M_r$, colour u−r, mean age (log(age/yr)), metallicity (Z/Zsun), stellar mass (log(M/Msun)), Hubble type, morphological class, type of the Halpha profile, light concentration, effective concentration, central concentration, and ellipticity

# a classificação da emissão Halfa nesse trabalho considera se o máximo da emissão está no centro (C) ou não (E, de extended) e, no primeiro caso, se a galáxia é 'early-type' (CE) ou 'late-type' (CL) de acordo com sua concentração cr

#para ver a linha 10:
dados[10,]

##     Galaxy     Mr  umr logt  logZ logMstar Hubble_type Morphological_class
## 10 NGC0496 -21.12 2.32 8.74 -0.31    10.84          Sc               Slate
##    Halpha_profile   cr ceHalpha ccHalpha  eps
## 10             CL 2.03     0.47     0.43 0.34

uma forma de acessar as variáveis de um arquivo: nome_do_arquivo$nome_da variável

# magnitude absoluta na banda r do SDSS:
dados$Mr

##  [1] -18.69 -20.81 -20.75 -21.51 -21.11 -21.86 -21.24 -20.75 -20.79 -21.12
## [11] -21.88 -21.12 -21.46 -22.17 -19.17 -22.15 -21.45 -21.20 -21.60 -20.40
## [21] -20.24 -20.56 -21.17 -22.16 -19.32 -19.96 -20.76 -21.71 -21.32 -20.38
## [31] -20.44 -22.68 -21.36 -19.34 -20.28 -20.60 -22.03 -20.72 -19.50 -22.04
## [41] -21.90 -21.96 -22.51 -21.43 -21.14 -21.29 -20.72 -20.02 -22.01 -22.94
## [51] -22.28 -21.60 -23.11 -22.76 -21.61 -22.11 -21.87 -22.35 -22.07 -21.47
## [61] -21.88 -21.16 -21.31 -19.12 -21.11 -21.28 -21.14 -20.29 -18.37 -20.57
## [71] -21.72 -21.91 -19.34 -22.41 -18.73 -22.09 -20.35 -22.66 -22.09 -18.90
## [81] -22.32 -20.60 -22.49 -21.26 -19.92 -20.28

#para a galáxia #3:
dados$Mr[3]

## [1] -20.75

# classificação de Hubble para as galáxias de 5 a 10:
dados$Hubble_type[5:10]

## [1] Sbc Sb  Sb  Sc  Sbc Sc 
## Levels: E S0 Sa Sb Sbc Sc Sd

outra forma de considerar variáveis: dando nomes explicitamente

Nome <- dados[,1]          #nome da galáxia
Mabs_r <- dados[,2]        #magnitude absoluta na banda r
umr <- dados[,3]           #cor u-r
Age <- dados[,4]           #log da idade média ponderada em luz em Gyr
Mass_star <- dados[,6]     #log da massa estelar em M_sun
Hubble <- as.character(dados[,7])        # tipo de Hubble
Morph_class <- dados[,8]   #classe morfológica
Profile <- dados[,9]       #tipo do perfil de Halfa
CC_Ha <- dados[,12]        #concentração central de Halfa

# exemplo: alguns dados da galaxia 10
Nome[10]

## [1] NGC0496
## 86 Levels: IC0776 IC1256 IC1683 IC4566 NGC0001 NGC0036 NGC0171 ... UGC12816

Mabs_r[10]

## [1] -21.12

Hubble[10]

## [1] "Sc"

CC_Ha[10]

## [1] 0.43

O comando attach permite acessar os elementos de uma base de dados diretamente pelos nomes das variáveis:

attach(dados)

## The following object is masked _by_ .GlobalEnv:
## 
##     umr

# nome da terceira galáxia
Galaxy[3]

## [1] IC1683
## 86 Levels: IC0776 IC1256 IC1683 IC4566 NGC0001 NGC0036 NGC0171 ... UGC12816

# log da idade média (em anos) da população estelar da segunda galáxia 
logt[2]

## [1] 9.08

II. Estatística descritiva e gráficos

1. plot de uma função

Vamos ver como plotar uma função em R:

# variável x: sequência de 0 a 10 com separação de 0.1
x = seq(0,10,0.1)
y = exp(x*sin(x)^2)
# preparando um plot com 1 linha e 2 colunas:
par(mfrow=c(1,2))
plot(x,y)
plot(x,y,type='l')

Há vários parâmetros para tornar os gráficos mais informativos:

O comando pch() define o tipo de símbolo que será usado na figura:

símbolos para pch

lwd: largura de linha (o default é 1)

type = ‘l’: plota uma linha juntando os pontos
há outras opções: faça ?plot

log = ‘x’,‘y’,‘xy’: escala logarítmica

cex: tamanho dos símbolos

cex.axis, cex.lab, cex.main: magnificação em relação a cex das anotações dos eixos, dos labels x e y e do main (o título do plot)

lty: tipo de linha

lwd: largura da linha (default = 1)

e muito mais… veja:

https://www.statmethods.net/advgraphs/parameters.html

Alguns exemplos:

par(mfrow=c(2,2))
# variável x: sequência de 0 a 10 com separação de 0.1
x = seq(0,10,0.1)
y = exp(x*sin(x)^2)

# tipo de linha
plot(x,y,xlab='x',ylab='y',type='l',col='chocolate1',col.lab='blue',lty=3)


plot(x,y,xlab='x',ylab='y',type='l',col='chocolate1',col.lab='blue',lty=6)

# aumento da largura da linha e dos labels
plot(x,y,xlab='x',ylab='y',type='l',col='chocolate1',lwd=2,cex.lab=1.5,col.lab='blue')

# para por o eixo y em log + mudando o tipo de linha
plot(x,y,xlab='x',ylab='y',type='h',log='y',col='chocolate1',lwd=2,cex.lab=1.5,col.lab='blue')

2. Estatística descritiva de um vetor

O objetivo da estatística descritiva é permitir que o analista conheça melhor os dados, em particular sua distribuição.

Vamos considerar a distribuição de valores de uma variável. Vamos exemplificar com a massa estelar das galáxias que consta da tabela acima.

Vou aproveitar e introduzir a biblioteca latex2exp, que permite usar latex nos labels de plots:

# labels usando latex
library(latex2exp)
# lembre-se: se uma biblioteca não está instalada rode
# install.packages('nome da biblioteca')
# e depois library(nome da biblioteca)

Para representar a distribuição de variáveis univariadas, as figuras mais comuns são o histograma e o box-plot:

#para plotar 2 figuras uma ao lado da outra:
par(mfrow=c(1,2))

hist(Mass_star,main="distribuição de massa estelar",breaks=20,xlab=TeX('$log (M_*/M_{sun})$'), ylab= 'frequência',col='salmon')

boxplot(Mass_star,col="red",main='distribuição de massa estelar',xlab=TeX('$log (M_*/M_{sun})$'))

# função quantil
quantile(Mass_star)

##     0%    25%    50%    75%   100% 
##  9.590 10.595 11.000 11.270 11.810

O box-plot mostra algumas estatísticas de uma variável: o mínimo, o percentil 25%, a mediana, o percentil 75% e o máximo.

Note que pelos valores da função quantil você pode ter uma ideia sobre a forma da distribuição.

Vamos ver como são definidas as estatísticas mais comuns em R:

Estatísticas de localização: média e mediana \[ \mu = E(x) = \int_{-\infty}^\infty x P(x) dx \] \[ \int_{-\infty}^{x_{med}} P(x) dx = 1/2 = \int_{x_{med}}^\infty P(x) dx \]

mean(Mass_star)

## [1] 10.89012

median(Mass_star)

## [1] 11

Estatísticas de largura: variância, desvio padrão, MAD, largura gaussiana equivalente \[ V = \int (x-\mu)^2 P(x) dx\] \[\sigma = \sqrt{V}\] \[ MAD = \int |x-\mu| P(x) dx\] \[ \sigma_G = 0.7413 (q_{75} - q_{25}) \]

# variância:
var(Mass_star)

## [1] 0.2705706

# desvio padrão:
sd(Mass_star)

## [1] 0.520164

# note que a variância é igual ao quadrado do desvio padrão:
sd(Mass_star)^2

## [1] 0.2705706

# MAD: maximum absolute deviation
mad(Mass_star)

## [1] 0.511497

# distância intraquartil normalizada pela gaussiana:
sig_G = 0.7413*(quantile(Mass_star,0.75,names = FALSE) - quantile(Mass_star,0.25,names = FALSE))
sig_G

## [1] 0.5003775

Compare o desvio padrão, MAD e sig_G: em que situações você acha que estas estatísticas vão ser muito diferentes entre si?

Estatísticas de forma: kurtosis e skewness \[ \Sigma = \int \Bigg(\frac{x-\mu}{\sigma}\Bigg)^3 P(x) dx\] \[ K = \int \Bigg(\frac{x-\mu}{\sigma}\Bigg)^4 P(x)dx\]

# vamos usar uma biblioteca específica:
library(moments)
# lembre-se: se uma biblioteca não está instalada rode
# install.packages('nome da biblioteca')
# e depois library(nome da biblioteca)

kurtosis(Mass_star)

## [1] 2.860328

skewness(Mass_star)

## [1] -0.6454713

Como estes valores se comparam com os de uma gaussiana? Para uma distribuição gaussiana a curtose é 3 e a distorção é zero.

3. Histogramas

Uma das formas mais comuns de representar a distribuição de uma variável é por histogramas.

Vamos ilustrar com a concentração central de Halfa:

hist(CC_Ha)

#definindo o número de bins (breaks) e a cor (col):
# note os labels: main, xlab, ylab
hist(CC_Ha,main="concentração central de Halfa",breaks=20,xlab='C(Halfa)', ylab= 'frequência',col='salmon')

#letras gregas na figura
hist(CC_Ha,main=expression(paste("concentração central de H",alpha)),breaks=20,xlab=expression(paste('C(H',alpha,')')),col='red')

#introduzindo uma legenda
 legend(0.4,10,legend='galáxias do CALIFA', col='black', text.col='springgreen3',cex=1.2)

# use a opção bty='n' para remover a caixa da legenda:
hist(Age,main=expression(paste("idade média ",tau)),breaks=20,xlab=expression(paste('log ',tau,'  (anos)')),ylab='número',col='maroon1')

legend(7.8,8,legend='dados do starlight',  text.col='slateblue',cex=1.2,bty='n')

# procure no Google por "Colors in R" para ver as opções de cores disponíveis.

Com R dá para fazer gráficos muito bonitos! Veja exemplos em

https://www.r-graph-gallery.com/index.html

Vamos começar carregando uma biblioteca. Se ela não estiver instalada execute antes o comando abaixo:

install.packages(“RColorBrewer”)

# carregando a biblioteca
library(RColorBrewer)

# preparando um plot com 2 linhas e 3 colunas:
par(mfrow=c(2,3))

# histogramas com várias paletas de cores
hist(CC_Ha,breaks=10, col=brewer.pal(3,"Set3"),main="Set3 3 colors")
hist(CC_Ha,breaks=3 ,col=brewer.pal(3,"Set2"),main="Set2 3 colors")
hist(CC_Ha,breaks=7, col=brewer.pal(3,"Set1"),main="Set1 3 colors")
hist(CC_Ha,breaks= 2, col=brewer.pal(8,"Set3"),main="Set3 8 colors")
hist(CC_Ha,col=brewer.pal(8,"Greys"),main="Greys 8 colors")
hist(CC_Ha,col=brewer.pal(8,"Greens"),main="Greens 8 colors")

4. Box Plots

Box Plot: úteis para descrever amostras univariadas

O box-plot mostra algumas estatísticas de uma variável: o mínimo, o percentil 25% ($Q_1$), a mediana ($Q_2$), o percentil 75% ($Q_3$) e o máximo.

Um jeito simples de identificar possíveis outliers (numa distribuição gaussiana) é verificando se há dados fora de $\Delta = \pm 1.5 \times$ IQR, onde IQR, o interquartile range, é o intervalo entre os percentis 25% e 75%.

Para uma distribuição gaussiana, $\Delta$ corresponde a 2.7$\sigma$ (onde $\sigma$ é o desvio padrão da gaussiana)

para exemplos de vários tipos de boxplots, veja https://www.r-graph-gallery.com/boxplot.html

A função summary() é muito útil para se ter uma ideia da distribuição dos dados:

summary(Mass_star)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    9.59   10.60   11.00   10.89   11.27   11.81

par(mfrow=c(2,2))
boxplot(Mass_star,col="red",main='stellar mass')
boxplot(Mass_star~Profile,col="red")
boxplot(CC_Ha~Profile,col=heat.colors(4))
boxplot(CC_Ha~Morph_class,col=topo.colors(3),main='Halpha central concentration')

5. Scatter Plots

A principal maneira de se visualizar relações entre duas variáveis é com scatter plots:

# para colocar um plot ao  lado de outro
par(mfrow=c(1,2))    #plot com uma linha e 2 colunas
plot(Mabs_r,Mass_star,xlab='Mr',ylab=' log (Mstar/Msun)',main='um exemplo')

# pch define a fonte dos pontos a serem plotados
# veja como introduzir uma letra grega num label de plot
plot(Mabs_r,Age,xlab='Mr',ylab=expression(paste('log ',tau,'  (yr)')),col='red',pch=20)

Um tipo de figura bem informativa é o scatter plot colorido por tipo

# aqui o plot corresponde a uma única figura
par(mfrow=c(1,1))
# plot colorido de acordo com a classe morfológica
summary(Morph_class)

##    ES0 Searly  Slate 
##     30     23     33

plot(Age,Mass_star,col=Morph_class,xlab='log (age/yr)',ylab='log (Mstar/Msun)', pch=16, cex=2, cex.lab=1.5) 
legend('topleft',legend=c('E/S0','Searly','Slate'),col=c('black','red','green'),  text.col = c('black','red','green'))

Adicionando certa transparência num scatter plot:

plot(Mabs_r, Mass_star, pch = 19, col = rgb(0, 0, 0, 0.15),cex=1.5)

Scatter plots usando hexbin:

# labels usando latex
library(latex2exp)

library(hexbin)
library(RColorBrewer)
# labels nos eixos x e y usando latex
a=hexbin(Mabs_r,Mass_star,xbins=10,xlab=TeX('$M_r$'),ylab=TeX('$log (M_* /M_{sun})$'))
plot(a)

outra forma do mesmo tipo de plot:

rf <- colorRampPalette(rev(brewer.pal(12,'Set3')))
hexbinplot(as.numeric(Mabs_r)~as.numeric(Mass_star),data=dados, colramp=rf,xlab=TeX('$M_r$'),ylab=TeX('$log (M_* /M_{sun})$'))

6. Pizza

vamos considerar os tipos de perfil de Halfa (uma variável categórica): CE (central early), CL (central late) e EX (extended):

par(mfrow=c(2,3))

pie(table(Profile), col=FALSE)

pie(table(Profile), col=rainbow(4))

pie(table(Profile), col=heat.colors(4))

pie(table(Profile), col=terrain.colors(4))

pie(table(Profile), col=topo.colors(4))

pie(table(Profile), col=cm.colors(4))

7. Plots de densidade

# a função smoothscatter está na biblioteca KernSmooth
library(KernSmooth)

## KernSmooth 2.23 loaded
## Copyright M. P. Wand 1997-2009

smoothScatter(Mabs_r,Mass_star)

# a função kde está na biblioteca MASS
library(MASS)
#crio uma matriz de pontos usando KDE: kernel density estimator
res = kde2d(Mabs_r,Mass_star,  h=0.5,n=500)
#ploto usando a função image
image(res, col=grey(13:0/15), xlab='Mr',ylab='M_*',main = 'função image')

8. Plots 3-D

Vamos agora ilustrar algumas formas de se plotar pontos em 3D. Vamos exemplificar com o espaço (Mabs_r, Mass_star, (u-r)):

library("plot3D")
points3D(Mabs_r,Mass_star,umr,xlab = "Mr",ylab = "M*",zlab = "(u-r)")

library("scatterplot3d")
scatterplot3d(Mabs_r,Mass_star,umr,main="3D Scatter Plot",xlab = "Mr",ylab = "M*",zlab = "(u-r)", pch = 16, color="steelblue", angle = 85)

# esta rotina é iterativa! use ela num terminal rodando R
library(car)

## Loading required package: carData

scatter3d(Mabs_r,Mass_star,umr,point.col = "blue", surface=FALSE,xlab = "Mr",ylab = "M*",zlab = "(u-r)")

## Loading required namespace: rgl

## This version of bslib is designed to work with shiny version 1.6.0 or higher.

## Loading required namespace: mgcv

## Warning in par3d(userMatrix = structure(c(1, 0, 0, 0, 0, 0.342020143325668, :
## font family "sans" not found, using "bitmap"

Todas essas funções têm muitas opções. Use ?nome para saber mais sobre cada uma.

III- Exercícios desta semana

parte 1:

Escolha uma semente de números aleatórios igual a seu número USP.

Gere um vetor v com 20 números aleatórios com distribuição uniforme entre 1 e 5.

Calcule algumas estatísticas do vetor v com as funções min(), max(), median(), sd() e sum(). Use ?nome_da_função se preciso.

Calcule $2*exp(-log(Número USP))$.

Construa as matrizes \[ A = \begin{bmatrix} 2 & 4 & 7 \\ 1 & 5 & 8 \\ 3 & 6 & 9 \end{bmatrix} \] e \[ B = \begin{bmatrix} 1 & 3 & 3 \\ 4 & 2 & 2 \\ 3 & 2 & 3 \end{bmatrix} \] e calcule $A^TB$, $A^{-1}B^{-1}$ e $(AB)^{-1}$.

Simule 100 números aleatórios distribuidos uniformemente entre 0 e 1. Use um loop para verificar quantos valores são maiores que 0.5.

Use vetorização para verificar quantos valores são maiores que 0.5.

Crie uma função que calcula a função sigmoide, $1/(1+e^{−x})$, e calcule sua integral entre -5 e +5.

parte 2:

O arquivo berlind2006tab3M20.txt contém propriedades de grupos e aglomerados de galáxias determinados por Berlind et al. (2006), usando dados do SDSS.

Leia o arquivo berlind2006tab3M20.txt. Verifique as variáveis que ele contém.

A riqueza dá o número de galáxias (acima de uma certa magnitude) em cada estrutura. Examine a distribuição da riqueza nessa amostra; calcule estatísticas dessa distribuição.

Verifique visualmente se a riqueza se correlaciona a dispersão de velocidades sigma. Compare a visualização usando coordenadas lineares e logarítmicas para essas quantidades (considere apenas objetos com sigma > 0).

Verifique a relação entre a riqueza e a cor g-r. Como se compara a cor de sistemas ricos com a de sistemas pobres? Defina sistemas ricos e pobres como os com riqueza acima e abaixo da riqueza mediana, respectivamente e calcule a cor mediana nos dois casos.

Apêndice 1: Outras introduções ao R

a introdução do CRAN (The Comprehensive R Archive Network):
https://cran.r-project.org/doc/manuals/R-intro.html

o grupo de Astroestatística da Penn State (https://sites.psu.edu/astrostatistics/) é um dos mais importantes da atualidade - não deixem de olhar

tutorial do Eric Feigelson (Penn State):
https://sites.psu.edu/astrostatistics/files/2021/02/Tutorial_1_Intro_R.html

curso da Carnegie Mellon:
https://www.stat.cmu.edu/~cshalizi/statcomp/14/

CRAN Task Views: informações sobre bibliotecas de vários tipos CRAN task views aim to provide some guidance which packages on CRAN are relevant for tasks related to a certain topic. https://cran.r-project.org/web/views/

ache uma biblioteca ou uma função
it works basically like a filter on google searches, to return only R-relevant content.
https://rseek.org/

novidades no r-bloggers: https://www.r-bloggers.com/

Apêndice 2: Figuras com ggplot

referências: https://www.curso-r.com/material/ggplot/

https://www.cedricscherer.com/2019/08/05/a-ggplot2-tutorial-for-beautiful-plotting-in-r/

A biblioteca ggplot2 implementa o conceito de gráfico como o mapeamento dos dados a partir de atributos estéticos (posição, cor, forma, tamanho) de objetos geométricos (pontos, linhas, barras, caixas). Nesta biblioteca os principais aspectos de um gráfico (dados, sistema de coordenadas, rótulos e anotações) são divididos em camadas, construídas uma a uma na elaboração do gráfico.

Vamos considerar a relação Mr x M*:

library(ggplot2)

ggplot(dados) +
  geom_point(aes(x = Mr, y = logMstar))

Note que:

a primeira camada é dada pela função ggplot() e recebe um data frame;

a segunda camada é dada pela função geom_point(), especificando a forma geométrica utilizada no mapeamento das observações;

as camadas são unidas com um +;

o mapeamento na função geom_point() recebe a função aes(), responsável por descrever como as variáveis serão mapeadas nos aspectos visuais da forma geométrica escolhida, no caso, pontos.

Vamos introduzir outra camada, com o ajuste linear entre as variáveis: em R o ajuste linear é dado pela função lm(), de linear model:

fit <- lm(dados$logMstar ~ dados$Mr)
fit

## 
## Call:
## lm(formula = dados$logMstar ~ dados$Mr)
## 
## Coefficients:
## (Intercept)     dados$Mr  
##      0.7350      -0.4803

# y =  a+bx
a <- fit$coef[1]
b <- fit$coef[2]

incluindo esta reta no plot:

ggplot(dados) +
  geom_point(mapping = aes(x = Mr, y = logMstar)) +
  geom_abline(intercept = a, slope = b, color = "red")

Os atributos estéticos mais usados são:

color=: altera a cor de formas que não têm área (pontos e retas).

fill=: altera a cor de formas com área (barras, caixas, densidades, áreas).

size=: altera o tamanho de formas.

type=: altera o tipo da forma, geralmente usada para pontos.

linetype=: altera o tipo da linha.

Por exemplo, para alterar a cor dos pontos:

ggplot(dados) +
  geom_point(aes(x = Mr, y = logMstar), color = "blue")

Formas geométricas são definidas com a função geom. As mais utilizadas são

geom_point(): gera gráficos de dispersão transformando pares (x,y)

 geom_line - para linhas definidas por pares (x,y).

 geom_abline - para retas definidas por um intercepto e uma inclinação.

```
 geom_hline - para retas horizontais.
```

```
 geom_bar - para barras.
```

```
geom_histogram - para histogramas.
```

```
 geom_boxplot - para boxplots.
```

```
geom_density - para densidades.
```

```
 geom_area - para áreas.
```

Exemplo: histograma

ggplot(dados) + 
  geom_histogram(aes(x = Mr), color = "black", fill = "white")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Exemplo: box plot

note o uso do comando pipe %$>% da biblioteca dplyr

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:car':
## 
##     recode

## The following object is masked from 'package:MASS':
## 
##     select

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

dados %>% 
  filter(Halpha_profile %in% c("CL", "CE", "EX")) %>%
ggplot() + 
  geom_boxplot(aes(x = Halpha_profile, y = logt))

Os gráficos podem ser construídos em camadas: por exemplo, plotando pontos e adicionando uma reta (no exemplo abaixo dentro do ggplot)

ggplot(dados) +
  geom_point(aes(x = Mr, y = logMstar)) +
  geom_smooth(aes(x = Mr, y = logMstar), se = FALSE, method = "lm")

## `geom_smooth()` using formula 'y ~ x'

Podemos economizar código especificando o aes() diretamente na função ggplot():

ggplot(dados, aes( x = Mr, y = logMstar)) +
  geom_point() +
  geom_smooth(se = FALSE, method = "lm")

## `geom_smooth()` using formula 'y ~ x'

Pode-se usar o aes() dentro dos geoms para especificar mapeamentos específicos. Por exemplo, colorir os pontos com outro parâmetro:

ggplot(dados, aes( x = Mr, y = logMstar)) +
  geom_point(aes(color = logt)) +
  geom_smooth(se = FALSE, method = "lm")

## `geom_smooth()` using formula 'y ~ x'

Podemos replicar um gráfico para cada tipo de variável:

dados %>%
  filter(Morphological_class %in% c("Slate", "Searly", "ES0")) %>%
  ggplot() +
  geom_point(aes(x = Mr, y = logMstar)) +
  facet_wrap(~Morphological_class, nrow = 3)

ou (repare no nrow e no ncol):

dados %>%
  filter(Morphological_class %in% c("Slate", "Searly", "ES0")) %>%
  ggplot() +
  geom_point(aes(x = Mr, y = logMstar)) +
  facet_wrap(~Morphological_class, ncol = 3)

Podemos juntar gráficos de tipos diferentes. Aqui vamos usar a biblioteca patchwork:

library(patchwork)

## 
## Attaching package: 'patchwork'

## The following object is masked from 'package:MASS':
## 
##     area

p1 <- dados %>% 
  filter(Morphological_class == "Slate") %>%
  group_by(Mr) %>% 
  summarise(log_Mstar = mean(logMstar, na.rm = TRUE)) %>% 
  ggplot() +
  geom_line(aes(x = Mr, y = log_Mstar))

p2 <- dados %>% 
  mutate(c = ceHalpha - ccHalpha) %>% 
  filter(Morphological_class == "Slate") %>% 
  ggplot() +
  geom_histogram(
    aes(x = c),
    fill = "lightblue", 
    color = "darkblue", 
    binwidth = 0.20
  )

p1 + p2

Estes são alguns exemplos com ggplot. Este pacote tem muita versatilidade e permite fazer gráficos sofisticados.

aula1

lsj

fevereiro, 2023

aula de hoje:

I. introdução ao R

II. estatística descritiva e gráficos

III. exercícios desta semana

apêndice 1: outras introduções ao R

apêndice 2: figuras com o ggplot

I. Introdução ao R

Luis Argerich, professor de Data Science na Universidad de Buenos Aires (UBA):

R is very different than other programming languages, you don’t learn R, you just use it.

This is very important and is something that most programmers will find puzzling. You don’t read a book about R control structures, functions and so on and start doing things to learn the language.

para instalar o R no computador de sua casa: visite http://cran.r-project.org/

1. R como linguagem

* comandos básicos: pode-se rodar R de várias formas:

* comando para criar um ambiente para rodar R em um terminal:

* comando para sair do ambiente R:

* comando para entrar no rstudio: digite no terminal

* para rodar um script em R da linha de comando do linux:

* para rodar R com um jupyter notebook:

* ajuda sobre uma função (no console):

2. R como calculadora:

R: operadores

3. sequências e repetições

4. vetores

5. matrizes

operadores matriciais:

arrays multi-dimensionais

6. funções

definindo funções

7. loops

8. combinação de variáveis em tabelas

9. a função list()

10. A função merge()

11. métodos numéricos

12. pacotes e bibliotecas

13. leitura de um arquivo de dados:

II. Estatística descritiva e gráficos

1. plot de uma função

2. Estatística descritiva de um vetor

3. Histogramas

4. Box Plots

5. Scatter Plots

6. Pizza

7. Plots de densidade

8. Plots 3-D

III- Exercícios desta semana

parte 1:

parte 2:

Apêndice 1: Outras introduções ao R

Apêndice 2: Figuras com ggplot