Estatística Básica - Bussab e Morettin (2017)
Comandos R para análises estatísticas
Capítulo 3: Medidas-Resumo
3.1 Medidas de Posição
Exemplo 3.1
Para calcular a média, utilize a função mean:
mean(x, na.rm = FALSE, ...)
# na.rm: Se este argumento for TRUE, então os valores NA são desconsiderados para o calculo da média.
mean(tab2_1$n_filhos,na.rm=T)## [1] 1.65Exemplo 3.2
Não há, no R-base uma função específica para o cáculo da moda. No entanto, é possível calculá-la facilmente a partir da função table ou da função hist que utilizamos no capítulo anterior. Relembre que para construir a Figura 2.7 utilizamos o seguinte código:
fig27<-hist(tab2_1$salario, breaks = seq(4,24,by=4),right=FALSE, plot=F)Para calcular a moda baseada nos cortes do histograma, construíremos uma função moda:
moda<-function(x,breaks,...){
h<-hist(x,breaks=breaks,plot=F,...)
h$mids[h$counts==max(h$counts)]
}Para utilizá-la, faremos:
moda(tab2_1$salario, breaks = seq(4,24,by=4),right=FALSE)## [1] 10Uma forma alternativa de construir esta função é considerar o seguinte código:
# Utilize o comando a seguir para construir a função `moda`
moda2<-function(x){return(names(sort(-table(as.vector(x))))[1])}No entanto, nesta função, o cálculo da moda não será realizado de acordo com os cortes/pontos médios do histograma. É uma função mais indicada se a variável de interesse for nominal ou discreta.
Para calcular a mediana, utilize a função median:
median(tab2_1$salario)## [1] 10.165Exemplo 3.2 (continuação)
Utilizando a função moda2 acima, calculemos:
moda2(tab2_1$grau_instrucao)## [1] "ensino médio"Já para o cálculo da mediana, precisamos transformar a variável grau_instrucao em uma variável numerica:
tab2_1$grau_instr_num<-1*c(tab2_1$grau_instrucao=="ensino fundamental")+2*c(tab2_1$grau_instrucao=="ensino médio")+3*c(tab2_1$grau_instrucao=="superior")Agora, calculando a mediana, temos:
median(tab2_1$grau_instr_num)## [1] 2Ou, se prefirmos que apareça o nome da categoria:
levels(as.factor(tab2_1$grau_instrucao))[median(tab2_1$grau_instr_num)]## [1] "ensino médio"Para calcular a “média aparada” ou “truncada” desprezando-se, por exemplo, as 1% maiores e 1% menores observações, podemos fazer:
#quebras de linha apenas lustrativas para facilitar a leitura
mean(tab2_1$salario[
tab2_1$salario>quantile(tab2_1$salario,probs=0.01)
&
tab2_1$salario<quantile(tab2_1$salario,probs=0.99)
])## [1] 10.974#Note que este valor é menor que o valor da média da variável sem a "truncagem":
mean(tab2_1$salario)## [1] 11.1223.2 Medidas de Dispersão
Exemplo 3.3
O R calcula a variância e desvio padrão sempre considerando o estimador não viesado, isto é, considerando N-1 como denominador. Para calcular a variância e desvio padrão considerando o denominador com N, faremos como se segue:
n<-sum(!is.na(tab2_1$n_filhos))
#Variância:
var(tab2_1$n_filhos,na.rm=T)*(n-1)/n## [1] 1.5275#Desvio padrão:
sd(tab2_1$n_filhos,na.rm=T)*(n-1)/n## [1] 1.2046#Desvio médio:
mean(abs(tab2_1$n_filhos[!is.na(tab2_1$n_filhos)]-mean(tab2_1$n_filhos,na.rm=T)))## [1] 0.985A função is.na(x) verifica, para cada objeto do vetor x, quais deles são missing e devolve TRUE ou FALSE como resposta ao passo que a função ! gera o complementar do objeto lógico sucedente.
# Exemplo:
x<-c(0,1,NA,NA,0,-2,pi,NA)
is.na(x)## [1] FALSE FALSE TRUE TRUE FALSE FALSE FALSE TRUE!is.na(x)## [1] TRUE TRUE FALSE FALSE TRUE TRUE TRUE FALSEExemplo 3.4
Eventualmente, para evitar os passos adicionais que foram inseridos no item anterior para o cálculo da variância amostral, podemos criar uma nova função do R, considerando a função var como base, mas com os passos adicionais. Logo:
varp<-function(x, pop=TRUE, na.rm=TRUE,...){
if (pop==TRUE & na.rm==TRUE){ # calcula var_n desconsiderando os missings
n<-sum(!is.na(x))
return(var(x,na.rm=na.rm,...)*(n-1)/n)
}
else{
if(na.rm==FALSE){
n<-length(x)
return(var(x,na.rm=na.rm,...)*(n-1)/n)
}
}
return(var(x,na.rm=na.rm,...))
}Da mesma forma, podemos construir uma função para o cálculo do desvio médio:
dm<-function(x,...){
return(mean(abs(x[!is.na(x)]-mean(x,...))))
}A partir daí, pode-se usar as funções dm e varp da mesma maneira que uma função já previamente declarada, como a seguir:
mean(tab2_1$salario,na.rm=T)## [1] 11.122dm(tab2_1$salario)## [1] 3.7309#n<-sum(!is.na(tab2_1$n_filhos))
varp(tab2_1$salario,na.rm=T)## [1] 20.463.3 Quantis Empíricos
Exemplo 3.5
Vamos declarar o vetor x e calcular seus quantis:
x<-c(15,5,3,8,10,2,7,11,12)
median(x)## [1] 8mean(x)## [1] 8.1111quantile(x,probs=c(.25,.50,.75),type=5) # "type 6" indica o mesmo método de cálculo que o SPSS e Minitab## 25% 50% 75%
## 4.50 8.00 11.25#Agora, adicionando a observação 67 ao vetor x, temos:
x<-c(sort(x),67) # A função `sort` ordena o vetor
x## [1] 2 3 5 7 8 10 11 12 15 67median(x)## [1] 9mean(x)## [1] 14quantile(x,probs=c(.25,.50,.75),type=5) ## 25% 50% 75%
## 5 9 12# 20% das observações a esquerda
quantile(x,probs=0.2)## 20%
## 4.6Exemplo 3.6
Lembre-se que, para a construção do histograma, utilizamos o seguinte código:
fig27<-hist(tab2_1$salario, breaks = seq(4,24,by=4),right=FALSE,probability = T,plot=F)
aux<-with(fig27, 100 * density* diff(breaks)[1])
labs <- paste(round(aux), "%", sep="")
#quebras de linha apenas ilustrativas para facilitar a leitura
plot(fig27,
freq = FALSE, labels = labs,
xlab="Salário",
ylab="",
col="darkgrey",
border="white",
yaxt="n",
xlim=c(0,24), xaxp=c(0,24,6),
ylim=c(0,.1), main="")
median(tab2_1$salario)## [1] 10.165quantile(tab2_1$salario,probs=c(.20,.5,.95,.75))## 20% 50% 95% 75%
## 7.390 10.165 18.913 14.060# Distancia inter-quartis:
q<-quantile(tab2_1$salario,probs=c(.25,.5,.75)) # armazenando as informações de quantis em um objeto `q`
q## 25% 50% 75%
## 7.5525 10.1650 14.0600dq=q[3]-q[1]
print(paste("Distância inter-quartis (d) = ", dq))## [1] "Distância inter-quartis (d) = 6.5075"Figuras 3.2 e 3.3
Para se ter o que se chama de Esquema dos cinco números, podemos utilizar a função summary do R:
# Utilizando os dados do exemplo 3.5 :
x<-c(15,5,3,8,10,2,7,11,12)
summary(x)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.00 5.00 8.00 8.11 11.00 15.00Exemplo 3.7
Note que a forma para calcular as medidas resumo no R é a mesma que no S-Plus:
cd_municipios<-read.table("cd-municipios.csv",h=T,skip=4,sep=";",dec=",")
# Medidas resumo para a base de dados dos municípios (n=30)
summary(cd_municipios$populacao,digits = 4)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 46.30 64.48 84.30 145.40 134.30 988.80# Medidas resumo para a base de dados dos municípios, desconsiderando São Paulo e Rio de Janeiro (n=28)
summary(cd_municipios$populacao[-c(1,2)],digits = 4) ## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 46.30 63.48 82.05 100.60 122.00 224.603.4 Box plots
Exemplo 3.8
Primeiro, precisamos restringir a análise aos 15 maiores municípios:
munic<-cd_municipios[order(cd_municipios$populacao,decreasing = TRUE),]
munic## n municipio populacao
## 1 1 São Paulo(SP) 988.8
## 2 2 Rio de Janeiro(RJ) 556.9
## 3 3 Salvador(BA) 224.6
## 4 4 Belo Horizonte(MG) 210.9
## 5 5 Fortaleza(CE) 201.5
## 6 6 Brasília(DF) 187.7
## 7 7 Curitiba(PR) 151.6
## 8 8 Recife(PE) 135.8
## 9 9 Porto Alegre(RS) 129.8
## 10 10 Manaus(AM) 119.4
## 11 11 Belém(PA) 116.0
## 12 12 Goiânia(GO) 102.3
## 13 13 Guarulhos(SP) 101.8
## 14 14 Campinas(SP) 92.4
## 15 15 São Gonçalo(RJ) 84.7
## 16 16 Nova Iguaçu(RJ) 83.9
## 17 17 São Luis(MA) 80.2
## 18 18 Maceió(AL) 74.7
## 19 19 Duque de Caxias(RJ) 72.7
## 20 20 São Bernardo do Campo(SP) 68.4
## 21 21 Natal(RN) 66.8
## 22 22 Teresina(PI) 66.8
## 23 23 Osasco(SP) 63.7
## 24 24 Santo André(SP) 62.8
## 25 25 Campo Grande(MS) 61.9
## 26 26 João Pessoa(PB) 56.2
## 27 27 Jaboatão(PE) 54.1
## 28 28 Contagem(MG) 50.3
## 29 29 São José dos Campos(SP) 49.7
## 30 30 Ribeirão Preto(SP) 46.3munic<-cd_municipios[1:15,]
munic## n municipio populacao
## 1 1 São Paulo(SP) 988.8
## 2 2 Rio de Janeiro(RJ) 556.9
## 3 3 Salvador(BA) 224.6
## 4 4 Belo Horizonte(MG) 210.9
## 5 5 Fortaleza(CE) 201.5
## 6 6 Brasília(DF) 187.7
## 7 7 Curitiba(PR) 151.6
## 8 8 Recife(PE) 135.8
## 9 9 Porto Alegre(RS) 129.8
## 10 10 Manaus(AM) 119.4
## 11 11 Belém(PA) 116.0
## 12 12 Goiânia(GO) 102.3
## 13 13 Guarulhos(SP) 101.8
## 14 14 Campinas(SP) 92.4
## 15 15 São Gonçalo(RJ) 84.7# Medidas resumo para a base de dados dos 15 municípios mais populosos (n=15)
summary(munic$populacao)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 84.7 109.0 136.0 227.0 206.0 989.0Figuras 3.6 e 3.7
Para construir o box plot, utilize a função boxplot.
boxplot(munic$populacao,
pch="*", # tipo de marcador dos outliers
col="lightblue", # cor do preenchimento do box plot
border="darkgrey", # cor da linha do box plot
boxwex=0.3 # Tamanho da caixa
)
# Os comandos 'text'a seguir imprimem no box plot os nomes das cidades dos 4 pontos destacados
text(x=cd_municipios$populacao[1],label=cd_municipios$municipio[1],pos=4,cex=0.7) # Máximo
text(x=cd_municipios$populacao[2],label=cd_municipios$municipio[2],pos=4,cex=0.7) # 2a. maior observacao
text(x=cd_municipios$populacao[3],label=cd_municipios$municipio[3],pos=4,cex=0.7) # 3a. maior observacao
text(x=cd_municipios$populacao[15],label=cd_municipios$municipio[15],pos=4,cex=0.7) # 15a. observacao 
Figura 3.6: Box plot para os quinze maiores municípios do Brasil“.
#Para imprimir o box plot desconsiderando aqueles pontos que foram considerados outliers de acordo com seu critério, faça:
boxplot(munic$populacao, col="lightblue", border="darkgrey", boxwex=0.3, outline=FALSE)
title("Box plot para os quinze maiores \n municípios do Brasil, desconsiderando outliers")
O tipo do marcador, informação no argumento pch, dentre outros como lty(tipo da linha)lwd, cex etc. é um argumento padrão para todos os gráficos no R. Veja ?par para maiores detalhes.
Tipos de marcadores para pch: 
3.5 Gráficos de Simetria
Exemplo 3.9
x<-c(0.5,2.3,4,6.4,8,15.3,13.5,12,9.8)
# Mediana de x
median(x)## [1] 8# Calculando os ui e vi
u<-sort((median(x)-x)[(median(x)-x)>0],decreasing = T)
v<-sort((x-median(x))[(x-median(x))>0],decreasing = T)
u## [1] 7.5 5.7 4.0 1.6v## [1] 7.3 5.5 4.0 1.8Figura 3.11
u<-median(cd_municipios$populacao)-cd_municipios$populacao
v<-cd_municipios$populacao-median(cd_municipios$populacao)
plot(sort(u),sort(v), pch=19, xlab="ui", ylab="vi",col="darkblue",xlim=c(0,max(u)),ylim=c(0,max(v)))
#title("Figura 3.11: Gráfico de simetria para o CD-Municípios.")
abline(0,1)
Figura 3.11: Gráfico de simetria para o CD-Municípios.
Transformações
Exemplo 3.10 e Figura 3.12
Para as Figura 3.12 e 3.13, criaremos um objeto chamado list, que é basicamente uma classe na qual o objeto contém outros objetos em seu conteúdo, podendo esses objetos serem de classes diferentes.
xp<-list() # declarando xp como uma lista vazia
xp[[1]]<-log(cd_municipios$populacao) # Transformação log
xp[[2]]<-(cd_municipios$populacao)^(1/4) # Transformação p=1/4
xp[[3]]<-(cd_municipios$populacao)^(1/2) # Transformação p=1/2
xp[[4]]<-(cd_municipios$populacao)^(1/3) # Transformação p=1/3par(mfrow=c(2,2))
hist(xp[[1]], main="log", ylab="", xlab="", col="darkgrey", border="white",cex.axis=0.8,cex.main=0.8)
hist(xp[[2]], main="p=1/4", ylab="", xlab="", col="darkgrey", border="white",cex.axis=0.8,cex.main=0.8)
hist(xp[[3]], main="p=1/2", ylab="", xlab="", col="darkgrey", border="white",cex.axis=0.8,cex.main=0.8)
hist(xp[[4]], main="p=1/3", ylab="", xlab="", col="darkgrey", border="white",cex.axis=0.8,cex.main=0.8)
Figura 3.12: Histogramas para os dados transformados. CD-Municípios.
Figura 3.13
boxplot(xp,pch="-", col="lightblue", border="darkgrey", boxwex=0.3, names=c("p=0","p=1/4","p=1/2","p=1/3"))
Figura 3.13: Box plots para os dados transformados. CD-Municípios. R.
3.7 Exemplos Computacionais
Exemplo 2.10 (continuação)
Como a tabela gerada pelo summary do R contém menos informações que a tabela descritiva do Minitab, adicionemos as informações faltantes de modo a termos tabelas equivalentes.
#Podemos construir uma nova função, digamos summary2, que conterá mais informação do que a summary atual:
summary2<-function(x,limtrim=c(.01,.99), pop=FALSE, na.rm=TRUE, ...){
# x : argumento numerico para o qual calcularemos as estatísticas descritivas
# limtrim: Limites para a média truncada
# pop: se esta variável é TRUE, entao utiliza-se 'n' como denominador da estimativa da variância. caso contrário, 'n-1'.
# na.rm: ação para os outliers
newsumm<-matrix(NA,nrow=dim(as.matrix(x))[2],ncol=11)
x<-as.matrix(x)
newsumm[,1]<-apply(x,2,function(...){return(dim(x)[1])})
newsumm[,2:7]<-t(apply(as.data.frame(x),2,summary))
newsumm[,8]<-apply(x,2,function(...){mean(x[x>quantile(x,probs=0.01) &
x<quantile(x,probs=0.99)])})
if (pop==TRUE & na.rm==TRUE){ # calcula var_n desconsiderando os missings
n<-sum(!is.na(x))
newsumm[,9]<-var(x,na.rm=na.rm)*(n-1)/n
}
else
{
if(pop==TRUE & na.rm==FALSE ){
n<-length(x)
newsumm[,9]<-var(x)*(n-1)/n
}
else {
if(pop==FALSE & na.rm==TRUE){
n<-sum(!is.na(x))
newsumm[,9]<-var(x,na.rm=TRUE)
}
else
{
n<-length(x)
newsumm[,9]<-var(x)
}
}
}
newsumm[,1]<-n
newsumm[,10]<-sqrt(newsumm[,9])
newsumm[,11]<-sqrt(newsumm[,10]/newsumm[,1])
colnames(newsumm)<-c("N", "Min.","1st Qu.","Median","Mean","3rd Qu.","Max.","Tr Mean","Var","StDev","SE Mean")
return(t(newsumm))
}
cd_notas<-read.table("cd-notas.csv",header=TRUE,skip=4,sep=";",dec=",") # Leitura dos dados
summary2(cd_notas$nota)## [,1]
## N 100.00000
## Min. 1.50000
## 1st Qu. 4.88000
## Median 6.00000
## Mean 5.92000
## 3rd Qu. 7.12000
## Max. 10.00000
## Tr Mean 5.84375
## Var 3.28472
## StDev 1.81238
## SE Mean 0.13462Figura 3.14
boxplot(cd_notas$nota,pch="-", col="lightblue", border="darkgrey")
Figura 3.14: Box plot para o CD-Notas. R.
Figura 3.15
u<-median(cd_notas$nota)-cd_notas$nota
v<-cd_notas$nota-median(cd_notas$nota)
plot(sort(u),sort(v), pch=19, xlab="ui", ylab="vi",col="darkblue",xlim=c(0,max(u)),ylim=c(0,max(v)))
title("Figura 3.15: Gráfico de simetria para o CD-Notas.")
abline(0,1)
Exemplo 2.11 (continuação)
cd_poluicao<-read.table("cd-poluicao.csv",header=TRUE,skip=8,sep=";",dec=",")
# Quadro 3.4 Medidas descritivas para temperaturas. R.
summary(cd_poluicao$temp)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 12.3 16.0 17.7 17.2 18.6 21.0Figura 3.16
boxplot(cd_poluicao$temp, pch="-", col="lightblue", border="darkgrey")
Figura 3.16: Box plot para as temperaturas de São Paulo. CD-Poluição. R.
Figura 3.17
u<-median(cd_poluicao$temp)-cd_poluicao$temp
v<-cd_poluicao$temp-median(cd_poluicao$temp)
plot(sort(u),sort(v), pch=19, xlab="ui", ylab="vi",col="darkblue",xlim=c(0,max(u)),ylim=c(0,max(v)))
abline(0,1)
Figura 3.17: Gráfico de simetria para as temperaturas de São Paulo. CD-Poluição.
Figura 3.18
x<-c(15,5,3,8,10,2,7,11,12)
plot(ecdf(x),xlim=c(0,16),ylim=c(0,1),main="",ylab="",pch=20,cex=0.7) # 'ecdf' é a função que calcula a função distribuição empírica
lines(y = c(0.5,0.5),x=c(-1,8),col="grey")
lines(y = c(0,0.5),x=c(8,8),col="grey")
text(-.5,0.5,"p=9/18",pos=4,col="darkblue")
text(8,0,"8",pos=3,col="darkblue")
#title("Figura 3.18: Funções de distribuição empírica(F_e) \n e f.d.e. alisada (Fs_e) para o Exemplo 3.5")
par(new=T) # Utilize este comando para sobrepor o gráfico a seguir no gráfico anterior
# Construiremos a seguir uma função 'scdf', que será a ecdf acima suavisada
scdf<-function(x){
f<-t<-table(x)
st<-sum(t[1:length(x)])
for(i in 1:length(x))
f[i]<-sum(t[1:i])/st - (t[i]/st)/2
return(f)
}
par(new=T)
plot(scdf(x),pch=20,cex=0.7,lty=2, type="b",xlab="",ylab="",col="grey",xaxt="n",yaxt="n",xlim=c(0,16),ylim=c(0,1))
legend(12,.4,c("F_e","Fs_e"),lty=c(1,2),col=c("black","grey"),lwd=2,cex=0.8)
Figura 3.18: Funções de distribuição empírica(F_e) e f.d.e. alisada (Fs_e) para o Exemplo 3.5.