Teste para a diferença de médias entre duas populações

Conjunto de Dados

O conjunto de dados usado para este exemplo contém informação sobre indivíduos do sexo masculino e feminino, tais como o nível de colesterol no sangue. Iremos focar-nos nessa variável e na variável sexo apenas.

Colestrol=data.frame(Base$CHOL,Base$Sexo)
colnames(Colestrol)=c("Colesterol","Sexo")
head(Colestrol)
##   Colesterol Sexo
## 1       3.23    1
## 2       4.80    1
## 3       5.20    1
## 4       4.74    1
## 5       4.32    1
## 6       6.05    1

Teste de Hipóteses

Inicialmente vamos comparar, por observação, a distribuição amostral do nível de colesterol no sangue, correspondente a cada um dos sexos. O gráfico seguinte fornce essa informação (a azul os homens e a cor-de-rosa as mulheres).

ggplot(HCV,aes(x=CHOL,fill=Sex))+geom_density(alpha=0.3)+theme_minimal()+
  labs(x="Colesterol",caption=c("Mulheres","Homens"))

À partida, as duas populações parecem ter características bastante semelhantes contudo, vamos realizar um teste de hipótese para avaliar a veracidade dessa suposição.

t.test(x=Mulheres,y=Homens,var.equal=T,mu=0)
## 
##  Two Sample t-test
## 
## data:  Mulheres and Homens
## t = 0.73801, df = 587, p-value = 0.4608
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.1173238  0.2585734
## sample estimates:
## mean of x mean of y 
##  5.434867  5.364242

Obtivemos um valor-p superior a alpha, pelo que não rejeitamos a hipótese nula, e por isso concluimos que não existem diferenças significativas entre as duas populações (homens e mulheres), relativamente ao nível de colesterol no sangue.

MediasAmos=tapply(Colestrol$Colesterol, Colestrol$Sexo, mean)

Diferenca=diff(tapply(Colestrol$Colesterol, Colestrol$Sexo, mean))

MedAmos=-0.07062483 # diferença das médias amostrais

O nosso objetivo é agora aplicar o bootstrap neste contexto e comparar os resultados.

set.seed(4)
N=10000
amostra=Base$CHOL
am=replicate(N, diff(tapply(sample(amostra, replace = TRUE), Colestrol$Sexo, mean)))
hist(am, main = "Distribuição amostral",col = "lightskyblue",xlab = "Diferença de Médias")
abline(v=MedAmos,col=2)

O valor-p calcula-se tomando a proporção de valores da estatística, calculados para cada reamostra, que foi superior à diferença de médias original.

sum(am>=MedAmos)/N
## [1] 0.7712

Intervalos de Confiança

Considerando apenas a população dos homens, vamos utilizar o Bootstrap para obter intervalos de confiança para a média do nível de colesterol. A função boot do pacote boot permite-nos obter os intervalos de confiança normal, percentil e básico, abordados no relatório.

library(boot)

meanfun <- function(Homens,i){
  d <- Homens[i]
  return(mean(d))   
}

set.seed(46)
bo=boot(data=Homens,statistic = meanfun ,R=2000)
boot.ci(bo,type=c("norm","basic","perc"))
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 2000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = bo, type = c("norm", "basic", "perc"))
## 
## Intervals : 
## Level      Normal              Basic              Percentile     
## 95%   ( 5.241,  5.483 )   ( 5.240,  5.486 )   ( 5.243,  5.489 )  
## Calculations and Intervals on Original Scale

Regressão linear via Bootstrap

O bootstrap de um modelo de regressão linear fornece uma visão de quão variantes são os parâmetros do modelo, ou seja, concede uma visão da variação aleatória existente nos coeficientes de regressão. Como a reamostragem de bootstrap usa um grande número de subamostras, o modelo criado pode ser complexo e computacionalmente intenso.

Conjunto de dados

O conjunto de dados usado possuí duas variáveis, YearsExperience and Salary. A regressão linear usa como variável explicativa YearsEperience e como variável objetivo Salary.

summary(Salary_Data)
##  YearsExperience      Salary      
##  Min.   : 1.100   Min.   : 37731  
##  1st Qu.: 3.200   1st Qu.: 56721  
##  Median : 4.700   Median : 65237  
##  Mean   : 5.313   Mean   : 76003  
##  3rd Qu.: 7.700   3rd Qu.:100545  
##  Max.   :10.500   Max.   :122391

Vamos inicialmente ver como estão distribuídos os dados originais.

hist(Salary_Data$Salary)

hist(Salary_Data$YearsExperience)

rsq_function

rsq_function <- function(formula, data, indices) {
d <- data[indices,] #allows boot to select sample
fit <- lm(formula, data=d) #fit regression model
return(summary(fit)$r.square) #return R-squared of model
}

Para calcular o erro quadrado de cada modelo foi criada uma função para ser chamada pelo Bootstrap. Nesta função usamos a função standard lm para a modelagem, e usamos o seu parametro r.squared para registar o erro quadrado a cada iteração.

Bootstrap Não Paramétrico

rsq_function <- function(formula, data, indices) {
d <- data[indices,] #allows boot to select sample
fit <- lm(formula, data=d) #fit regression model
return(summary(fit)$r.square) #return R-squared of model
}

reps <- boot(data=Salary_Data, statistic=rsq_function, R=2000, formula=Salary~YearsExperience)

reps
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = Salary_Data, statistic = rsq_function, R = 2000, 
##     formula = Salary ~ YearsExperience)
## 
## 
## Bootstrap Statistics :
##      original        bias    std. error
## t1* 0.9569567 -0.0004034265  0.01276198
plot(reps)

Podemos imediatamente destes resultados observar uma melhoria muito significativa no erro. Sendo que o R-Squared(erro quadrado) estimado para este modelo de regressão é 0.9569567, enquanto que com o erro padrão para esta simulação, com auxilio de bootstrap não paramétrico, é 0.01247134. Podemos deste valor de erro padrão concluir que as amostras são representativas do conjunto de dados.

Nos gráficos gerados pela função bootstrap, t representa a estatística que estamos a analisar, neste caso o erro quadrado. Podemos também observar em ambos os gráficos uma linha tracejada, esta linha indica a estatística dos dados originais. Este gráficos reforçam a vantagem do bootstrap não paramétrico em aplicações de regressão linear, em particular em aplicações com uma amostragem original reduzida.

Bootstrap Paramétrico

rsq_function <- function(formula, data, indices) {
d <- data[indices,] #allows boot to select sample
fit <- lm(formula, data=d) #fit regression model
return(summary(fit)$r.square) #return R-squared of model
}

reps <- boot(data=Salary_Data,sim = "parametric", statistic=rsq_function, R=2000, formula=Salary~YearsExperience)

reps
## 
## PARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = Salary_Data, statistic = rsq_function, R = 2000, 
##     sim = "parametric", formula = Salary ~ YearsExperience)
## 
## 
## Bootstrap Statistics :
##      original  bias    std. error
## t1* 0.9569567       0           0
plot(reps)
## [1] "All values of t* are equal to  0.956956664143508"

Após analisar este resultado podemos materializar uma das desvantagens do Bootstrap paramétrico. No histograma acima representado pelo Bootstrap não paramétrico, podemos observar a distribuição da nossa estatística alvo, o erro quadrado. Podemos concluir desse mesmo gráfico que esta estatística não possui uma distribuição normal. A função boot, da biblioteca boot do R, aqui usada assume por defeito para bootstraping paramétrico uma distribuição normal para a estatística a ser estudada. Dado isto o Bootstrap não consegue melhorar a estatística alvo, resultando num erro quadrado constante sempre igual ao dos dados originais.