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

apêndice 3: gráficos tipo pizza

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/

CRAN: The Comprehensive R Archive Network

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

# letras minúsculas e maiúsculas importam: x e' diferente de X
X <- x+x     

# para ver os valores de uma variável:
x

## [1] 10

## [1] 20

# 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)
round(y,3)

## [1] 16.487

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

# note que, enquanto os elementos de uma função são envolvidos por parênteses, (...), os elementos de um vetor o são por chaves, [...]
# 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

6. funções

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: valor da função para x=-1, a=3 e b=0.1
cauchy(-1,3,0.1)

## [1] 0.001988194

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

plot de uma função

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

# variável x: sequência de -10 a 13 com separação de 0.1
x = seq(-10,13,0.1)
y = cauchy(x,3,1)
# preparando um plot com 1 linha e 2 colunas:
par(mfrow=c(1,2))
# a função plot aceita muitos parâmetros
plot(x,y)
plot(x,y,type='l')

# para fazer uma figura sozinha:
par(mfrow=c(1,1))
# incluindo mais parâmetros:
plot(x,cauchy(x,3,1),type='l',xlab='x',ylab='f(x)',col='salmon',lwd=4)
# incluindo uma linha vertical em x=3:
abline(v=3,lwd=2)
# incluindo uma legenda:
legend('topleft',legend='função de Cauchy',text.col='blue',  bty='n',cex=1.5)

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

col: cor dos pontos ou linhas

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

funções pré-definidas

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

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

reprodutibilidade (importante!):

reprodutibilidade - 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 uniforme 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

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. 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
integrate(f, 0, 1)$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

# I é uma variável tipo lista
# o valor numérico da integral é dado por
I$value

## [1] 2.506628

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

## [1] 1

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

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

# note que esta tabela contém tanto dados reais quanto categóricos

# 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"

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"

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"

# 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. Estatística descritiva de um vetor

O objetivo da estatística descritiva é permitir um conhecimento melhor dos dados, em particular sua distribuição.

Vamos exemplificar com a massa estelar das galáxias da tabela.

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.

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

Estatísticas de posição: média e mediana

media = mean(Mass_star)
media

## [1] 10.89012

mediana = median(Mass_star)

par(mfrow=c(1,1))
hist(Mass_star,main="distribuição de massa estelar",breaks=20,xlab=TeX('$log (M_*/M_{sun})$'), ylab= 'frequência',col='salmon')
abline(v = media,col='blue',lwd=3)
abline(v = mediana,col='green',lwd=3)

Estatísticas de largura: variância, desvio padrão, MAD, largura gaussiana equivalente
($\sigma_G = 0.7413 (q_{75} - q_{25})$)

# variância:
var(Mass_star)

## [1] 0.2705706

# desvio padrão:
sigma = sd(Mass_star)
sigma

## [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 = mad(Mass_star)
MAD

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

# visualização
par(mfrow=c(1,1))
hist(Mass_star,main="distribuição de massa estelar",breaks=20,xlab=TeX('$log (M_*/M_{sun})$'), ylab= 'frequência',col='salmon')
abline(v = median(Mass_star),col='green',lwd=4)
segments(mediana-sigma, 5, mediana+sigma, 5, col= 'black',lwd=4)
segments(mediana-MAD, 4, mediana+MAD, 4, col= 'magenta',lwd=4)
segments(mediana-sig_G, 3, mediana+sig_G, 3, col= 'blue',lwd=4)
legend('topleft',legend=c('sigma','MAD','sig_G'),text.col=c('black','magenta','blue'),  bty='n')

Estatísticas de forma: curtose e assimetria (kurtosis e skewness)

curtose:

assimetria:

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

2. Box Plots

úteis para a estatística descritiva de amostras univariadas

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.

A função summary() é muito útil para se obter estas estatísticas:

summary(Mass_star)

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

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

Para uma distribuição gaussiana, $\Delta$ corresponde a aproximadamente 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

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

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

##    Length     Class      Mode 
##        86 character character

plot(Age,Mass_star,col=factor(Morph_class),xlab='log (age/yr)',ylab='log (Mstar/Msun)', pch=16, cex=2, cex.lab=1.5) 
legend("topleft",
       legend = levels(factor(Morph_class)),
       pch = 19,
       col = factor(levels(factor(Morph_class))))

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
xl=TeX('$M_r$')
xl

##    LaTeX: $M_r$ 
## plotmath: M[r]

#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})$'))

4. 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')

5. 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)",phi=10,theta=70,colvar=Mass_star)

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

## Loading required namespace: mgcv

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

III- Exercícios desta semana

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. Quantos são maiores que 3?

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

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

Leia o arquivo novais.txt e analise os dados da coluna logt, que dá a idade média (ponderada em luz) da população estelar de uma galáxia. Exiba a distribuição desses dados tanto usando histograma quanto box-plot. Calcule os valores das estatísticas descritivas dessa quantidade. Existe evidência de “outliers”?

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/

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

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

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 Mr e logMstar: em R o ajuste linear é dado pela função lm(), de linear model:

fit <- lm(dados$logMstar ~ dados$Mr)
# resultado:
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.

Apêndice 3: gráficos de 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))

aula1

lsj

março, 2025

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

apêndice 3: gráficos tipo pizza

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:

6. funções

definindo funções

plot de uma função

funções pré-definidas

7. loops

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

9. métodos numéricos

10. pacotes e bibliotecas

11. leitura de um arquivo de dados:

II. Estatística descritiva e gráficos

1. Estatística descritiva de um vetor

2. Box Plots

3. Scatter Plots

4. Plots de densidade

5. Plots 3-D

III- Exercícios desta semana

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

Apêndice 2: Figuras com ggplot

Apêndice 3: gráficos de pizza