Funções do R para teste de hipóteses III
Mauro Campos
15 de maio de 2023
Referências para leitura deste material:
- P. Morettin & W. Bussab. Estatística Básica. 6a ed. São Paulo: Saraiva, 2010.
- M. Magalhães & AC. Lima. Noções de Probabolidade e Estatística Básica. 5a ed. São Paulo: EDUSP, 2002.
- M. Campos. Introdução ao R. Vitória, 2018. Disponível em: https://maurocmcampos.github.io/ec/html/intror.html.
- M. Campos. Funções do R para estatística descritiva. Vitória, 2023. Disponível em: https://maurocmcampos.github.io/ec/html/edR.html.
- M. Campos. Funções do R para teste de hipóteses I. Vitória, 2023. Disponível em: https://maurocmcampos.github.io/ec/html/th1R.html.
1 Testes qui-quadrado
1.1 Teste de aderência
Problema. Seja \(X\) uma variável aleatória (va.a.) que representa o número de impactos anteriores a falha de um equipamento eletrônico. Uma amostra aleatória simples de 80 observações de \(X\) foi obtida e os dados estão apresentados no script abaixo. Deseja-se realizar o seguinte teste:
\[ \begin{cases} H_{0}: & \mbox{$X$ segue o modelo geométrico com parâmetro $p=0.4$} \\ H_{1}: & \mbox{$X$ segue outra distribuição.} \end{cases} \]
A função chisq.test() pode ser usada para realizar este
teste de aderência.
df0 <- read.table(file="impactos.txt",header=TRUE)
str(df0)
## 'data.frame': 80 obs. of 1 variable:
## $ x: int 0 1 0 1 0 0 1 1 0 0 ...
#
# Dados
#
df0$x
## [1] 0 1 0 1 0 0 1 1 0 0 1 1 0 0 5 1 0 1 5 1 7 1 0 2 3 1 6 4 0 3 0 2 0 0 2 0 3 0
## [39] 1 0 2 1 1 1 1 2 1 4 4 0 0 4 1 0 4 3 0 2 0 0 1 2 0 1 1 0 0 0 1 1 1 2 0 2 1 2
## [77] 1 0 0 3
table(df0$x)
##
## 0 1 2 3 4 5 6 7
## 30 26 10 5 5 2 1 1
classes <- names(table(df0$x))
freq <- as.vector(table(df0$x))
freq.r <- round(100*as.vector(prop.table(freq)),2)
freq.a <- cumsum(freq.r)
df <- data.frame(classe=c(classes,"Total"),
frequencia.absoluta=c(freq,sum(freq)),
freq.relativa=c(freq.r,sum(freq.r)),
freq.acumulada=c(freq.a,NaN))
df
## classe frequencia.absoluta freq.relativa freq.acumulada
## 1 0 30 37.50 37.50
## 2 1 26 32.50 70.00
## 3 2 10 12.50 82.50
## 4 3 5 6.25 88.75
## 5 4 5 6.25 95.00
## 6 5 2 2.50 97.50
## 7 6 1 1.25 98.75
## 8 7 1 1.25 100.00
## 9 Total 80 100.00 NaN
freq.o <- df$frequencia.absoluta[-dim(df)[1]]
vec.p <- c(dgeom(0:6,prob=0.4),pgeom(6,prob=0.4,lower.tail=FALSE))
freq.e <- sum(freq.o)*vec.p
data.frame(classe=c("0","1","2","3","4","5","6",">6"),freq.o,freq.e,vec.p)
## classe freq.o freq.e vec.p
## 1 0 30 32.000000 0.4000000
## 2 1 26 19.200000 0.2400000
## 3 2 10 11.520000 0.1440000
## 4 3 5 6.912000 0.0864000
## 5 4 5 4.147200 0.0518400
## 6 5 2 2.488320 0.0311040
## 7 6 1 1.492992 0.0186624
## 8 >6 1 2.239488 0.0279936
ni <- c(freq.o[1:4],sum(freq.o[5:8]));n <- sum(ni)
prob <- c(dgeom(0:3,prob=0.4),pgeom(3,prob=0.4,lower.tail=FALSE))
data.frame(i=1:5,ni,npi=n*prob,prob)
## i ni npi prob
## 1 1 30 32.000 0.4000
## 2 2 26 19.200 0.2400
## 3 3 10 11.520 0.1440
## 4 4 5 6.912 0.0864
## 5 5 9 10.368 0.1296
out <- chisq.test(x=ni,p=prob);out
##
## Chi-squared test for given probabilities
##
## data: ni
## X-squared = 3.4433, df = 4, p-value = 0.4866
names(out)
## [1] "statistic" "parameter" "p.value" "method" "data.name" "observed"
## [7] "expected" "residuals" "stdres"
out$statistic
## X-squared
## 3.443287
out$parameter
## df
## 4
out$p.value
## [1] 0.4865539
out$method
## [1] "Chi-squared test for given probabilities"
out$data.name
## [1] "ni"
out$observed
## [1] 30 26 10 5 9
out$expected
## [1] 32.000 19.200 11.520 6.912 10.368
q <- qchisq(p=0.05,df=out$parameter,lower.tail=FALSE);q
## [1] 9.487729Não rejeitamos \(H_{0}\), ou seja, existem evidências experimentais que \(X\) segue o modelo geométrico com \(p=0.4\).
Problema. Seja \(X\) uma va.a. que representa a porcentagem de cinzas contidas no carvão produzido por uma certa empresa. Uma amostra aleatória simples de 250 observações de \(X\) foi obtida e os dados estão apresentados no script abaixo. Deseja-se realizar o seguinte teste:
\[ \begin{cases} H_{0}: & \mbox{$X$ segue o modelo normal} \\ H_{1}: & \mbox{$X$ segue outra distribuição.} \end{cases} \]
df0 <- read.table(file="cinzas.txt",header=TRUE)
str(df0)
## 'data.frame': 250 obs. of 1 variable:
## $ x: num 13 14.6 15.5 15 13.4 ...
#
# Dados
#
head(df0)
## x
## 1 13.029
## 2 14.553
## 3 15.548
## 4 15.010
## 5 13.390
## 6 16.161
tail(df0)
## x
## 245 15.252
## 246 15.402
## 247 13.278
## 248 11.789
## 249 12.788
## 250 12.015
xd <- cut(df0$x,breaks=seq(9.5,19.5,1),right=FALSE)
table(xd)
## xd
## [9.5,10.5) [10.5,11.5) [11.5,12.5) [12.5,13.5) [13.5,14.5) [14.5,15.5)
## 2 5 16 42 69 51
## [15.5,16.5) [16.5,17.5) [17.5,18.5) [18.5,19.5)
## 32 23 9 1
classes <- names(table(xd))
freq <- as.vector(table(xd))
freq.r <- round(100*as.vector(prop.table(freq)),2)
freq.a <- cumsum(freq.r)
df <- data.frame(classe=c(classes,"Total"),
frequencia.absoluta=c(freq,sum(freq)),
freq.relativa=c(freq.r,sum(freq.r)),
freq.acumulada=c(freq.a,NaN))
df
## classe frequencia.absoluta freq.relativa freq.acumulada
## 1 [9.5,10.5) 2 0.8 0.8
## 2 [10.5,11.5) 5 2.0 2.8
## 3 [11.5,12.5) 16 6.4 9.2
## 4 [12.5,13.5) 42 16.8 26.0
## 5 [13.5,14.5) 69 27.6 53.6
## 6 [14.5,15.5) 51 20.4 74.0
## 7 [15.5,16.5) 32 12.8 86.8
## 8 [16.5,17.5) 23 9.2 96.0
## 9 [17.5,18.5) 9 3.6 99.6
## 10 [18.5,19.5) 1 0.4 100.0
## 11 Total 250 100.0 NaN
dim(df)
## [1] 11 4
freq.o <- df$frequencia.absoluta[-dim(df)[1]]
aux <- seq(10,19,1)
x.bar <- sum(freq.o*aux)/sum(freq.o);x.bar
## [1] 14.512
s2 <- sum(freq.o*(aux-x.bar)^(2))/(sum(freq.o)-1)
s <- sqrt(s2);s
## [1] 1.643368
vec.p <- NULL
r <- length(aux)
for(k in 1:r){
if(k==1){vec.p[k] <- pnorm(aux[k]+0.5,x.bar,s)}
if((k>1)&(k<r)){vec.p[k] <- pnorm(aux[k]+0.5,x.bar,s)-pnorm(aux[k]-0.5,x.bar,s)}
if(k==r){vec.p[k] <- 1-sum(vec.p[1:(k-1)])}
}
freq.e <- sum(freq.o)*vec.p
data.frame(classe=classes,freq.o,freq.e,vec.p)
## classe freq.o freq.e vec.p
## 1 [9.5,10.5) 2 1.829171 0.007316685
## 2 [10.5,11.5) 5 6.524473 0.026097891
## 3 [11.5,12.5) 16 19.250539 0.077002156
## 4 [12.5,13.5) 42 39.648433 0.158593733
## 5 [13.5,14.5) 69 57.019114 0.228076454
## 6 [14.5,15.5) 51 57.265283 0.229061130
## 7 [15.5,16.5) 32 40.164249 0.160656996
## 8 [16.5,17.5) 23 19.669870 0.078679479
## 9 [17.5,18.5) 9 6.724386 0.026897545
## 10 [18.5,19.5) 1 1.904483 0.007617930
ni <- c(sum(freq.o[1:2]),freq.o[3:8],sum(freq.o[9:10]));n <- sum(ni)
prob <- c(sum(vec.p[1:2]),vec.p[3:8],sum(vec.p[9:10]))
data.frame(i=1:8,ni,npi=n*prob,prob)
## i ni npi prob
## 1 1 7 8.353644 0.03341458
## 2 2 16 19.250539 0.07700216
## 3 3 42 39.648433 0.15859373
## 4 4 69 57.019114 0.22807645
## 5 5 51 57.265283 0.22906113
## 6 6 32 40.164249 0.16065700
## 7 7 23 19.669870 0.07867948
## 8 8 10 8.628869 0.03451548
out <- chisq.test(x=ni,p=prob)
out
##
## Chi-squared test for given probabilities
##
## data: ni
## X-squared = 6.5518, df = 7, p-value = 0.477
out$observed
## [1] 7 16 42 69 51 32 23 10
out$expected
## [1] 8.353644 19.250539 39.648433 57.019114 57.265283 40.164249 19.669870
## [8] 8.628869
pvalor <- as.numeric(pchisq(q=out$statistic,df=out$parameter,lower.tail=FALSE))
pvalor
## [1] 0.4769844
q <- qchisq(p=0.05,df=out$parameter,lower.tail=FALSE);q
## [1] 14.06714
pvalor.corrigido <- as.numeric(pchisq(q=out$statistic,df=out$parameter-2,lower.tail=FALSE))
pvalor.corrigido
## [1] 0.2561615
q <- qchisq(p=0.05,df=out$parameter-2,lower.tail=FALSE);q
## [1] 11.0705Não rejeitamos \(H_{0}\), ou seja, existem evidências experimentais que \(X\) segue o modelo normal com \(\mu=14.5\) e \(\sigma^{2}=2.7\).
1.2 Teste de homogeneidade
Problema. Estamos interessados em saber se a preferência por certo tipo de filme se altera em função do estado civil. Selecionou-se 100 indivíduos em cada uma das seguintes categorias: solteiro, casado, divorciado e viúvo. Para cada categoria levantou-se a preferência dos 100 indivíduos em relação aos seguintes tipos de filmes, a saber: policial, comédia e romance. Os resultados estão apresentados no seguinte script. Se \(X\) é uma variável categórica representando o estado civil do indivíduo e \(Y\) é outra variável categórica representando sua preferência por tipo de filme deseja-se verificar se a distribuição condicional da preferência por tipo de filme dado o estado civil é homogênea com respeito às diferentes categorias de estado civil. Ou seja, deseja-se realizar o seguinte teste:
\[ \begin{cases} H_{0}: & \mbox{A distribuição de $Y$ dado $[X=x]$ é homogênea com respeito a $x$} \\ H_{1}: & \mbox{A distribuição de $Y$ dado $[X=x]$ não é homogênea com respeito a $x$.} \end{cases} \]
A função chisq.test() pode ser usada para realizar este
teste de homogeneidade.
df <- read.table(file="filmes.txt",header=TRUE,stringsAsFactors=FALSE)
df$x <- factor(df$x,levels=c("Solteiro","Casado","Divorciado","viuvo"))
df$y <- factor(df$y,levels=c("Policial","Comedia","Romance"))
str(df)
## 'data.frame': 400 obs. of 2 variables:
## $ x: Factor w/ 4 levels "Solteiro","Casado",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ y: Factor w/ 3 levels "Policial","Comedia",..: 1 3 2 2 1 1 3 3 1 1 ...
#
# Dados
#
head(df)
## x y
## 1 Solteiro Policial
## 2 Solteiro Romance
## 3 Solteiro Comedia
## 4 Solteiro Comedia
## 5 Solteiro Policial
## 6 Solteiro Policial
tail(df)
## x y
## 395 viuvo Comedia
## 396 viuvo Policial
## 397 viuvo Romance
## 398 viuvo Policial
## 399 viuvo Comedia
## 400 viuvo Comedia
tab <- table(df)
addmargins(tab)
## y
## x Policial Comedia Romance Sum
## Solteiro 45 25 30 100
## Casado 26 41 33 100
## Divorciado 36 33 31 100
## viuvo 28 38 34 100
## Sum 135 137 128 400
out <- chisq.test(x=tab);out
##
## Pearson's Chi-squared test
##
## data: tab
## X-squared = 11.256, df = 6, p-value = 0.08077
names(out)
## [1] "statistic" "parameter" "p.value" "method" "data.name" "observed"
## [7] "expected" "residuals" "stdres"
out$statistic
## X-squared
## 11.25643
out$parameter
## df
## 6
out$p.value
## [1] 0.08076717
out$method
## [1] "Pearson's Chi-squared test"
out$data.name
## [1] "tab"
out$observed
## y
## x Policial Comedia Romance
## Solteiro 45 25 30
## Casado 26 41 33
## Divorciado 36 33 31
## viuvo 28 38 34
out$expected
## y
## x Policial Comedia Romance
## Solteiro 33.75 34.25 32
## Casado 33.75 34.25 32
## Divorciado 33.75 34.25 32
## viuvo 33.75 34.25 32
q <- qchisq(p=0.05,df=out$parameter,lower.tail=FALSE);q
## [1] 12.59159Não rejeitamos \(H_{0}\), ou seja, a preferência pelo tipo de filme é a mesma, não importando o estado civil.
#
# Visualização dos dados
#
round(100*prop.table(tab,1),1)
## y
## x Policial Comedia Romance
## Solteiro 45 25 30
## Casado 26 41 33
## Divorciado 36 33 31
## viuvo 28 38 34
cores <- c("blue","red","green")
legenda <- c("P","C","R")
barplot(t(prop.table(tab,1)),col=cores,ylab="Frequência relativa condicional",ylim=c(0,1))
legend(0.30,0.98,legend=legenda,fill=cores)
dev.off()
## null device
## 1Problema. Em aparelhos dentários são usados grampos de dois tipos: um modelo em T e outro circunferencial C. Deseja-se verificar se a resistência à remoção de grampos do tipo T é a mesma dos grampos do tipo C. Foram usados 40 corpos de provas (dente-grampo), sendo 20 corpos (selecionados aleatoriamente) destinados para o tipo T e os outros 20 corpos restantes destinados para o tipo C, com 5 leituras de resistência à remoção para cada corpo de prova, totalizando 100 observações para cada tipo de grampo. O grupo de controle é aquele constituído pelos grampos do tipo T e o tratamento é aquele constituído pelos grampos do tipo C. Se \(X\) é uma variável categórica representando o tipo de grampo e \(Y\) é uma variável numérica representando a resistência à remoção do grampo (medida em kg), deseja-se verificar se a distribuição condicional da resistência à remoção do grampo dado seu tipo é homogênea com respeito aos dois tipo de grampo estudados. Ou seja, deseja-se realizar o seguinte teste:
\[ \begin{cases} H_{0}: & \mbox{A distribuição de $Y$ dado $[X=x]$ é homogênea com respeito a $x$} \\ H_{1}: & \mbox{A distribuição de $Y$ dado $[X=x]$ não é homogênea com respeito a $x$.} \end{cases} \]
df <- read.table(file="grampos.txt",header=TRUE,stringsAsFactors=FALSE)
df$x <- factor(df$x,levels=c("T","C"))
yd <- cut(df$y,breaks=seq(0.4,2.8,0.6),right=FALSE)
classes <- levels(yd);classes
## [1] "[0.4,1)" "[1,1.6)" "[1.6,2.2)" "[2.2,2.8)"
df$yd <- yd
class(df)
## [1] "data.frame"
str(df)
## 'data.frame': 200 obs. of 3 variables:
## $ x : Factor w/ 2 levels "T","C": 1 1 1 1 1 1 1 1 1 1 ...
## $ y : num 0.672 0.714 1.262 1.56 0.877 ...
## $ yd: Factor w/ 4 levels "[0.4,1)","[1,1.6)",..: 1 1 2 2 1 2 2 2 2 1 ...
#
# Dados
#
head(df)
## x y yd
## 1 T 0.6716 [0.4,1)
## 2 T 0.7144 [0.4,1)
## 3 T 1.2622 [1,1.6)
## 4 T 1.5595 [1,1.6)
## 5 T 0.8769 [0.4,1)
## 6 T 1.0647 [1,1.6)
tail(df)
## x y yd
## 195 C 1.0157 [1,1.6)
## 196 C 1.5322 [1,1.6)
## 197 C 0.8795 [0.4,1)
## 198 C 1.1393 [1,1.6)
## 199 C 1.9992 [1.6,2.2)
## 200 C 0.8686 [0.4,1)
tab <- table(df[,c(1,3)])
colnames(tab) <- c("Baixa","Moderada","Alta","Extra-alta")
dimnames(tab) <- list(Grampo=c("T","C"),Resistencia=c("Baixa","Moderada","Alta","Extra-alta"))
tab
## Resistencia
## Grampo Baixa Moderada Alta Extra-alta
## T 29 60 9 2
## C 37 44 13 6
addmargins(tab)
## Resistencia
## Grampo Baixa Moderada Alta Extra-alta Sum
## T 29 60 9 2 100
## C 37 44 13 6 100
## Sum 66 104 22 8 200
out <- chisq.test(tab)
## Warning in chisq.test(tab): Chi-squared approximation may be incorrect
out
##
## Pearson's Chi-squared test
##
## data: tab
## X-squared = 6.1585, df = 3, p-value = 0.1041
out$observed
## Resistencia
## Grampo Baixa Moderada Alta Extra-alta
## T 29 60 9 2
## C 37 44 13 6
out$expected
## Resistencia
## Grampo Baixa Moderada Alta Extra-alta
## T 33 52 11 4
## C 33 52 11 4
q <- qchisq(p=0.05,df=out$parameter,lower.tail=FALSE);q
## [1] 7.814728É importante ressaltar duas observações:
Observação: Para realização do teste, a variável numérica \(Y\) foi categorizada em 4 faixas ou categorias: resistência baixa (\(B=(0.4,1.0]\) kg), resistência moderada (\(M=(1.0,1.6]\) kg), resistência alta (\(A=(1.6,2.2]\) kg) e resistência extra-alta (\(EA=(2.2,2.8]\) kg).
Observação: O teste foi realizado e a decisão é não rejeitamos \(H_{0}\), ou seja, a distribuição da resistência à remoção do grampo é a mesma, não importando o tipo do grampo. Porém, um aviso foi apresentado informando que a aproximação qui-quadrado pode estar incorreta para este teste. Isto se deve ao fato, que algumas frequências esperadas são menores que 5 (confira acima). Existem duas opções para solucionar este problema:
- Solicitar que o p-valor do teste seja calculado via sumulação Monte Carlo.
- Agregar categorias para que todas as frequências esperadas sejam maiores do que 5.
#
# Opção 1
#
chisq.test(tab,simulate.p.value=TRUE,B=2000)
##
## Pearson's Chi-squared test with simulated p-value (based on 2000
## replicates)
##
## data: tab
## X-squared = 6.1585, df = NA, p-value = 0.1059#
# Opção 2
#
# As categorias A e EA foram agregadas e renomeadas em conjunto por A (Alta)
#
tab.modificada <- matrix(c(29,60,11,37,44,19),nrow=2,ncol=3,byrow=TRUE)
rownames(tab.modificada) <- c("T","C")
colnames(tab.modificada) <- c("Baixa","Moderada","Alta")
dimnames(tab.modificada) <- list(Grampo=c("T","C"),Resistencia=c("Baixa","Moderada","Alta"))
tab.modificada
## Resistencia
## Grampo Baixa Moderada Alta
## T 29 60 11
## C 37 44 19
out <- chisq.test(x=tab.modificada);out
##
## Pearson's Chi-squared test
##
## data: tab.modificada
## X-squared = 5.5646, df = 2, p-value = 0.0619
out$observed
## Resistencia
## Grampo Baixa Moderada Alta
## T 29 60 11
## C 37 44 19
out$expected
## Resistencia
## Grampo Baixa Moderada Alta
## T 33 52 15
## C 33 52 15
q <- qchisq(p=0.05,df=out$parameter,lower.tail=FALSE);q
## [1] 5.991465Após implementar as duas opções propostas, a decisão do teste não se alterou. Não rejeitamos \(H_{0}\), ou seja, a distribuição da resistência à remoção do grampo é a mesma, não importando o tipo do grampo.
#
# Visualização dos dados
#
round(100*prop.table(tab,1),1)
## Resistencia
## Grampo Baixa Moderada Alta Extra-alta
## T 29 60 9 2
## C 37 44 13 6
cores <- c("blue","red","green","yellow")
legenda <- c("B","M","A","EA")
barplot(t(prop.table(tab,1)),col=cores,ylab="Frequência relativa condicional",ylim=c(0,1))
legend(0.30,0.80,legend=legenda,fill=cores)
dev.off()
## null device
## 11.3 Teste de independência
Problema. Um exame com questões das disciplinas física e matemática foi aplicado a 528 estudantes do ensino médio. As notas obtidas por cada estudante em física e matemática foram registradas e classificadas nas seguintes categorias: alta, média e baixa. Os dados estão apresentados no script abaixo. Se \(X\) e \(Y\) são variáveis categóricas descrevendo as notas em física e matemática, deseja-se testar se existe dependência entre as notas dessas duas disciplinas. Ou seja, deseja-se realizar o seguinte teste:
\[ \begin{cases} H_{0}: & \mbox{$X$ e $Y$ são independentes} \\ H_{1}: & \mbox{$X$ e $Y$ não são independentes.} \end{cases} \]
A função chisq.test() pode ser usada para realizar este
teste de independência.
df <- read.table(file="notas.txt",header=TRUE,stringsAsFactors=FALSE)
df$x <- factor(df$x,levels=c("alta","media","baixa"))
df$y <- factor(df$y,levels=c("alta","media","baixa"))
str(df)
## 'data.frame': 528 obs. of 2 variables:
## $ x: Factor w/ 3 levels "alta","media",..: 3 1 1 2 3 2 2 3 1 2 ...
## $ y: Factor w/ 3 levels "alta","media",..: 3 1 2 1 1 1 2 3 2 2 ...
head(df)
## x y
## 528 baixa baixa
## 37 alta alta
## 116 alta media
## 170 media alta
## 389 baixa alta
## 178 media alta
tail(df)
## x y
## 351 media baixa
## 488 baixa baixa
## 155 media alta
## 95 alta media
## 297 media media
## 174 media alta
#
# Dados
#
tab <- table(df)
tab
## y
## x alta media baixa
## alta 56 71 12
## media 47 163 38
## baixa 14 42 85
addmargins(tab)
## y
## x alta media baixa Sum
## alta 56 71 12 139
## media 47 163 38 248
## baixa 14 42 85 141
## Sum 117 276 135 528
freq <- prop.table(tab);100*round(freq,4)
## y
## x alta media baixa
## alta 10.61 13.45 2.27
## media 8.90 30.87 7.20
## baixa 2.65 7.95 16.10
out <- chisq.test(x=tab);out
##
## Pearson's Chi-squared test
##
## data: tab
## X-squared = 145.78, df = 4, p-value < 2.2e-16
pvalor.exato <- pchisq(q=out$statistic,df=out$parameter,lower.tail=FALSE)
pvalor.exato
## X-squared
## 1.631603e-30Rejeitamos \(H_{0}\), ou seja, as notas de física e matemática não são independentes.
2 Outros testes de aderência
2.1 Teste de Kolmogorov-Smirnov
Problema. Considere os dados abaixo (veja o script), que supomos ser uma amostra de tamanho 30 de uma normal com média 10 e variância 25. Desejamos testar esta afirmação.
dados <- c(1.04,1.73,3.93,4.44,6.37,6.51,
7.61,7.64,8.18,8.48,8.57,8.65,
9.71,9.87,9.95,10.01,10.52,10.69,
11.72,12.17,12.61,12.98,13.03,13.16,
14.11,14.60,14.64,14.75,16.68,22.14)
summary(dados)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.040 7.775 9.980 10.216 13.018 22.140
var(dados)
## [1] 19.9545
sd(dados)
## [1] 4.467046par(mfrow=c(1,2))
hist(dados,breaks=seq(0,25,5),freq=FALSE,right=FALSE,
xlim=c(-5,25),ylim=c(0,0.1),col="lightblue",
xlab="x",ylab="Densidade",main="")
curve(dnorm(x,mean=10,sd=5),-5,25,col="red",add=TRUE)
boxplot(dados,xlab="x",col="lightgoldenrod")
par(mfrow=c(1,2))
plot(ecdf(dados),xlim=c(-5,25),main="")
curve(pnorm(x,mean=10,sd=5),-5,25,col="red",add=TRUE)
qqnorm(dados,type="p",pch=16,main="",xlab="Quantis teóricos",ylab="Quantis amostrais")
qqline(dados,probs=c(0.25,0.75),col="red")
dev.off()
## null device
## 1dados
## [1] 1.04 1.73 3.93 4.44 6.37 6.51 7.61 7.64 8.18 8.48 8.57 8.65
## [13] 9.71 9.87 9.95 10.01 10.52 10.69 11.72 12.17 12.61 12.98 13.03 13.16
## [25] 14.11 14.60 14.64 14.75 16.68 22.14
out <- ks.test(x=dados,y="pnorm",mean=10,sd=5)
out
##
## One-sample Kolmogorov-Smirnov test
##
## data: dados
## D = 0.11633, p-value = 0.769
## alternative hypothesis: two-sided
names(out)
## [1] "statistic" "p.value" "alternative" "method" "data.name"2.2 Teste de Shapiro-Wilk para normalidade
Problema. Considere novamente os dados abaixo (veja o script), que supomos ser uma amostra de tamanho 30 de uma normal com média 10 e variância 25. Desejamos testar esta afirmação de normalidade para os dados.
dados <- c(1.04,1.73,3.93,4.44,6.37,6.51,
7.61,7.64,8.18,8.48,8.57,8.65,
9.71,9.87,9.95,10.01,10.52,10.69,
11.72,12.17,12.61,12.98,13.03,13.16,
14.11,14.60,14.64,14.75,16.68,22.14)
dados
## [1] 1.04 1.73 3.93 4.44 6.37 6.51 7.61 7.64 8.18 8.48 8.57 8.65
## [13] 9.71 9.87 9.95 10.01 10.52 10.69 11.72 12.17 12.61 12.98 13.03 13.16
## [25] 14.11 14.60 14.64 14.75 16.68 22.14
summary(dados)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.040 7.775 9.980 10.216 13.018 22.140
var(dados)
## [1] 19.9545
sd(dados)
## [1] 4.467046
out <- shapiro.test(x=dados)
out
##
## Shapiro-Wilk normality test
##
## data: dados
## W = 0.97813, p-value = 0.7739
names(out)
## [1] "statistic" "p.value" "method" "data.name"Problema. O gerente de um restaurante está interessado na distribuição do tempo de entrega de pedidos realizados pelo aplicativo do restaurante. Para tal estudo, uma amostra aleatória simples de 32 tempos de entrega (em minutos) foram registrados (veja os dados no script abaixo). Teste a hipótese de que esses tempos formam uma amostra aleatória simples de uma distribuição normal. Posteriormente, teste a hipótese que o tempo médio de entrega de pedidos pelo aplicativo é igual a 20 minutos versus a hipótese que esse tempo médio é diferente de 20 minutos. Por fim, proporcione uma estimativa intervalar com 95% de confiança para o tempo médio de entrega de pedidos pelo aplicativo.
tempos <- c(4.24,9.94,11.31,11.84,13.89,19.97,21.18,22.15,
22.44,24.69,25.20,25.36,25.62,26.96,27.25,27.42,
29.98,32.54,32.98,34.88,35.34,35.44,35.57,36.29,
37.03,38.13,41.88,44.16,48.89,49.78,51.63,52.27)
tempos
## [1] 4.24 9.94 11.31 11.84 13.89 19.97 21.18 22.15 22.44 24.69 25.20 25.36
## [13] 25.62 26.96 27.25 27.42 29.98 32.54 32.98 34.88 35.34 35.44 35.57 36.29
## [25] 37.03 38.13 41.88 44.16 48.89 49.78 51.63 52.27
summary(tempos)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4.24 22.37 28.70 29.88 36.48 52.27
var(tempos)
## [1] 153.3435
sd(tempos)
## [1] 12.38319
tempos.d <- cut(tempos,breaks=seq(0,60,10),right=FALSE)
classes <- levels(tempos.d)
freq <- as.vector(table(tempos.d))
freq.r <- round(100*as.vector(prop.table(freq)),2)
freq.a <- cumsum(freq.r)
data.frame(classe=c(classes,"Total"),
frequencia=c(freq,sum(freq)),
freq.relativa=c(freq.r,sum(freq.r)),
freq.acumulada=c(freq.a,NaN))
## classe frequencia freq.relativa freq.acumulada
## 1 [0,10) 2 6.25 6.25
## 2 [10,20) 4 12.50 18.75
## 3 [20,30) 11 34.38 53.13
## 4 [30,40) 9 28.12 81.25
## 5 [40,50) 4 12.50 93.75
## 6 [50,60) 2 6.25 100.00
## 7 Total 32 100.00 NaNpar(mfrow=c(1,3))
hist(tempos,breaks=seq(0,60,10),freq=TRUE,right=FALSE,
col="lightblue",xlab="Tempo de entrega",ylab="Frequência",main="")
hist(tempos,breaks=seq(0,60,10),freq=FALSE,right=FALSE,
col="lightblue",xlab="Tempo de entrega",ylab="Densidade",main="")
curve(dnorm(x,mean=mean(tempos),sd=sd(tempos)),-10,70,col="red",add=TRUE)
abline(v=mean(tempos),col="red")
boxplot(tempos,ylab="Tempo de entrega",col="lightgoldenrod")
par(mfrow=c(1,2))
plot(ecdf(tempos),xlim=c(-10,70),main="")
curve(pnorm(x,mean(tempos),sd=sd(tempos)),-10,70,col="red",add=TRUE)
qqnorm(tempos,type="p",pch=16,main="",xlab="Quantis teóricos",ylab="Quantis amostrais")
qqline(tempos,probs=c(0.25,0.75),col="red")
dev.off()
## null device
## 1#
# Testes de hipóteses
#
tempos
## [1] 4.24 9.94 11.31 11.84 13.89 19.97 21.18 22.15 22.44 24.69 25.20 25.36
## [13] 25.62 26.96 27.25 27.42 29.98 32.54 32.98 34.88 35.34 35.44 35.57 36.29
## [25] 37.03 38.13 41.88 44.16 48.89 49.78 51.63 52.27
out1 <- shapiro.test(x=tempos) # teste de normalidade
out1
##
## Shapiro-Wilk normality test
##
## data: tempos
## W = 0.97573, p-value = 0.6696
out2 <- t.test(x=tempos,mu=20,alternative="two.sided",conf.level=0.95)
out2
##
## One Sample t-test
##
## data: tempos
## t = 4.5146, df = 31, p-value = 8.588e-05
## alternative hypothesis: true mean is not equal to 20
## 95 percent confidence interval:
## 25.41819 34.34743
## sample estimates:
## mean of x
## 29.88281
#
# Estimativa intervalar
#
out2$conf.int
## [1] 25.41819 34.34743
## attr(,"conf.level")
## [1] 0.95Problema. Um material radioativo emite partículas alfa ao longo do tempo. Deseja-se entender a distribuição dos tempos entre duas emissões sucessivas de partículas alfa deste material. Para tal estudo, uma amostra aleatória simples de 36 tempos entre duas emissões sucessivas foram registrados (veja os dados no script abaixo). Teste a hipótese que esses tempos aderem ao modelo normal.
temp <- c(0.02,0.04,0.05,0.06,0.07,0.08,0.09,0.10,0.12,
0.22,0.28,0.33,0.34,0.35,0.36,0.43,0.45,0.46,
0.61,0.62,0.74,0.80,0.81,0.85,0.90,0.93,1.02,
1.09,1.17,1.18,1.40,1.53,1.54,1.81,2.02,2.10)
summary(temp)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0200 0.1950 0.5350 0.6936 1.0375 2.1000
var(temp)
## [1] 0.353178
sd(temp)
## [1] 0.5942878
#
# Visualização dos dados
#
par(mfrow=c(1,2))
hist(temp,breaks=seq(0.0,2.5,0.5),freq=FALSE,
xlab="Tempo entre emissões sucessivas",
ylab="Densidade",main="")
boxplot(temp,ylab="Tempo entre emissões sucessivas",
col="lightgoldenrod")
dev.off()
## null device
## 1dados
## [1] 1.04 1.73 3.93 4.44 6.37 6.51 7.61 7.64 8.18 8.48 8.57 8.65
## [13] 9.71 9.87 9.95 10.01 10.52 10.69 11.72 12.17 12.61 12.98 13.03 13.16
## [25] 14.11 14.60 14.64 14.75 16.68 22.14
shapiro.test(x=temp)
##
## Shapiro-Wilk normality test
##
## data: temp
## W = 0.90421, p-value = 0.004458
#
# A hipótese de normalidade é rejeitada pelo
# teste de Shapiro-Wilk. Este resultado já
# era esperado a partir da visualização dos
# dados.
#
# A hipótese que os dados aderem ao modelo
# exponencial com taxa 1.45 é aceita pelo
# teste de Kolmogorov-Smirnov.
#
ks.test(x=temp,y="pexp",rate=1.45)
##
## One-sample Kolmogorov-Smirnov test
##
## data: temp
## D = 0.10318, p-value = 0.8008
## alternative hypothesis: two-sided3 Equivalências entre testes para uma mesma situação
3.1 Teste qui-quadrado e o teste-z: testes para proporção
Problema. Uma amostra de 500 habitantes de uma cidade com idade igual ou superior a 16 anos foi selecionada para participar de uma pesquisa de opinião. Cada pessoa na amostra respondeu se concordava ou não com a solução proposta pela prefeitura da cidade para resolver um problema municipal que está prejudicando seriamente a população. As opções de resposta para a pergunta são: sim (se a pessoa é favorável a solução proposta) ou não (se a pessoa não é favorável a solução proposta). Das 500 pessoas ouvidas, 300 responderam sim e 200 responderam não. Com base nos dados, deseja-se testar se a proporção da população em questão (habitantes da cidade com idade igual ou superior a 16 anos) que é favorável a solução proposta pela prefeitura é igual a 1/2 ou é diferente de 1/2.
n0 <- 200;n1 <- 300
n <- n0+n1
respostas <- c(rep("nao",n0),rep("sim",n1))
respostas <- factor(respostas,levels=c("nao","sim"))
df <- data.frame(Indivíduos=1:n,Resposta=respostas)
str(df)
## 'data.frame': 500 obs. of 2 variables:
## $ Indivíduos: int 1 2 3 4 5 6 7 8 9 10 ...
## $ Resposta : Factor w/ 2 levels "nao","sim": 1 1 1 1 1 1 1 1 1 1 ...
freq <- as.vector(table(df$Resposta))
prop <- as.vector(prop.table(table(df$Resposta)))
tabela <- data.frame(Classe=c("Não","Sim"),Frequencia=freq,Proporções=prop)
tabela
## Classe Frequencia Proporções
## 1 Não 200 0.4
## 2 Sim 300 0.6
prop.test(x=n1,n=n,p=0.5,alternative="two.sided",conf.level=0.95,correct=FALSE)
##
## 1-sample proportions test without continuity correction
##
## data: n1 out of n, null probability 0.5
## X-squared = 20, df = 1, p-value = 7.744e-06
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.5564541 0.6420210
## sample estimates:
## p
## 0.6
chisq.test(x=c(n0,n1),p=c(0.5,0.5),correct=FALSE)
##
## Chi-squared test for given probabilities
##
## data: c(n0, n1)
## X-squared = 20, df = 1, p-value = 7.744e-06
#
# Note que z^(2)=X-squared do teste qui-quadrado
#
p.hat <- n1/n;p.hat
## [1] 0.6
z <- (p.hat-0.5)/sqrt((0.5*(1-0.5))/n);z
## [1] 4.472136
z^(2)
## [1] 203.2 Teste qui-quadrado e o teste-z: testes para comparação de proporções
Problema. Para o lançamento da nova embalagem do sabonete SEBO a divisão de criação estuda duas propostas, a saber: (A) amarela com letras vermelhas ou (B) preta com letras douradas. Eles acreditam que a proposta A chama a atenção em pelo menos 5% a mais do que a proposta B. Para verificar a validade de tal informação conduziu-se o seguinte experimento: em cada um de dois supermercados semelhantes foram colocados sabonetes com cada tipo de embalagem, e a clientes selecionados aleatoriamente foi perguntado se tinham notado o sabonete e que descrevessem qual a embalagem. Os resultados estão apresentados no script abaixo.
tabela <- matrix(c(168,180,232,420),nrow=2,ncol=2,byrow=FALSE)
rownames(tabela) <- c("A","B")
colnames(tabela) <- c("Sim","Não")
dimnames(tabela) <- list(Proposta=c("A","B"),Percepeção=c("Sim","Não"))
tabela
## Percepeção
## Proposta Sim Não
## A 168 232
## B 180 420
addmargins(tabela)
## Percepeção
## Proposta Sim Não Sum
## A 168 232 400
## B 180 420 600
## Sum 348 652 1000
prop.test(x=tabela,alternative="two.sided",conf.level=0.95,correct=FALSE)
##
## 2-sample test for equality of proportions without continuity
## correction
##
## data: tabela
## X-squared = 15.232, df = 1, p-value = 9.51e-05
## alternative hypothesis: two.sided
## 95 percent confidence interval:
## 0.05930437 0.18069563
## sample estimates:
## prop 1 prop 2
## 0.42 0.30
chisq.test(x=tabela,correct=FALSE)
##
## Pearson's Chi-squared test
##
## data: tabela
## X-squared = 15.232, df = 1, p-value = 9.51e-05
p1.hat <- 168/400;p1.hat
## [1] 0.42
p2.hat <- 180/600;p2.hat
## [1] 0.3
p <- 348/1000;p
## [1] 0.348
z <- (p1.hat-p2.hat)/sqrt((p*(1-p))*((1/400)+(1/600)));z
## [1] 3.902774
z^(2)
## [1] 15.23165
#
# Note que z^(2)=X-squared do teste qui-quadrado
# Resultado análogo vale mesmo quando se aplica
# a correção de continuidade.
#
prop.test(x=tabela,alternative="two.sided",conf.level=0.95,correct=TRUE)
##
## 2-sample test for equality of proportions with continuity correction
##
## data: tabela
## X-squared = 14.707, df = 1, p-value = 0.0001256
## alternative hypothesis: two.sided
## 95 percent confidence interval:
## 0.05722104 0.18277896
## sample estimates:
## prop 1 prop 2
## 0.42 0.30
chisq.test(x=tabela,correct=TRUE)
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: tabela
## X-squared = 14.707, df = 1, p-value = 0.0001256