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.2Entramos 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] 20t.m <- mean(tempo)
t.m## [1] 4.745t.v <- var(tempo)
t.v## [1] 0.9920789A 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.211157Podemos 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.382186Escrever 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.745Nesta 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.285058Portanto, 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.285058Podemos 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.285058t.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.745Carregue os dados do endereço abaixo
url <- "http://www.leg.ufpr.br/~fernandomayer/dados/crabs.csv"E faça o que se pede abaixo: