1 Breve introdução ao Rstudio e como instalar pacotes no Rstudio
- 1.1 O que é o Rstudio
- 1.2 Instalando pacotes
2 Estatística Descritiva
- 2.1 Medidas Resumo
3 Estatística Inferencial
4 Regressão Linear Simples e Múltipla
- 4.1 Formulação do Modelo
- 4.2 Exemplo prático

1 Breve introdução ao Rstudio e como instalar pacotes no Rstudio

1.1 O que é o Rstudio

O RStudio é um ambiente de desenvolvimento integrado (IDE) para R. Ele inclui um console(onde as saídas dos comandos são apresentadas), editor de sintaxe que suporta execução direta de código, bem como ferramentas para plotagem, histórico, depuração e gerenciamento de espaço de trabalho.

1.2 Instalando pacotes

Se sua necessidade de usar o R requer um pacote específico você pode seguir as instruções abaixo para instalar o pacote.

Abra o Rstudio
Clique na aba Package no seção do canto inferior direito e clique em install. A seguinte caixa de dialogo irá aparecer.
Na caixa de diálogo, escreva o nome do pacote que deseja instalar. Depois clique em install e espere o pacote ser instalado.
Para carregar o pacote, vá até a mesma aba do passo 1 e na caixa de procura digite o nome do pacote. Após achar o pacote clique na caixa para selecioná-lo e o pacote será carregado.

2 Estatística Descritiva

É uma das áreas mais importantes da estatística. Ela deve ser considera como ferramenta inicial para análise de qualquer conjunto de dados. Tem como objetivo sumarizar e descrever os dados.

Vamos utilizar o software R para iniciarmos o nosso estudo em estatística descritiva básica.

2.1 Medidas Resumo

Uma medida de tendência central (também conhecida como medidas de centro ou localização central) é uma medida resumo que tenta descrever um conjunto inteiro de dados com um único valor que representa o centro de sua distribuição.

Existem três principais medidas de tendência central: a moda, a mediana e a média. Cada uma dessas medidas descreve uma indicação diferente do valor típico ou central na distribuição dos dados.

Essas medidas são simples de serem calculadas dentro do R. Vamos utilizar o conjunto de dados abaixo para aplicar as funções presentes no R para calcular estas medidas.

conc	rate	state
0.02	76	treated
0.02	47	treated
0.06	97	treated
0.06	107	treated
0.11	123	treated
0.11	139	treated
0.22	159	treated
0.22	152	treated
0.56	191	treated
0.56	201	treated
1.10	207	treated
1.10	200	treated
0.02	67	untreated
0.02	51	untreated
0.06	84	untreated
0.06	86	untreated
0.11	98	untreated
0.11	115	untreated
0.22	131	untreated
0.22	124	untreated
0.56	144	untreated
0.56	158	untreated
1.10	160	untreated

Os dados tratam de um experimento para investigar a velocidade de uma determinada reação e a concentração do substrato em uma reação enzimática envolvendo células tratadas e não tratadas com puromicina.

O conjunto de dados apresenta 3 variáveis:

conc
rate
state

conc e rate são variáveis quantitativas. Já a variável state trata-se de uma variável qualitativa. Visto isso, vamos importar o conjunto de dados.

2.1.1 Importando dados no R

Neste caso, o conjunto de dados já está disponível na própria instalação do R. Para usá-lo basta utilizar o comando abaixo:

dados <- Puromycin
head(dados, 10) # Apresentando no console as 10 primeiras linhas do dataset

##    conc rate   state
## 1  0.02   76 treated
## 2  0.02   47 treated
## 3  0.06   97 treated
## 4  0.06  107 treated
## 5  0.11  123 treated
## 6  0.11  139 treated
## 7  0.22  159 treated
## 8  0.22  152 treated
## 9  0.56  191 treated
## 10 0.56  201 treated

Para importar arquivos que não estejam na base de dados do R, use o seguinte comando:

meus_dados <- read.csv("Digite o caminho do arquivo", sep = ",", dec = ".")

Os parâmetros sep e dec são respectivamente o tipo de parador usado no conjunto de dados para separar cada observação(vírgula) e o separador decimal para as variáveis quantitativas(ponto).

Dado que já importamos os dados. Vamos calcular algumas medidas de tendência central e dispersão.

2.1.2 Funções para medidas descritivas

Vamos calcular a média, mediana e amplitude(mínimo e máximo) das variáveis quantitativas do conjunto de dados. Para isso vamos usar os comandos abaixo:

# MEDIDAS DESCRITIVAS(GERAL) PARA A VARIAVEL CONC
mean(dados$conc)    # MEDIA

## [1] 0.3121739

