Nesta sessão vamos verificar como utilizar o R para obter intervalos de confiança para parâmetros de distribuições de probabilidade.
Para fins didáticos, mostrando os recursos do R, vamos mostrar três possíveis soluções:
Considere o seguinte problema:
O tempo de reação de um novo medicamento pode ser considerado como tendo distribuição Normal e deseja-se fazer inferência sobre a média que é desconhecida obtendo um intervalo de confiança. Vinte pacientes foram sorteados e tiveram seu tempo de reação anotado. Os dados foram os seguintes (em minutos):
2.9 3.4 3.5 4.1 4.6 4.7 4.5 3.8 5.3 4.9
4.8 5.7 5.8 5.0 3.4 5.9 6.3 4.6 5.5 6.2
Entramos com os dados com o comando
tempo <- c(2.9, 3.4, 3.5, 4.1, 4.6, 4.7, 4.5, 3.8, 5.3, 4.9,
4.8, 5.7, 5.8, 5.0, 3.4, 5.9, 6.3, 4.6, 5.5, 6.2)
Sabemos que o intervalo de confiança para média de uma distribuição normal com variância desconhecida, para uma amostra de tamanho \(n\) é dado por:
\[ \left(\bar{x} - t_t \sqrt{\frac{s^2}{n}} \quad, \quad \bar{x} + t_t \sqrt{\frac{s^2}{n}} \right) \]
onde \(t_t\) é o quantil de ordem \(1-\alpha/2\) da distribuição \(t\) de Student, com \(n-1\) graus de liberdade.
Vamos agora obter a resposta das três formas diferentes mencionadas acima.
Nos comandos a seguir calculamos o tamanho da amostra, a média e a variância amostral.
n <- length(tempo)
n
## [1] 20
t.m <- mean(tempo)
t.m
## [1] 4.745
t.v <- var(tempo)
t.v
## [1] 0.9920789
A seguir montamos o intervalo utilizando os quantis da distribuição \(t\), para obter um IC a 95% de confiança.
t.ic <- t.m + qt(c(0.025, 0.975), df = n - 1) * sqrt(t.v/length(tempo))
t.ic
## [1] 4.278843 5.211157
Podemos generalizar a solução acima agrupando os comandos em uma função. Nos comandos abaixo, primeiro definimos a função, inclusive adicionando um argumento conf
que permitirá que o usuário especifique o nível de confiança desejado
ic.m <- function(x, conf = 0.95){
n <- length(x)
media <- mean(x)
variancia <- var(x)
tinf <- (1 - conf)/2
tsup <- 1 - (1-conf)/2
quantis <- qt(c(tinf, tsup), df = n - 1)
ic <- media + quantis * sqrt(variancia/n)
return(ic)
}
A seguir, podemos aplicar nossa função ao objeto tempo
, calculando os intervalos com 95% e 99% de confiança
## IC com 95% de confiança
ic.m(tempo)
## [1] 4.278843 5.211157
## IC com 99% de confiança
ic.m(tempo, conf = 0.99)
## [1] 4.107814 5.382186
Escrever uma função é particularmente útil quando um procedimento será utilizado várias vezes.
t.test()
Mostramos as soluções acima para ilustrar a flexibilidade e o uso do programa. Entretanto não precisamos fazer isto na maioria das vezes porque o R já vem com várias funções para procedimentos estatísticos já escritas. Neste caso a função t.test()
realiza o teste \(t\) para a média, mas como resultado ela também retorna o intervalo de confiança. Note que as saídas abaixo coincidem com os obtidos anteriormente.
## IC com 95% de confiança
t.test(tempo)
##
## One Sample t-test
##
## data: tempo
## t = 21.305, df = 19, p-value = 1.006e-14
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 4.278843 5.211157
## sample estimates:
## mean of x
## 4.745
## IC com 99% de confiança
t.test(tempo, conf.level = 0.99)
##
## One Sample t-test
##
## data: tempo
## t = 21.305, df = 19, p-value = 1.006e-14
## alternative hypothesis: true mean is not equal to 0
## 99 percent confidence interval:
## 4.107814 5.382186
## sample estimates:
## mean of x
## 4.745
Nesta sessão vamos verificar como utilizar o R para fazer testes de hipótese sobre parâmetros de distribuições para as quais os resultados são bem conhecidos.
Os comandos e cálculos são bastante parecidos com os vistos em intervalos de confiança e isto nem poderia ser diferente visto que intervalos de confiança e testes de hipótese são relacionados.
Assim como fizemos com intervalos de confiança, aqui sempre que possível, e para fins didáticos mostrando os recursos do R, vamos mostrar três possíveis soluções:
Considerando o problema anterior, do tempo de reação de um medicamento, faça um teste de hipótse para avaliar se o tempo médio é de 4,5 minutos.
Com isso, já temos as hipótese nula e alternativa para a verdadeira média \(\mu\) do tempo de reação:
\[ \begin{align} \text{H}_0: \mu = 4,5 \\ \text{H}_1: \mu \neq 4,5 \\ \end{align} \]
Novamente vamos resolver esse problema de três maneiras diferentes.
Primeiro espacificamos o nível de significância \(\alpha = 0,05\), e achamos o valor crítico de \(t\). Depois calculamos a estatística de teste, e o p-valor associado.
## Usando alpha = 0,05, achamos o valor crítico de t com
alpha <- 0.05
tcrit <- abs(qt(alpha/2, df = n - 1))
tcrit
## [1] 2.093024
## Estatística de teste
tcalc <- (mean(tempo) - 4.5)/(sd(tempo)/sqrt(n))
tcalc
## [1] 1.100039
## tcalc > tcrit?
tcalc > tcrit
## [1] FALSE
## Calculano o p-valor
## Note que essa probabilidade deve ser multiplicada por 2, pois a
## função pt() só calcula uma cauda, e o teste é bicaudal
pt(tcalc, df = n - 1, lower.tail = FALSE) * 2
## [1] 0.285058
Portanto, por esse resultado, a hipótese nula não pode ser rejeitada.
Podemos agregar os comandos acima e algumas definições para criar nossa própria função para fazer o teste \(t\). Siga os passos nos procedimentos gerais descritos acima e tente criar sua própria função.
## Função para fazer o teste t
teste.t <- function(x, mu, conf = 0.95){
n <- length(x)
media <- mean(x)
desvpad <- sd(x)
tcrit <- abs(qt((1 - conf)/2, df = n - 1))
tcalc <- (media - mu)/(desvpad/sqrt(n))
pvalor <- pt(tcalc, df = n - 1, lower.tail = FALSE) * 2
saida <- list("t critico" = tcrit,
"t calculado" = tcalc,
"tcalc > tcrit" = tcalc > tcrit,
"p-valor" = pvalor)
return(saida)
}
Aplicando a função, temos como resultado
teste.t(x = tempo, mu = 4.5)
## $`t critico`
## [1] 2.093024
##
## $`t calculado`
## [1] 1.100039
##
## $`tcalc > tcrit`
## [1] FALSE
##
## $`p-valor`
## [1] 0.285058
Podemos alterar o nível de significância \(\alpha\) para 0,01 por exemplo
teste.t(x = tempo, mu = 4.5, conf = 0.99)
## $`t critico`
## [1] 2.860935
##
## $`t calculado`
## [1] 1.100039
##
## $`tcalc > tcrit`
## [1] FALSE
##
## $`p-valor`
## [1] 0.285058
t.test()
Como vimos anteriormente, a função t.test()
realiza o teste \(t\), da mesma forma que programamos acima. Veja os resultados da saída dessa função e compare com os que você calculou acima.
## Teste t com alpha = 0,05
t.test(tempo, mu = 4.5)
##
## One Sample t-test
##
## data: tempo
## t = 1.1, df = 19, p-value = 0.2851
## alternative hypothesis: true mean is not equal to 4.5
## 95 percent confidence interval:
## 4.278843 5.211157
## sample estimates:
## mean of x
## 4.745
## Teste t com alpha = 0,01
t.test(tempo, mu = 4.5, conf.level = 0.99)
##
## One Sample t-test
##
## data: tempo
## t = 1.1, df = 19, p-value = 0.2851
## alternative hypothesis: true mean is not equal to 4.5
## 99 percent confidence interval:
## 4.107814 5.382186
## sample estimates:
## mean of x
## 4.745
Carregue os dados do endereço abaixo
url <- "http://www.leg.ufpr.br/~fernandomayer/dados/crabs.csv"
E faça o que se pede abaixo: