aula2

aula de hoje:

I. Dados e histogramas

II. Estatística Descritiva

III. Bootstrap

IV. Análise de Correlação

V. Exercícios da Semana

---

---

I. Dados e histogramas

I.1 leitura dos 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

vamos dar nomes para as variáveis 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
Z <- dados[,5]             # log da metalicidade (fração da massa em metais) 
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

I.2 distribuição de variáveis

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

# função summary
summary(Mass_star)

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

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.

as funções quantile() (em seu default) e summary() também apresentam estas quantidades

I.3 histogramas

o histograma é o procedimento mais comum de visualização da distribuição de densidade de uma variável aleatória: ele mostra a frequência da variável em intervalos discretos de valores, geralmente constantes, que denominaremos bins

considere o problema: temos um conjunto de N dados que queremos representar com um histograma: qual é o tamanho “ideal” do bin?

existe um conjunto de regras para determinar a largura do bin: seja \(\sigma_s\) o desvio padrão da amostra e \(\Delta\) a largura do bin; então:

• regra de Scott: \[ \Delta = \frac{3.5 \sigma_s}{N^{1/3}}\]

• regra de Freedman-Diaconis: \[ \Delta = \frac{2(q_{75} - q_{25})}{N^{1/3}}\] onde \(q_{25}\) e \(q_{75}\) são quartis 25% e 75% da distribuição

• regra de Silverman: \[ \Delta = \frac{1.06 \sigma_s}{N^{1/5}} \]

• regra de Knuth: (https://arxiv.org/pdf/physics/0605197.pdf) - método bayesiano assumindo uma verossimilhança multinomial e um prior não-informativo
  métodos bayesianos dão a distribuição de probabilidade (denominada posterior) do número de bins M com contagens \(n_k\) em cada bin: \[ P(M|D) = N \log M + \log \Bigg[ \Gamma \Bigg(\frac{M}{2} \Bigg) \Bigg] - M \log \Bigg[ \Gamma \Bigg(\frac{1}{2} \Bigg) \Bigg] - \] \[ -\log \Bigg[ \Gamma \Bigg(\Gamma(N+\frac{M}{2} \Bigg) \Bigg] + \sum_{k=1}^M \log \Bigg[ \Gamma \Bigg(n_k+\frac{1}{2} \Bigg) \Bigg] \] \(M\) é obtido maximizando-se o posterior

Vamos exemplificar com a distribuição de massa estelar das galáxias da amostra que estamos considerando:

# dados:
par(mfrow = c(1,2))
# vamos usar nesta figura 40 bins
hist(Mass_star,xlab=TeX('$log (M_*/M_{sun})$'), 
col='salmon',main='',breaks=40)

# distribuições normalizadas: histograma com freq = FALSE
# a integral do histograma é igual a 1
hist(Mass_star,xlab=TeX('$log (M_*/M_{sun})$'), 
col='salmon',main='',breaks=40,freq=FALSE)

número “ótimo” de bins de acordo com as várias regras:

# regra de Scott:
sigma = sd(Mass_star)
# tamanho do bin
bin_scott = 3.5*sigma/length(Mass_star)^(1/3)
bin_scott

[1] 0.412454

# número de bins
ldata = max(Mass_star)-min(Mass_star)
n_scott=ldata/bin_scott+1
# para arredondar:
round(n_scott)

[1] 6

# regra de Freedman-Diaconis
# tamanho do bin
bin_fd = 2*(quantile(Mass_star,0.75,names = FALSE) 
- quantile(Mass_star,0.25,names = FALSE))/
  length(Mass_star)^(1/3)
bin_fd

[1] 0.3058447

# número de bins
n_fd = ldata/bin_fd+1
round(n_fd)

[1] 8

# regra de Silverman:
# tamanho do bin
bin_silverman = 1.06*sigma/length(Mass_star)^(1/5)
bin_silverman

[1] 0.2262281

# número de bins
n_silverman = ldata/bin_silverman+1
round(n_silverman)

[1] 11

# regra de Knuth: teste com até 100 bins

# regra de Knuth
# https://arxiv.org/pdf/physics/0605197.pdf
# optBINS finds the optimal number of bins for a 
# one-dimensionaldata set using the posterior 
# probability for the number of bins
# This algorithm uses a brute-force search trying 
# every possible bin number in the given range.  
# This can of course be improved.
# Generalization to multidimensional data sets is straightforward.
## Usage:
#           optBINS = function(data,maxM)
# Where:
#           data is a (1,N) vector of data points
#           maxM is the maximum number of bins to consider
## Ref: K.H. Knuth. 2012. Optimal data-based binning for histograms
# and histogram-based probability density models, Entropy.
optBINS = function(data,maxM){

N = length(data)

# Simply loop through the different numbers of bins
# and compute the posterior probability for each.

logp = rep(0,N)

for(M in 2:maxM){

# Bin the data (equal width bins here)
n = table(cut(data, breaks = M))

soma = 0
for(k in 1:M){
  soma = soma + lgamma(n[k]+0.5)}
  logp[M] = N*log(M) + lgamma(M/2) - lgamma(N+M/2) - M*lgamma(1/2)+ soma
}
return(logp)
}

a = optBINS(Mass_star,100)

Mknuth = which(a == max(a))
Mknuth

[1] 5

# distribuição de probabilidades de M
par(mfrow = c(1,1))
plot(seq(2,100,1),a[2:100],type='l',xlab='M', ylab='log(prob)',col='red',lwd=3)
abline(v=Mknuth,lwd=2)

# número de bins dos vários métodos:
c(round(n_scott),round(n_fd),round(n_silverman),round(Mknuth))

[1]  6  8 11  5

# visualização
par(mfrow = c(2,2))
plot(cut(Mass_star, breaks = n_scott),col='red')
legend('topleft', 'Scott',cex=0.9,bty='n')
plot(cut(Mass_star, breaks = n_fd),col='green')
legend('topleft', 'Freedman-Diaconis',cex=0.9,bty='n')
plot(cut(Mass_star, breaks = n_silverman),col='darkgoldenrod')
legend('topleft', 'Silverman',cex=0.9,bty='n')
plot(cut(Mass_star, breaks = Mknuth),col='blue')
legend('topleft', 'Knuth',cex=0.9,bty='n')

note que a escolha do tamanho dos bins não é óbvia pois, a rigor, o resultado deve depender da forma intrínseca da distribuição, que não se conhece!

---

II. Estatística descritiva de um vetor

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

[1] 11

par(mfrow=c(1,1))
hist(Mass_star,main="distribuição de massa 
estelar",breaks=8,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=8,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.

Para verificar isso, vamos fazer uma simulação:

# vamos gerar 100 amostras de uma gaussiana de média 0 e desvio padrão 1:

# reprodutibilidade
set.seed(123)   

# a função rnorm gera números aleatórios com distribuição gaussiana
x = rnorm(1000)

# distribuição de x:
hist(x)

# curtose e distorção:
kurtosis(x)

[1] 2.925747

skewness(x)

[1] 0.06529391

as diferenças em relação aos valores esperados são devido ao uso de uma amostra finita.

---

III. Bootstrap

bootstrap: útil para estimar incertezas em estatísticas

técnica baseada na reamostragem (com substituição) dos dados

exemplo- estimativa do erro da mediana das massas estelares:

# erro na mediana de Mass_star por bootstrap
x = Mass_star

# vamos fazer nsim simulações de bootstrap
nsim = 100000

# vetor para armazenar o valor mediano da reamostragem:
medianas = rep(NA, nsim)

# reamostragem nsim vezes:
for (i in 1:nsim) {
# Reamostragem com reposição
  xl = sample(x, replace = T) 

# mediana dos dados reamostrados
  medianas[i] = median(xl) 
}

median(medianas)

[1] 11

sd(medianas)

[1] 0.06705462

q=quantile(medianas, probs = c(0.05, 0.5, 0.95))
q

   5%   50%   95% 
10.85 11.00 11.08 

# visualização 
hist(medianas,xlab='mediana de x', col='salmon',
main='distribuição das medianas no bootstrap')
abline(v=q[1],lwd=2)
abline(v=q[2],lwd=2)
abline(v=q[3],lwd=2)
text(10.7,30000,'quantis 0.05, 0.5, 0.95',col='blue')

---

IV. Análise de Correlação

Vamos considerar o uso dos coeficientes de correlação para analisar a relação entre duas variáveis.

Funções úteis:

• cor(): correlação entre dois vetores, todas as colunas de um data.frame ou dois data.frames

• cor.test(): testa a significância de uma correlação; a hipótese nula é que não há correlação

Vamos ilustrar estas funções com dados de parâmetros de galáxias em Novais & Sodré (2018). Será que a idade média (Age) e a metalicidade (Z) das populações estelares das galáxias da amostra estão correlacionadas?

# a correlação entre dois vetores pode ser calculada com a função cor()
# faça ?cor para saber mais sobre esta função

# o parâmetro method permite escolher a 
# estatística de correlação e covariância, onde as 
# opções são "pearson" (default), "spearman" ou 
# "kendall"
# a função cor() dá diretamente o valor do 
# coeficiente de correlação:
cor(Age,Z,method = 'pearson')

[1] 0.4555887

cor(Age,Z,method = 'spearman')

[1] 0.5200067

cor(Age,Z,method = 'kendall')

[1] 0.3568582

plot(Age,Z,xlab='log t (Ganos)',ylab='log Z/Zsol',
main='relação idade x metalicidade',pch=20,
col='salmon')

Note como diferentes estatísticas produzem resultados diferentes para o coeficiente de correlação.

uma função correlata é a cor.test(), que além do coeficiente de correlação, dá outras estatísticas, como o “valor-p”, que é pequeno se a correlação é significativa:

cor.test(Age,Z,method = 'spearman')

Warning in cor.test.default(Age, Z, method = "spearman"): Cannot compute exact
p-value with ties

    Spearman's rank correlation rho

data:  Age and Z
S = 50877, p-value = 2.879e-07
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.5200067 

No caso, existe correlação, ela não é muito forte (os coeficientes de correlação não são muito altos) mas é significativa (valor-p pequeno).

Vamos agora avaliar a correlação entre a magnitude absoluta na banda r e a massa estelar:

cor(Mabs_r,Mass_star,method = 'spearman')

[1] -0.9823558

plot(Mabs_r,Mass_star,xlab='Mr',ylab='log Mstar/Msol',main='',pch=20,col='salmon')

essa (anti)correlação é muito forte!

Pode-se calcular várias correlações ao mesmo tempo:

library(corrplot)

corrplot 0.95 loaded

library(RColorBrewer)
M <-cor(dados[,c(2:6)],method = 'spearman')  
M

                 Mr        umr       logt       logZ   logMstar
Mr        1.0000000 -0.5844048 -0.6181165 -0.4360995 -0.9823558
umr      -0.5844048  1.0000000  0.8135370  0.4571195  0.6258605
logt     -0.6181165  0.8135370  1.0000000  0.5200067  0.6826351
logZ     -0.4360995  0.4571195  0.5200067  1.0000000  0.4785962
logMstar -0.9823558  0.6258605  0.6826351  0.4785962  1.0000000

A função corrplot da biblioteca de mesmo nome permite plotar a significância de várias variáveis ao mesmo tempo:

library(corrplot)
library(RColorBrewer)
M <-cor(dados[,c(2:6,10:13)],method = 'pearson')  
corrplot(M, type="upper", order="hclust",
         col=brewer.pal(n=8, name="RdYlBu"))        

---

V. Exercícios desta semana

1. Leia o arquivo novais.txt e analise os dados da coluna logt (Age), que dá a idade média (ponderada em luz) da população estelar de uma galáxia. Exiba a distribuição desses dados usando box-plot.

2. Calcule os valores das estatísticas descritivas de logt. Comente os resultados.

3. Faça histogramas de logt usando as várias regras acima e comente os resultados.

4. Inicialize o gerador de números aleatórios com seu número USP. Use bootstrap para calcular o erro da mediana de logt.

5. Analise a correlação entre logt e Mass_star usando Pearson e Spearman; comente os resultados.