median(dados$conc)  # MEDIANA

## [1] 0.11

range(dados$conc)   # AMPLITUDE

## [1] 0.02 1.10

min(dados$conc)     # MINIMO

## [1] 0.02

max(dados$conc)     # MAXIMO

## [1] 1.1

var(dados$conc)     # VARIANCIA

## [1] 0.1318632

sd(dados$conc)      # DESVIO PADRAO

## [1] 0.3631298

quantile(dados$conc, c(0.25,0.5,0.75)) # QUARTIS

##  25%  50%  75% 
## 0.06 0.11 0.56

# MEDIDAS DESCRITIVAS(GERAL) PARA A VARIAVEL RATE
mean(dados$rate)    # MEDIA

## [1] 126.8261

median(dados$rate)  # MEDIANA

## [1] 124

range(dados$rate)   # AMPLITUDE

## [1]  47 207

min(dados$rate)     # MINIMO

## [1] 47

max(dados$rate)     # MAXIMO

## [1] 207

var(dados$rate)     # VARIANCIA

## [1] 2257.514

sd(dados$rate)      # DESVIO PADRAO

## [1] 47.5133

quantile(dados$rate, c(0.25,0.5,0.75)) # QUARTIS

##   25%   50%   75% 
##  91.5 124.0 158.5

Para a variável categórica, podemos calcular a frequência de cada classe já que média e mediana não fazem sentido para esse tipo de informação. Para isso, use o comando abaixo:

table(dados$state)

## 
##   treated untreated 
##        12        11

Agora vamos calcular as mesmas estatísticas agora considerado cada um dos grupos definidos pela variável state. Para isso usaremos o seguinte comando:

tapply(dados$conc, dados$state, mean)

##   treated untreated 
## 0.3450000 0.2763636

Observe a ordem a qual cada parâmetro foi adicionado na função. Primeiro inserimos a variável de interesse(conc). Em seguite, inserimos a variável que define os grupos(state). E por último a medida descritiva(no caso do exemplo acima foi a média).

Podemos repetir isso para outras medidas descritivas

tapply(dados$conc, dados$state, median)

##   treated untreated 
##     0.165     0.110

tapply(dados$conc, dados$state, range)

## $treated
## [1] 0.02 1.10
## 
## $untreated
## [1] 0.02 1.10

tapply(dados$conc, dados$state, min)

##   treated untreated 
##      0.02      0.02

tapply(dados$conc, dados$state, max)

##   treated untreated 
##       1.1       1.1

tapply(dados$conc, dados$state, var)

##   treated untreated 
## 0.1589000 0.1126055

tapply(dados$conc, dados$state, sd)

##   treated untreated 
## 0.3986226 0.3355674

Outra importante estatística a ser calculada é o nível de associação linear entre duas variáveis. Para medirmos isso usamos o coeficiente de correlação. Para isso usaremos o comando abaixo:

cor(dados[,-3], method = "spearman")

##           conc      rate
## conc 1.0000000 0.9520149
## rate 0.9520149 1.0000000

Observe que a saída da função nos retorna uma matriz. Chamamos esta matriz de matriz de correlação. Nela apresentam-se as correlações entre as variáveis do conjunto de dados usado. Perceba também que os dados foram colocados dentro da função de forma diferente. Neste caso, retiramos a variável state do conjunto de dados(utilizando [,-3], ou seja, removendo a 3 coluna do conjunto de dados) visto para o cálculo do coeficiente de correlação, não deve ser usada variáveis de natureza qualitativa. Observe também o parâmetro method. Nele é definido o tipo de coeficiente de correlação que se deseja calcular(Pearson, Spearman e Kendall). Mais informações sobre as difernças e em que situações usar esses coeficientes, clique AQUI.

Tente fazer o mesmo para a variável rate!

2.1.3 Descritiva por gráficos

Muitas vezes a visualização de como os dados se apresentam em algum tipo de gráfico torna o entendimento de seu comportamento mais fácil. Um gráfico bastante usado em estatística é o boxplot.

Vamos inicialmente fazer um boxplot para a variável conc para cada um dos grupos(treated e untreated). Para isso, vamos utilizar o comando abaixo:

boxplot(dados$conc ~ dados$state, col="dodgerblue", main="Boxplot para variável conc")

A função boxplot recebe iniciamente como parâmetro a variável de interesse(conc) e em seguida a variável que define os grupos(state). O parâmetros col aplicar uma cor a escolha do usuário para os boxplots e o parâmetro main dá título ao boxplot.Link para uma tabela de cores AQUI!

Observe pelo tamanho das caixas que a variabilidade do grupo treated é maior do que a do grupo untreated. Resultado esse que verificamos mais acima utilizando a função tapply() juntamente com a função var() para calcularmos a variância desta variável para cada um dos grupos.

Replicando os comandos apresentados acima, vamos fazer o mesmo procedimento para a variável rate.

boxplot(dados$rate ~ dados$state, col="dodgerblue", main="Boxplot para variável rate")

Novamente voltando os resultados das medidas descritivas. Obtivemos um coeficiente de correlação de Pearson entre as variáveis conc e rate de 0.8055. Isso representa uma associação positiva de caráter forte entre as duas variáveis. Podemos verificar esta associação através de um ** gráfico de dispersão** ou scatter plot.

Para verificarmos esse comportamento graficamente, usaremos o comando abaixo:

plot(dados$conc, dados$rate, col="firebrick", pch=16)

Observe o comportamento crescente da variável rate em relação ao aumento da variável conc. Assim temos uma confirmação visual do nível de associação entre essas duas variáveis(como visto mais acima com o coeficiente de correlação).

Outro gráfico relevante em uma análise descritiva é o histograma. Com ele é possível verificar comportamentos como assimetria e curtose nos dados(Clique nos nomes para mais informações!).

Nativamente no R, não há função para calcularmos os coeficientes de assimetria e curtose mas podemos instalar o pacote fBasics que fornece funções para o cálculo desses coeficientes.

Na seção 1.2 vimos como instalar pacotes. Porém podemos usar um comando para instalar e carregar pacotes no R. Vamos instalar o pacote fBasics utilizando este comando.

install.packages("fBasics") # INSTALA O PACOTE
library(fBasics) # CARREGA O PACOTE PARA USO

Com o pacote carregado vamos utilizar as funções skewness() e kurtosis() para calcular respectivamente a assimetria e a curtose da variável conc do conjunto de dados.

skewness(dados$conc)

## [1] 1.205576
## attr(,"method")
## [1] "moment"

kurtosis(dados$conc)

## [1] 0.05829217
## attr(,"method")
## [1] "excess"

Observe que a saída da função retornou outras informações. Estas informações são referentes a como foram feito os cálculos para assimetria e curtose. Para o coeficiente de assimetria usou-se os momentos amostrais. E para curtose também foram usado os momentos amostrais porém o valor apresentado é o excesso de curtose em relação a distribuição Normal(Curtose da Normal é 3), isto é, a variável conc tem 0.05829217 a mais de curtose que a distribuição Normal.

Tente reproduzir tudo que foi feito até o momento e tirar conclusões sobre a variável rate deste conjunto de dados.

2.1.4 Tabelas de Contingência

Uma tabela de contingência calcula observações para mais de uma variável categórica. Suas linhas e colunas representam as variáveis categóricas do conjunto de dados. Vamos utilizar como exemplo prático a eleição para chefe de departamento entre dois candidatos. Uma pesquisa registrou o sexo e o voto de 100 eleitores escolhidos através de um esquema de amostragem aleatório e tabulou-se os dados conforme mostrado abaixo:

	Candidato A	Candidato B	TOTAL
HOMENS	28	20	48
MULHERES	39	13	52
TOTAL	67	33	100

Esta tabela de contingência correlaciona as respostas por sexo e voto. A contagem na intersecção da linha i e coluna j é identificada por \(n_{ij}\), e representa o número de observações que exibem essa combinação de níveis. Por exemplo, \(n_{1,2}\) exibe o número de entrevistados do sexo masculino que tenham votado para o Candidato B.

A tabela também inclui totais marginais para cada nível das variáveis. Os totais marginais para as linhas mostram que 52 dos entrevistados eram mulheres. Os totais marginais para colunas mostram que 67 entrevistados votaram no Candidato A. Além disso, o total global mostra que o tamanho amostral é de 100.

As tabelas de contingência também podem revelar associação entre as duas variáveis. O teste Qui-quadrado ou Teste Exato de Fisher podem determinar se as contagens observadas diferem significativamente das contagens esperadas sob a hipótese nula de não associação. Por exemplo, você pode testar se existe uma associação entre sexo e voto. Veremos isto mais a frente na seção 2.

As tabelas de contingência mais simples são tabelas com dois fatores que contam as respostas por duas variáveis. Você pode categorizar observações por três ou mais variáveis “cruzando” as informações que cada variável representa. No exemplo do voto, pode-se classificar as respostas, por exemplo, por situação empregatícia como é apresentado abaixo:

	Candidato A	Candidato B	TOTAL
HOMENS/EMPREGADOS	18	19	37
HOMENS/DESEMPREGADOS	10	1	11
MULHERES/EMPREGADAS	33	10	43
MULHERES/DESEMPREGADAS	6	3	9
TOTAL	67	33	100

Vamos produzir uma tabela de contigência para o banco de dados cancer.csv. Para isso vamos importar os dados conforme o comando abaixo:

dados <- read.csv("~/Trabalhos/Escola de Bioestatistica/cancer.csv")
head(dados, 10) # APRESENTANDO AS 10 PRIMEIRAS LINHAS DO CONJUNTO DE DADOS,

##    TREATMENT  STATE
## 1    Aspirin Cancer
## 2    Aspirin Cancer
## 3    Aspirin Cancer
## 4    Aspirin Cancer
## 5    Aspirin Cancer
## 6    Aspirin Cancer
## 7    Aspirin Cancer
## 8    Aspirin Cancer
## 9    Aspirin Cancer
## 10   Aspirin Cancer

                # PORÉM O CONJUNTO DE DADOS APRESENTA 39876 LINHAS.

Com os dados já importados, agora vamos criar nossa tabela de contigência para Tratamento \(\times\) Estado:

tabela <- table(dados$STATE, dados$TREATMENT)
tabela

##            
##             Aspirin Placebo
##   Cancer       1438    1427
##   No cancer   18496   18515

Vamos observar graficamente o resultado desta tabela de maneira gráfica. Para isso usaremos a função mosaicplot():

mosaicplot(t(tabela), col=c("firebrick", "dodgerblue"), cex.axis = 0.75, sub = "Tratamento", ylab = "Frequência Relativa", main = "")

Perceba que foi aplicada a função t() na tabela criada anteriormente. Esta função transpõe os dados. Ou seja, as linhas viram coluna e as colunas viram linha. O parâmetro col é a cor de cada variável, cex.axis é o tamanho da fonte dentro do gráfico, sub é a legenda do eixo X, ylab é a legenda do eixo Y, main é o título do gráfico que no exemplo foi deixado em branco.

2.1.4.1 Cálculo da Razão de chances

A Razão de Chances é usada para comparar as chances de dois eventos. A razão são as chances de sucesso dado que existe uma certa condição dividida pelas chances de sucesso dado que existe uma condição diferente. Você pode encontrar mais detalhes sobre isso AQUI.

Para o cálculo da razão de chances (Odds Ratio) iremos utilizar o pacote epitools. Instale e carregue este pacote da mesma forma que fizemos mais acima como pacote fBasics.

install.packages("epitools")
library(epitools)

Em alguns livros a fórmula para calcular a razão de chances usa o layout da tabela conforme o apresentado abaixo:

	TRATAMENTO	CONTROLE
DOENTE	\(X_{1,1}\)	\(X_{1,2}\)
NÃO-DOENTE	\(X_{2,1}\)	\(X_{2,2}\)

A tabela que criamos já se apresenta no layout apresentado assim. Então basta aplicarmos a função oddsratio().

oddsratio(tabela, method = "wald")

## $data
##            
##             Aspirin Placebo Total
##   Cancer       1438    1427  2865
##   No cancer   18496   18515 37011
##   Total       19934   19942 39876
## 
## $measure
##            odds ratio with 95% C.I.
##             estimate     lower    upper
##   Cancer    1.000000        NA       NA
##   No cancer 1.008744 0.9349043 1.088415
## 
## $p.value
##            two-sided
##             midp.exact fisher.exact chi.square
##   Cancer            NA           NA         NA
##   No cancer  0.8224348    0.8310911  0.8223986
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "Unconditional MLE & normal approximation (Wald) CI"

O parâmetro method é usado para indicar qual método será usado para calcular as razões e os intervalos de confiança. A função aceita 3 métodos: Median-Unbiased Estimation (midp), Estimação por Máxima Verossimilhança condicional (Fisher), Estimação por Máxima Verossimilhança não-condicional (Wald) e Small Sample Adjustment (Small).

2.1.4.2 Risco Relativo

Risco relativo é um termo estatístico usado para descrever o risco de um certo evento acontecer em um grupo contra o outro. Você pode encontrar mais detalhes sobre isso AQUI.

Para calcularmos o risco relativo, usaremos outra função contida no pacote epitools chamada riskratio(). Porém, precisamos alterar a nossa tabela para que ela nos retorne uma saída correta. O layout da tabela que a função aceita como entrada é o oposto da que geramos mais acima. Então vamos usar a função t() em nossa tabela para obtermos a tabela transposta. Para isso vamos criar outra tabela, como mostrado abaixo:

tabela_trans <- t(tabela)
tabela_trans

##          
##           Cancer No cancer
##   Aspirin   1438     18496
##   Placebo   1427     18515

Vamos comparar as duas tabelas para ratificarmos as diferenças:

tabela

##            
##             Aspirin Placebo
##   Cancer       1438    1427
##   No cancer   18496   18515

tabela_trans

##          
##           Cancer No cancer
##   Aspirin   1438     18496
##   Placebo   1427     18515

Observe que as linhas de tabela são as colunas da tabela_trans. Essa operação que acabamos de fazer chamasse transpor uma tabela.

Com a tabela já transposta podemos agora inserí-la na função riskratio() e obter os resultados.

riskratio(tabela_trans, method = "wald")

## $data
##          
##           Cancer No cancer Total
##   Aspirin   1438     18496 19934
##   Placebo   1427     18515 19942
##   Total     2865     37011 39876
## 
## $measure
##          risk ratio with 95% C.I.
##           estimate     lower    upper
##   Aspirin 1.000000        NA       NA
##   Placebo 1.000626 0.9951756 1.006106
## 
## $p.value
##          two-sided
##           midp.exact fisher.exact chi.square
##   Aspirin         NA           NA         NA
##   Placebo  0.8224348    0.8310911  0.8223986
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "Unconditional MLE & normal approximation (Wald) CI"

3 Estatística Inferencial

3.1 O que é estatística inferencial

A estatística inferencial é um dos dois principais ramos da estatística. Ela usa uma amostra aleatória de dados de uma população para descrever e fazer inferências sobre a população. Ela é valiosa quando analisar cada membro de uma população inteira não é conveniente ou possível. Por exemplo, medir a taxa de açucar no sangue de todos os brasileiros é algo impraticável. Você pode medir as taxas de açucar em uma amostra aleatória representativa da população brasileira. Você pode usar as informações da amostra para fazer generalizações sobre a população. Para informações mais detalhadas clique AQUI.

3.2 Teste \(\chi^2\) para independência

O teste de independência Qui-Quadrado é usado para descobrir se existe uma associação entre as variáveis em uma tabela de contingência construída à partir de observações da amostra. Para maiores detalhes sobre o teste clique AQUI.

Algumas condições devem ser atendidas para o uso deste teste. As condições são:

Preferencialmente deve-se ser aplicado em amostras grandes(\(n \geq 30\))
Independência entre as observações.
Não pode haver frequências inferiores a 1.

As hipóteses deve ser montada da seguinte forma:

\[H_0: \text{As variáveis não são associadas}\] \[H_1: \text{As variáveis são associadas}\]

O teste se dá através do cálculo da estatística:

\[\chi^2 = \sum_{i,j}\frac{(f_{ij}-e_{ij})^2}{e_{ij}}\]

Onde \(f_{ij}\) é o vetor com de contagens observadas de cada classe e \(e_{ij}\) é o vetor que representa a frequência esperada de cada classe.

Vamos agora utilizar o R para exemplificar este teste.

Existe um conjunto de dados no R chamado survey. A coluna Smoke armazena os estudantes que que possuem o hábito de fumar, enquanto a coluna Exer indica o nível de atividade física de cada um dos estudantes. Vamos carregar as informações e criar uma tabela de contingência para este conjunto de dados.

library(MASS)
dados <- survey
tabela_1 <- table(dados$Smoke, dados$Exer)
tabela_1

##        
##         Freq None Some
##   Heavy    7    1    3
##   Never   87   18   84
##   Occas   12    3    4
##   Regul    9    1    7

Agora vamos testar a hipótese de que o hábito de fumar entre os estudantes da nossa amostra não está associado ao seus nível de atividade física a um nível de significância de 5%.

Para isso basta aplicarmos a função chisq.test():

chisq.test(tabela_1)

## 
##  Pearson's Chi-squared test
## 
## data:  tabela_1
## X-squared = 5.4885, df = 6, p-value = 0.4828

O p-valor calculado pela função, 0.4828, é maior que o nível de significância especificado no teste. Logo não existem indícios para rejeitarmos a nossa hipótese nula de que o hábito de fumar dos estudantes não esteja associado a seus níveis de atividade física.

3.3 Teste t para uma amostra

O teste t de uma amostra é usado para comparar a média de uma amostra com uma média conhecida (ou teórica/hipotética). Para maiores detalhes sobre o teste, clique AQUI.

O teste é baseado na seguinte estatística:

\[t = \frac{m-\mu}{\dfrac{s}{\sqrt{n}}},\]

onde:

m: média amostral
n: tamanho amostral
s: desvio padrão amostral
\(\mu\): média teórica/hipotética

Algumas condições devem ser verificadas antes de se aplicar o teste:

A amostra deve ser suficientemente grande (\(n>30\))
Caso a amostra seja menos, deve-se verificar sua normalidade através de testes.

Vamos utilizar o conjunto de dados bottles.csv para exemplificarmos o teste utilizando o R. Novamente, vamos importar os dados para iniciarmos o teste:

dados <- read.csv("~/Trabalhos/Escola de Bioestatistica/bottles.csv")

Este conjunto de dados é referente ao seguinte problema. Em uma fábrica de produtos farmacêuticos, existe uma linha de produção de frascos de soro fisiológico com volume de 500ml. Desconfia-se de que a máquina responsável pelo preenchimento dos frascos esteja desregulada e que os frascos estejam saindo com menos de 500ml de soro. Uma amostra de 20 frascos foi coletada e seus volumes avaliados. Será que realmente os frascos estão saindo da linha de produção com menos de 500ml?

Vamos formular as hipóteses:

\[H_0: \text{Média do volume dos frascos é igual a 500ml}\] \[H_1: \text{Média do volume dos frascos é menor que 500ml}\]

Observe que temos apenas a observação do volume de 20 frascos. O que nos dá um tamanho amostral de 20. O que, a priori, fere um dos pré-requisitos para a aplicação do teste. Porém, neste caso, vamos verificar a normalidade dos dados através do teste de Shapiro-Wilks e com a função qqPlot() do pacote car(instale e carregue este pacote para usar a função). Para isso vamos usar o comando:

shapiro.test(dados$Volume)

## 
##  Shapiro-Wilk normality test
## 
## data:  dados$Volume
## W = 0.97329, p-value = 0.8223

qqPlot(dados$Volume)

## [1] 12 19

A um nível de significância de 5%, podemos considerar que os dados seguem distribuição Normal (p-valor > 0.05).

Agora vamos calcular o volume médio da amostra retirada.

mean(dados$Volume)

## [1] 491.5705

Aparentemente a máquina mostra-se desregulada mas não podemos tirar conclusões ainda. Vamos aplicar o teste-t para esta amostra utilizando a função t.test(). Vamos utilizar um nível de significância de 5% (ou nível de confiança de 95%).

t.test(dados$Volume, mu=500, alternative="less", conf.level=0.95)

## 
##  One Sample t-test
## 
## data:  dados$Volume
## t = -1.5205, df = 19, p-value = 0.07243
## alternative hypothesis: true mean is less than 500
## 95 percent confidence interval:
##      -Inf 501.1569
## sample estimates:
## mean of x 
##  491.5705

Utilizando-se um nível de significância de 5%, não temos indícios suficientes para rejeitar a hipótese nula de que o volume médio dos frascos produzidos seja 500ml.

3.4 Teste-t para duas amostras independentes

Teste-t para duas amostras independentes ou Teste-t não-pareado para duas amostras é auto-explicativo. Ele testa hipóteses sobre a média de duas amostras independentes. Para mais detalhes e informações sobre o teste, clique AQUI.

Neste teste temos que considerar alguns pontos:

Independencia entre os grupos.
A normalidade dos grupos.
A homogeneidade dos grupos (variâncias iguais).

Para verificar a normalidade temos o teste de Shapiro-Wilks e para verificar a variância dos grupos testemos o Teste-F.

Sabendo que as suposições de independencia e normalidade são atendidas, temos duas opções para verificar suposições sobre a média dos grupos:

Quando as variâncias dos grupos são iguais.
Quando as variâncias dos grupos não são iguais.

Para cada situação acima apresentada, temos diferenças no cálculo da estatística para o teste. Você pode verificar como as estatísticas são calculadas em cada caso neste link.

Vamos utilizar o exemplo onde queremos testar se a média do peso em um grupo composto por 9 homens e 9 mulheres são diferentes ou não. Assim, nossas hipóteses são:

\[H_0: \mu_{homem} = \mu_{mulher}\] \[H_0: \mu_{homem} \neq \mu_{mulher}\]

Vamos carregar os dados e verificar inicialmente a normalidade dos dados e a homogeneidade dos grupos, tudo isso a um nível de significância de 5% (visto que o quesito independência já está explicito no experimento. O peso de um homem não influencia no peso de uma mulher).

dados <- read.csv("~/Trabalhos/Escola de Bioestatistica/peso_sexo.csv")

Peso	Sexo
38.9	Masculino
61.2	Masculino
73.3	Masculino
21.8	Masculino
63.4	Masculino
64.6	Masculino
48.4	Masculino
48.8	Masculino
48.5	Masculino
67.8	Feminino
60.0	Feminino
63.4	Feminino
76.0	Feminino
89.4	Feminino
73.3	Feminino
67.3	Feminino
61.3	Feminino
62.4	Feminino

Varificando normalidade para os dois grupos.

tapply(dados$Peso, dados$Sexo, shapiro.test)

## $Feminino
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.86425, p-value = 0.1066
## 
## 
## $Masculino
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.94266, p-value = 0.6101

Podemos considerar que os grupos são oriundos de uma distribuição Normal.

Agora vamos verificar a homogeneidade dos grupos. Ou seja, se a variância nos grupos são iguais. Para isso usaremos a função var.test().

var.test(dados$Peso ~ dados$Sexo, data = dados)

## 
##  F test to compare two variances
## 
## data:  dados$Peso by dados$Sexo
## F = 0.36134, num df = 8, denom df = 8, p-value = 0.1714
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.08150656 1.60191315
## sample estimates:
## ratio of variances 
##          0.3613398

Como o p-valor > 0.05, então não temos indícios para rejeitar a hipótese nula de que as variâncias de cada grupo sejam iguais.

Finalmente, iremos aplicar o teste-t para duas amostras.

t.test(dados$Peso ~ dados$Sexo, data = dados, var.equal=TRUE)

## 
##  Two Sample t-test
## 
## data:  dados$Peso by dados$Sexo
## t = 2.7842, df = 16, p-value = 0.01327
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   4.029759 29.748019
## sample estimates:
##  mean in group Feminino mean in group Masculino 
##                68.98889                52.10000

Como o p-valor < 0.05, existem indícios suficientes para rejeitar a hipótese nula de que o peso médio dos homens é diferente do peso médio das mulheres.

4 Regressão Linear Simples e Múltipla

A regressão linear (ou modelo linear) é usada para prever uma variável quantitativa \(Y\) com base em uma ou várias(caso múltiplo) variáveis preditoras X. (James et al., 2014, P. Bruce e Bruce, 2017).

O objetivo é construir um modelo matemático que a variável resposta \(Y\) em função das variáveis preditoras (ou explicativas) \(X\). Dessa forma, uma vez construído um modelo estatísticamente significante, é possível usá-lo para fazer previsões futuras com base em novos valores de \(X\).

Quando você cria um modelo de regressão, é necessário avaliar o desempenho do modelo preditivo. Em outras palavras, você precisa avaliar o quão bem o modelo prevê o resultado de novos dados de teste que não foram usados para construir o modelo.

Duas métricas importantes são comumente usadas para avaliar o desempenho do modelo de regressão:

Erro quadrátio médio (EQM), que mede o erro de previsão do modelo. Corresponde à diferença média entre os valores observados e o valor previsto pelo modelo.
\(R^2\) ajustado, representando a correlação quadrática entre os valores observados e os valores previstos pelo modelo. Quanto maior o \(R^2\), melhor o modelo.

4.1 Formulação do Modelo

A fórmula matemática para os modelos de regressão linear simples é definida por:

\[y_i = \beta_0 + \beta_1 x_i + \epsilon_i, \\ i=\{1,2,...,n\},\]

onde n é a quantidade de observações.

Quando se tem mais de uma variável preditora, o modelo é definido por:

\[y_i = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + ... + \beta_j x_{ji} + \epsilon_{ji}, \\ i=\{1,2,...,n\},\\ j=\{1,2,..., m\},\]

onde n é a quantidade de observações e m é a quantidade variáveis preditoras.

Grande parte da teoria sobre modelos de regressão linear simples e múltiplos será omitida nesta seção, pois o intuito deste documento é gerar um guia para o uso de algumas ferramentas estatísticas sob o apoio computacional do software R. Para mais informações sobre toda a teoria que envolve os modelos de regressão linear simples e múltipla, verifique este livro.

4.2 Exemplo prático

O conjunto de dados publicidade.csv contêm o impacto de três mídias de publicidade (YouTube, Facebook e Jornal Convencional) nas vendas. Os dados são o orçamento de publicidade em milhares de dólares, juntamente com as vendas. O experimento de publicidade foi repetido 200 vezes.

Inicialmente, vamos importar e verificar os dados.

dados <- read.csv("~/Trabalhos/Escola de Bioestatistica/publicidade.csv")

youtube	facebook	newspaper	sales
276.12	45.36	83.04	26.52
53.40	47.16	54.12	12.48
20.64	55.08	83.16	11.16
181.80	49.56	70.20	22.20
216.96	12.96	70.08	15.48
10.44	58.68	90.00	8.64
69.00	39.36	28.20	14.16
144.24	23.52	13.92	15.84
10.32	2.52	1.20	5.76
239.76	3.12	25.44	12.72
79.32	6.96	29.04	10.32
257.64	28.80	4.80	20.88
28.56	42.12	79.08	11.04
117.00	9.12	8.64	11.64
244.92	39.48	55.20	22.80
234.48	57.24	63.48	26.88
81.36	43.92	136.80	15.00
337.68	47.52	66.96	29.28
83.04	24.60	21.96	13.56
176.76	28.68	22.92	17.52

