Comandos R para análises estatísticas

Capítulo 8: Variáveis Aleatórias Multidimensionais

8.1 Introdução

Comecemos com alguns comandos para tabelas de múltipla entrada. Lembremos que alguns comandos já foram apresentados no Capítulo 4. Um pacote interessante é o pacote tables, que constrói tabelas de dupla entrada de forma semelhante ao proc tabulate do SAS através do comando tabular(table,...). Neste comando, o argumento table é uma expressão que especifica as linhas e colunas da tabela desejada. A expressão divide as linhas das colunas através de ~, i.e., antes do sinal ~ especificaremos as linhas e depois dele, as colunas. os sinais + e * indicam se desejamos apenas adicionar as variáveis ou cruzá-las, respectivamente. Além disso, para as variáveis quantitativas de interesse, podemos atribuir também o cálculo de funções, como médias e desvios-padrão, por exemplo. Dessa forma, considerando a Tabela 2.1 do Capítulo 2, se quisermos construir uma tabela com estado civil versus grau de instrução e idade média(em anos), a sintaxe seria:

library(tables)
load("dados.RData")
attach(tab2_1)
tabular( estado_civil ~ grau_instrucao*idade_anos*(mean) )



estado_civil		 grau_instrucao
ensino fundamental
idade_anos
mean
ensino médio
idade_anos
mean
superior
idade_anos
mean
casado 35.80 34.42 37.00
solteiro 37.71 31.50 29.67

Ou seja, temos que, na Tabela 2.1, 5 as pessoas casadas com ensino fundamental tem, em média, 35.8 anos, as casadas com ensino fundamental, 34.42 e assim por diante. Se quisessemos que o grau de instrução ficasse nas linhas, cruzado com o estado civil e, além da média, quiséssemos também o desvio-padrão da idade, a sintaxe seria:

tabular( estado_civil*grau_instrucao ~ idade_anos*(mean+sd) )

estado_civil
grau_instrucao		 idade_anos
mean
sd
casado ensino fundamental 35.80 6.181
ensino médio 34.42 5.282
superior 37.00 9.539
solteiro ensino fundamental 37.71 7.910
ensino médio 31.50 7.918
superior 29.67 4.163

Por fim, para imprimir o total por linhas e colunas, basta adicionar o termo +1 à variável a qual se deseja o total:

tabular( (estado_civil + 1)*(grau_instrucao+1) ~ idade_anos*(mean+sd) )

estado_civil
grau_instrucao		 idade_anos
mean
sd
casado ensino fundamental 35.80 6.181
ensino médio 34.42 5.282
superior 37.00 9.539
All 35.15 5.896
solteiro ensino fundamental 37.71 7.910
ensino médio 31.50 7.918
superior 29.67 4.163
All 33.88 7.805
All ensino fundamental 36.92 6.999
ensino médio 33.44 6.205
superior 33.33 7.711
All 34.58 6.737

Uma observação importante sobre este pacote é o fato de que as variáveis não-quantitativas devem ser da classe factor. Assim, se a variável de interesse não for deste tipo, é necessário primeiramente transformá-la utlizando o comando factor(). Por exemplo, na tabela anterior, caso estado civil não fosse um fator, seria necessário executar o comando da seguinte forma:

tabular( (factor(estado_civil) + 1)*(grau_instrucao+1) ~ idade_anos*(mean+sd) )

Gráficos

Gráfico de dispersão com marginais:

Vamos construir uma função que utilizaremos para imprimir um gráfico de dispersão com as distribuições marginais. O argumento type definirá o tipo de gráfico que teremos nas marginais: “histogram”, “boxplot” ou “dispersion”.

