aula1
lsj
fevereiro, 2024
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] 4exp(-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.145faça ?log para ver as opções da função log
sin(2*pi/3)## [1] 0.8660254sqrt(19)## [1] 4.358899# expoentes
13^(1/5)## [1] 1.670278# atribuindo um valor a uma variável: <- ou =
x <- 10
x = 10
# x e' diferente de X
X <- x+x
# em 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.71442y <- x/exp(-x/X)
y## [1] 16.48721R: 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 24. 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 4v3 <- vetor1/vetor2^2
v3## [1] 0.25000000 -0.05555556 -1.75000000 3.00000000vetor1*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] 3min(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 45. 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 2operadores matriciais:
A <- matrix(c(1, 1, 2, 4),2) #matriz 2x2
A## [,1] [,2]
## [1,] 1 2
## [2,] 1 4B <- 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 16. 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.001988194x = 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))
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:
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:
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 1: 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 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 seedexemplo 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'))
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 6Se 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 1Podemos 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 1pode-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 649. 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:
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# I é uma variável tipo lista
# o valor numérico da integral é dado por
I$value/
# note que o valor exato dessa integral é sqrt(2*pi)
I$value/sqrt(2*pi)## [1] 0.398942310. 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.34uma 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.12Hubble[10]## [1] "Sc"CC_Ha[10]## [1] 0.43O 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.08II. Estatística descritiva e gráficos
1. 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 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.810O 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.89012mediana = 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')
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: 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.860328skewness(Mass_star)## [1] -0.6454713Como 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.81Um 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 characterplot(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-2009smoothScatter(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)")
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: carDatascatter3d(Mabs_r,Mass_star,umr,point.col = "blue", surface=FALSE,xlab = "Mr",ylab = "M*",zlab = "(u-r)")## Loading required namespace: rgl## Loading required namespace: mgcvTodas 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.
- Calcule algumas estatísticas do vetor v com as funções min(), max(), median(), sd() e summary().
Use ?nome_da_função se preciso.
- Calcule \(2 \times e^{-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}\) (usando produtos matriciais).
- 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.
- 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”?
- Usando dados de novais.txt, use scatter plots para examinar a relação entre logt com a cor (u-r) e a massa estelar. Compare (usando medianas) as idades médias e as cores das galáxias com massa estelar abaixo e acima da massa estelar mediana e comente os resultados.
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
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, uniondados %>%
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':
##
## areap1 <- 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))