Acima foram apresentadas as 20 primeiras linhas do conjunto de dados. Como dito anteriormente, o experimento foi realizado 200 vezes.

Inicialmente, vamos identificar as variáveis regressoras e a variável resposta:

Variável Resposta: sales
Variáveis Regressoras: youtube, facebook, newspaper

Outra verificação que devemos fazer é se existe realmente alguma relação linear entre as variáveis regressoras e a variável resposta. Para isso vamos usar o gráfico que apresenta essa correlação.

pairs(dados, col="blue", pch=16)

Observe que existe, mesmo que de maneira moderada para a variável newspaper, uma correlação linear entre sales com as demais variáveis. Assim é plausível usar um modelo de regressão linear neste conjunto de dados.

Inicialmente vamos escolher apenas a variável facebook para modelar uma regressão linear simples. Para isso usaremos o seguinte comando:

reg_simples <- lm(sales ~ facebook, data = dados)
summary(reg_simples)

## 
## Call:
## lm(formula = sales ~ facebook, data = dados)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -18.8766  -2.5589   0.9248   3.3330   9.8173 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 11.17397    0.67548  16.542   <2e-16 ***
## facebook     0.20250    0.02041   9.921   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.13 on 198 degrees of freedom
## Multiple R-squared:  0.332,  Adjusted R-squared:  0.3287 
## F-statistic: 98.42 on 1 and 198 DF,  p-value: < 2.2e-16

plot(dados$facebook, dados$sales, col="blue", pch=16)
abline(reg_simples, lwd=2, col="red")

Com as informações que obtivemos com a saída da função lm(), temos que o nosso modelo é:

\[\hat{y}_{sales} = 11.17397 + 0.20250 x_{facebook}\]

Observe também que os dois coeficientes são extremamente significativos neste modelo pois o p-valor calculado para os dois foi algo muito próximo de 0 (\(2.10^{-16}\)). O \(R^2\) ajustado calculado para o modelo foi de 0.3287. Em geral esse valor é baixo e pode ser explicado pelo espalhamento dos dados no gráfico dispersão entre a variável facebook e sales. Perceba que os dados estão bastante “espalhados” no gráfico, indicando uma variabilidade muito grande. Como o objetivo da reta regressora é passar o mais próximo de todos os pontos, quanto maior a dispersão dos dados menor será a adequabilidade do modelo.

Agora vamos fazer um modelo de regressão linear múltipla utilizando os mesmos dados.

reg_mult <- lm(sales ~ youtube + facebook + newspaper, data = dados)
summary(reg_mult)

## 
## Call:
## lm(formula = sales ~ youtube + facebook + newspaper, data = dados)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.5932  -1.0690   0.2902   1.4272   3.3951 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.526667   0.374290   9.422   <2e-16 ***
## youtube      0.045765   0.001395  32.809   <2e-16 ***
## facebook     0.188530   0.008611  21.893   <2e-16 ***
## newspaper   -0.001037   0.005871  -0.177     0.86    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.023 on 196 degrees of freedom
## Multiple R-squared:  0.8972, Adjusted R-squared:  0.8956 
## F-statistic: 570.3 on 3 and 196 DF,  p-value: < 2.2e-16

Observe agora que temos 3 variáveis preditoras, porém apenas 2 se mostram significativas no modelo (youtube, facebook - p-valor = \(2.10^{-16}\) e newspaper - p-valor = 0.86). Perceba também que o valor do \(R^2\) ajustado é alto (0.8956) o que é um indício de que nosso modelo está bem ajustado. Já que vimos que a variável newspaper não é significativa em nosso modelo vamos retirá-la e recalcular o modelo.

reg_mult <- lm(sales ~ youtube + facebook, data = dados)
summary(reg_mult)

## 
## Call:
## lm(formula = sales ~ youtube + facebook, data = dados)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.5572  -1.0502   0.2906   1.4049   3.3994 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.50532    0.35339   9.919   <2e-16 ***
## youtube      0.04575    0.00139  32.909   <2e-16 ***
## facebook     0.18799    0.00804  23.382   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.018 on 197 degrees of freedom
## Multiple R-squared:  0.8972, Adjusted R-squared:  0.8962 
## F-statistic: 859.6 on 2 and 197 DF,  p-value: < 2.2e-16

Perceba agora que o \(R^2\) ajustado mudou muito pouco. O que significa que a variável newspaper é de pouca relevância na construção do modelo.