scatter.marginal <- function(x, y, type="histogram",breaks=NA,col="darkred",pch=20,method="jitter",...){
  #  method="jitter" define o tipo de dispersão que teremos na marginal quando type for "dispersion".
  #         Outras opções são "stack", para empilhar os pontos e "overplot", para sobrepô-los.
  
  if(type=="histogram"){
    m=matrix(c(2,0,1,3), ncol=2, byrow=TRUE) # define as regiões a serem utilizadas no gráfico.
    layout(mat=m, widths=c(0.82,0.18), heights=c(0.18,0.82)) 
    if(is.na(breaks)){
      xhist = hist(x, plot=FALSE,breaks=min(max(5,floor(length(x)/7)),10)) # marginal x
      yhist = hist(y, plot=FALSE,breaks=min(max(5,floor(length(x)/7)),10)) # marginal y
    }
    else{
      xhist = hist(x, plot=FALSE,breaks=breaks) # marginal x
      yhist = hist(y, plot=FALSE,breaks=breaks) # marginal y
    }
    top = max(c(xhist$counts, yhist$counts)) # maximo entre os dois eixos
    par(mar=c(5,5,1,1)) 
    plot(x,y,pch=pch,col=col,...) # Grafico principal
    par(mar=c(0,5,1,1))
    barplot(xhist$counts, axes=FALSE, ylim=c(0, top), space=0, col="lightgray",border="white") # imprime o grafico da marginal x
    par(mar=c(5,0,1,1))
    barplot(yhist$counts, axes=FALSE, xlim=c(0, top), space=0, horiz=TRUE,  col="lightgray",border="white") # imprime o grafico da marginal y
  }
  
  if(type=="boxplot"){
    m=matrix(c(2,0,1,3), ncol=2, byrow=TRUE)
    layout(mat=m, widths=c(0.8,0.2), heights=c(0.18,0.82))
    par(mar=c(5,5,1,1))
    plot(x,y,pch=pch,col=col,...) # Grafico principal
    par(mar=c(0,5,1,1))
    boxplot(x, axes=FALSE, col="lightgray",pch=20, horizontal=TRUE) # marginal x
    par(mar=c(5,0,1,1))
    boxplot(y, axes=FALSE, col="lightgray",pch=20) # marginal y
  }
  
  if(type=="dispersion"){
    m=matrix(c(2,0,1,3), ncol=2, byrow=TRUE)
    layout(mat=m, widths=c(0.85,0.15), heights=c(0.18,0.82))
    par(mar=c(5,5,1,1))
    plot(x,y,pch=pch,col=col,...) # gráfico principal
    par(mar=c(0,5,1,1))
    stripchart(x, method = method, pch = 20, col="gray",vertical = FALSE,yaxt="n",xaxt="n") # marginal x
    par(mar=c(5,0,1,1))
    stripchart(y, method = method, pch = 20, col="gray",vertical = TRUE,yaxt="n",xaxt="n") # marginal y
  }
}

Para imprimir o gráfico utilizaremos a função que acabamos de construir:

load("dados.RData")
scatter.marginal(tab2_1$idade_anos, tab2_1$salario, ylab="salario", xlab="idade")

Figura 8.A: Gráfico de dispersão utilizando histogramas nas marginais

scatter.marginal(tab2_1$idade_anos, tab2_1$salario, ylab="salario", xlab="idade",type="boxplot")

Figura 8.B: Gráfico de dispersão utilizando boxplots nas marginais

scatter.marginal(tab2_1$idade_anos, tab2_1$salario, ylab="salario", xlab="idade",type="dispersion")

Figura 8.C: Gráfico de dispersão utilizando dispersão nas marginais

Pode-se também utilizar os pacotes ggplot2 e ggExtra para a construção destes gráficos. Os comandos para a construção das figuras a 8.A a 8.C são:

attach(tab2_1)
library(ggplot2)
library(ggExtra)

dispersao = ggplot(tab2_1, aes(idade_anos,salario)) + geom_point()  #Dispersao
ggMarginal(dispersao, type="histogram") # Marginais com histograma

dispersao = ggplot2::ggplot(tab2_1, ggplot2::aes(idade_anos, salario)) + ggplot2::geom_point() #Dispersao
ggMarginal(dispersao,type = "boxplot") # Marginais com boxplot

dispersao = ggplot2::ggplot(tab2_1, ggplot2::aes(idade_anos, salario)) + ggplot2::geom_point() #Dispersao

scatter <- ggplot(tab2_1, aes(idade_anos,salario))   +  geom_point() + #dispersão
         scale_x_continuous(limits=c(min(idade_anos),max(idade_anos))) + # marginal X
         scale_y_continuous(limits=c(min(salario),max(salario))) +  # marginal Y
         geom_rug(size=0.1) + theme_set(theme_minimal(base_size = 18)) 
scatter

É possivel também imprimir as marginais em formato de densidade:

ggMarginal(dispersao,type = "density") # Marginais com função densidade empírica

Gráficos de perfis - Gráficos de dispersão Cruzados

Outro gráfico útil para análise de v.a.’s multivariadas é o gráfico de dispersão em pares, ou gráficos de perfis. São úteis especialmente quando temos 3 ou mais variáveis e temos interesse em analisar o comportamento de uma em função das outras. Podemos construí-los com a função pairs(...). Ainda considerando a Tabela 2.1, Podemos cruzar as variáveis n_filhos, salario e idade_anos para identificar possíveis padrões de correlação entre elas:

pairs( ~ n_filhos + salario + idade_anos, pch=16, col="blue")

Note que os gráficos abaixo da diagonal principal são equivalentes os gráficos acima da diagonal principal, transpostos. Um opção útil para situações em que há muitas variáveis é imprimir apenas abaixo (ou acima) da diagonal principal, utilizando a opção upper.panel = NULL:

pairs( ~ n_filhos + salario + idade_anos, pch=16, col="blue",upper.panel = NULL)

OBS: para suprimir o os gráficos abaixo da diagonal, a opção é lower.panel=NULL.

Gráficos de curva de nível e mapas de calor

Outros gráficos úteis são os gráficos de curvas de níveis, que são utilizados quando temos 3 variáveis contínuas, sendo duas explicativas(digamos, \(X\) e \(Y\)) e 1 resposta(\(Z\)). Neles, as vairações na escala de \(Z\) são apresentadas em função de cortes projetados no plano \(X \times Y\). São semelhantes aos “mapas de calor”, que apresentam as diferenças entre os níveis por variações de tons de cor. O pacote graphics possui a função filled.contour, que produz tais mapas de calor. Mapas de calor também são conhecidos como gráficos de contorno.

No exemplo a seguir, com dados extraídos de help(filled.contour), apresentamos a altitude do vulcão Maunga Whau, na Nova Zelândia, em função de sua latitude e longitude. No gráfico da Figura 8.A temos as curvas de nível e na Figura 8.B, o mapa de calor.

contour(x=1:nrow(volcano), y=1:ncol(volcano), volcano, col="darkred",xlab = "Sul - Norte", ylab = "Leste - Oeste")
title("Mapa topográfico: Maunga Whau", font = 4) # Titulo do Gráfico

Figura 8.A: Gráfico de curvas de nível do vulcão Maunga Whau em função da latitude e longitude.

library(colorRamps) # Paleta de cores que contem a sequencia de cores `blue2red`.
filled.contour(x=1:nrow(volcano), y=1:ncol(volcano), z=volcano, color = blue2red,
    plot.title = title(main = "Mapa topográfico: Maunga Whau \n ",
    xlab = "Sul - Norte", ylab = "Leste - Oeste"),
    plot.axes = { axis(1, seq(100, 800, by = 100))
                  axis(2, seq(100, 600, by = 100)) },
    key.axes = axis(4, seq(90, 190, by = 10)))# outra sugestão de cor: `color = terrain.colors` ou não declarar esta opção para visualizar a cor padrão.

Figura 8.B: Mapa de calor do vulcão Maunga Whau em função da latitude e longitude.

8.4 Covariância amostral

Para calcular a covariância amostral entre duas variáveis X e Y, utilize o comando cov(X,Y). Considerando os dados de salário e idade em anos temos que a covariância entre elas é:

cov(salario,idade_anos)
## [1] 11.231
cov(idade_anos,salario)
## [1] 11.231

Note que a ordem de entrada das variáveis não influencia no valor da covariância.

8.8 Distribuição Normal Bidimensional

Para simular uma v.a. Normal multivariada, pode-se utilizar a função rmvnorm da bilioteca mvtnorm. Simulando então Um v.a. Normal bivariada X, com média \(\mu\) e variância \(\Sigma\) dadas por:

\[ \mu=\left(\begin{array}{c} 0 \\ 0 \end{array}\right)\\ \Sigma=\left(\begin{array}{cc} 1 & 0.5 \\ 0.5 & 1 \end{array}\right) \]

Teremos:

library(mvtnorm)
Sigma=cbind(c(1,0.6),c(0.6,1)) # matriz de variancias e covariancias
mu=c(0,0) # Vetor de médias
Z.sim<-rmvnorm(n = 1000, mean=mu, sigma=Sigma) # simula 100 observações de uma normal bivariada

Para calcular probabilidade conjunta, isto é, a F.D.A., utilizamos a função dmvnorm para construir uma função mais amigável para o usuário:

fdamvnorm<-function(upper,lower=c(-Inf,-Inf),mu=rep(0,times=length(upper)),Sigma=diag(length(upper))){
return(pmvnorm(lower=c(-Inf,-Inf),upper=upper,mean=mu,sigma=Sigma)[[1]]) # Calcula a integral dupla da fdp da normal multivariada entre no retangulo `lower X upper`.
}
fdamvnorm(upper=c(0,0))   # FDA considerando o terceiro quadrante
## [1] 0.25
fdamvnorm(upper=c(Inf,0)) # FDA considerando o terceiro e quarto quadrantes
## [1] 0.5
fdamvnorm(c(0,Inf))       # FDA considerando o segundo e terceiro quadrantes
## [1] 0.5
fdamvnorm(lower=c(-1.96,-1.96),upper=c(1.96,1.96))  #FDA considerando o quadrado  (-1.96,-1.96)X(1.96,1.96) 
## [1] 0.95063

A função densidade conjunta da distribuição normal multivariada é calculada com a função dmvnom, com a sintaxe muito semelhante à sintaxe da função rmvnorm. Se quisermos, por exemplo, o valor da curva da função no ponto (0,0), para \(\mu\) e \(\Sigma\) como definidos anteriormente, teremos que o valor é 0.19894:

Sigma=cbind(c(1,0.6),c(0.6,1)) # matriz de variancias e covariancias
mu=c(0,0) # Vetor de médias
dmvnorm(x=c(0,0), mean=mu, sigma=Sigma)
## [1] 0.19894

Figura 8.15

Precisamos calcular o valor da função densidade para uma região do plano \(XY\). Vamos então percorrer estes dois eixos entre os valores -4 e 4, amostrando a função densidade a cada 0.1:

xy<-seq(-4,4,0.1)
mu=c(0,0)
Sigma=diag(2)
Z1<-matrix(NA,length(xy),length(xy))
for(i in 1:length(xy)){
  for(j in 1:length(xy)){
    Z1[i,j]<-dmvnorm(x=c(xy[i],xy[j]), mean=mu, sigma=Sigma)
  }
}

Para construir o gráfico de dispersão 3D de uma v.a. normal multivariada, utilizaremos a função persp do pacote graphics. Basta então utilizar o comando abaixo:

persp(Z1, phi = 30, theta = -45,col=heat.colors(ncol(Z1)*nrow(Z1),.8), ticktype = "detailed",xlab="X")

# `phi` e `theta` são os parâmetros de deslocamento da perspectiva da figura.
# heat.colors(.) modifica a paleta de cores do gráfico
# box=FALSE elimina a caixa que envolve o gráfico
# scale =FALSE define que as variáveis não devem ser padronizadas  
# ticktype: tipo de eixos. Pode assumir valores "simple" ou "detailed"

Figura 8.15(a): f.d.p. de normais bidimensionais \(\mu_x=\mu_y=0,\sigma_x\sigma_y=1,\rho=0\);

Modificando agora o valor de \(\rho\) para \(\rho=0.6\) temos:

xy<-seq(-4,4,0.1)
mu=c(0,0)
Sigma=cbind(c(1,0.6),c(0.6,1))
Z2<-matrix(NA,length(xy),length(xy))
for(i in 1:length(xy)){
  for(j in 1:length(xy)){
    Z2[i,j]<-dmvnorm(x=c(xy[i],xy[j]), mean=mu, sigma=Sigma)
  }
}
persp(Z2, phi = 30, theta = -45,col=heat.colors(ncol(Z2)*nrow(Z2),.8), ticktype = "detailed", xlab="X")

Figura 8.15(b): f.d.p. de normais bidimensionais \(\mu_x=\mu_y=0,\sigma_x\sigma_y=1,\rho=0.6\);

Figura 8.17

Como vimos anteriormente, as curvas de nível são construídas a partir do comando contour:

par(mfrow=c(1,2))
contour(x=1:nrow(Z2), y=1:ncol(Z2), Z2, col="darkred",xlab = "X", ylab = "Y")
contour(x=1:nrow(Z1), y=1:ncol(Z1), Z1, col="darkred",xlab = "X", ylab = "Y")

Figura 8.17: Curvas de nível para a normal bidimensional

Capítulo Anterior | Introdução | Próximo Capítulo