A Inferência estatística tem como objetivo estudar generalizações sobre uma população através de evidências fornecidas por uma amostra retirada desta população.
E os dois principais tipos de inferência são:
Quantidades de interesse:
Veremos como minizar/maximizar funções, contruir intervalos de confiança e fazer testes de hipóteses!
Quando minimizar ou maximizar funções pode ser útil dentro da inferência estatística?
A função optim pode ser utilizada para encontrar os valores que minimizam uma função.
## Criando dados e fazendo um diagrama de dispersão.
dat=data.frame(x=c(1, 2, 3, 4, 5, 6),
y=c(1, 3, 5, 6, 8, 12))
with(dat, plot(x, y))
Minimizando soma dos resíduos ao quadrado:
## Função que calcula soma do quadrado dos resíduos
min.RSS <- function(data, par) {
with(data, sum((par[1] + par[2] * x - y)^2))
}
## Minimizando a soma do quadrado dos resíduos
# Note que eh necessario informar um vetor com valores iniciais
result <- optim(par = c(0, 1), min.RSS, data = dat)
result
## $par ## [1] -1.266846 2.028620 ## ## $value ## [1] 2.819048 ## ## $counts ## function gradient ## 89 NA ## ## $convergence ## [1] 0 ## ## $message ## NULL
## Comparando com a função pronta para ajustar um modelo de regressão linear lm(y ~ x, data = dat)
## ## Call: ## lm(formula = y ~ x, data = dat) ## ## Coefficients: ## (Intercept) x ## -1.267 2.029
## Gráfico de dispersão com reta ajustada plot(y ~ x, data = dat, main="Least square regression", pch=16, col=4) abline(a = result$par[1], b = result$par[2], col = "red", lwd=1.5)
fr <- function(x) { ## Rosenbrock Banana function
x1 <- x[1]
x2 <- x[2]
100 * (x2 - x1 * x1)^2 + (1 - x1)^2
}
grr <- function(x) { ## Gradient of 'fr'
x1 <- x[1]
x2 <- x[2]
c(-400 * x1 * (x2 - x1 * x1) - 2 * (1 - x1),
200 * (x2 - x1 * x1))
}
## Grafico da função que desejamos minimizar
x1 = x2 = seq(-2.5, 2.5, by=0.1)
a = expand.grid(x=x1, y=x2)
y = matrix(fr(a)[[1]], nrow=length(x1), ncol=length(x2))
persp(x1, x2, y, theta = 210, phi = 20, expand = 0.75, d=1, zlab="",
shade = 0.5, ticktype = "detailed", r=2)
# Gráfico Interativo #require(rgl) #x11() #op <- par(bg = "white") #open3d() #persp3d(x1, x2, y, expand = 0.25, zlab="", main=" ", # xlab=expression(nu[1]), ylab=expression(nu[2]), col=2) ## minimizando optim(c(-1.2, 1), fr)
## $par ## [1] 1.000260 1.000506 ## ## $value ## [1] 8.825241e-08 ## ## $counts ## function gradient ## 195 NA ## ## $convergence ## [1] 0 ## ## $message ## NULL
res <- optim(c(-1.2, 1), fr, grr, method = "BFGS") res
## $par ## [1] 1 1 ## ## $value ## [1] 9.594956e-18 ## ## $counts ## function gradient ## 110 43 ## ## $convergence ## [1] 0 ## ## $message ## NULL
contour(x1, x2, y, col="gray", nlevels = 200, xlim=c(0, 2), ylim=c(0, 2)) points(res$par[1], res$par[2], col=2)
Como podemos maximizar uma função?
Considere os dados abaixo seguindo distribuição de Poisson. Note que para maximizar bastar definir a função a ser maximizada com sinal negativo.
obs = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 17, 42, 43) freq = c(1392, 1711, 914, 468, 306, 192, 96, 56, 35, 17, 15, 6, 2, 2, 1, 1) x <- rep(obs, freq) plot(table(x), main="Count data")
## função de log-verossimilhança
fc = function(x, lamb){
sum(dpois(x, lamb, log=T))
}
# Minimizando
res = optim(par = 2, fc, x = x, control=list(fnscale=-1))
## Warning in optim(par = 2, fc, x = x, control = list(fnscale = -1)): one-dimensional optimization by Nelder-Mead is unreliable: ## use "Brent" or optimize() directly
res
## $par ## [1] 2.703516 ## ## $value ## [1] -9966.067 ## ## $counts ## function gradient ## 24 NA ## ## $convergence ## [1] 0 ## ## $message ## NULL
## Comparando com uma função pronta para encontrar a estimativas do parâmetro da Poisson library(MASS) fitdistr(x, "Poisson")
## lambda ## 2.70368239 ## (0.02277154)
Será que o resultado faz sentido?
lambda = y = seq(0.5, 6, by=0.01)
for(i in 1:length(lambda)){
y[i] = fc(x, lambda[i])
}
plot(lambda, y, type = "l", col = 3, lwd=2)
abline(v=res$par, col = 2)
Qual a utilidade?
A função uniroot pode ser utilizada para encontrar uma raiz de uma equação.
## Para uma reta
f <- function (x, a){ 3*x - a}
x = seq(0, 5, by=0.1)
plot(x, f(x, 2), type="l", lwd=2)
abline(h=0, lty=2, col=3)
res = uniroot(f, c(0, 10), tol = 0.0001, a = 2)
res
## $root ## [1] 0.6666667 ## ## $f.root ## [1] 0 ## ## $iter ## [1] 1 ## ## $init.it ## [1] NA ## ## $estim.prec ## [1] 9.333333
abline(v=res$root, lty = 3, col=2)
## Para uma parabola
f1 <- function (x, a, b, c ){ a*x^2 + b*x + c}
x = seq(-8, 5, by=0.1)
y = f1(x, a = 1, b = 3, c = -2)
plot(x, y, type="l", lwd=2)
abline(h=0, lty=2, col=3)
# Primeira raiz
res = uniroot(f1, lower = -0, upper = 1, a = 1, b = 3, c = -2)
res
## $root ## [1] 0.5615696 ## ## $f.root ## [1] 6.91674e-05 ## ## $iter ## [1] 4 ## ## $init.it ## [1] NA ## ## $estim.prec ## [1] 6.103516e-05
abline(v=res$root, lty = 3, col=2) # Segunda raiz res = uniroot(f1, lower = -8, upper = -2, a = 1, b = 3, c = -2) res
## $root ## [1] -3.561534 ## ## $f.root ## [1] -7.741396e-05 ## ## $iter ## [1] 8 ## ## $init.it ## [1] NA ## ## $estim.prec ## [1] 6.103516e-05
abline(v=res$root, lty = 3, col=2)
## Para polinômios podemos utilizar polyroot(c(-2,3,1))
## [1] 0.5615528+0i -3.5615528-0i
## Outro exemplo
plot(function(x){cos(x) - x}, from = -pi, to = pi)
raiz = uniroot(function(x) cos(x) - x, lower = -pi, upper = pi, tol = 1e-9)$root
abline(v=raiz, lty=2, col=2)
abline(h=0, lty=3, col=4)
Para encontrar múltiplas raizes podemos utilizar funções do pacote rootSolve.
require(rootSolve) # Usando uniroot.all- para raizes de uma equação fun <- function (x) cos(2*x)^3 curve(fun(x), 0, 10, main = "uniroot.all") All <- uniroot.all(fun, c(0, 10)) points(All, y = rep(0, length(All)), pch = 16, cex = 1, col=2)
## Usando o multiroot - para raizes de n equações - Exemplo 1
model <- function(x) c(F1 = x[1]^2+ x[2]^2 -1,
F2 = x[1]^2- x[2]^2 +0.5)
ss <- multiroot(f = model, start = c(1, 1))
ss
## $root ## [1] 0.5000000 0.8660254 ## ## $f.root ## F1 F2 ## 2.323138e-08 2.323308e-08 ## ## $iter ## [1] 5 ## ## $estim.precis ## [1] 2.323223e-08
## Exemplo 2
model <- function(x) c(F1 = x[1] + x[2] + x[3]^2 - 12,
F2 = x[1]^2 - x[2] + x[3] - 2,
F3 = 2 * x[1] - x[2]^2 + x[3] - 1 )
# first solution
ss <- multiroot(f = model, start = c(1, 1, 1))
ss
## $root ## [1] 1 2 3 ## ## $f.root ## F1 F2 F3 ## 3.087877e-10 4.794444e-09 -8.678146e-09 ## ## $iter ## [1] 6 ## ## $estim.precis ## [1] 4.593792e-09
# second solution; use different start values ss <- multiroot(model, c(0, 0, 0)) ss
## $root ## [1] -0.2337207 1.3531901 3.2985649 ## ## $f.root ## F1 F2 F3 ## 1.092413e-08 1.920978e-07 -4.850423e-08 ## ## $iter ## [1] 10 ## ## $estim.precis ## [1] 8.384205e-08
Suposições:
O intervalo de \(\gamma 100\%\) de confiança é dado por \[IC(\mu,\gamma)=\left[\bar{x}-z_{\gamma}\frac{\sigma}{\sqrt{n}},\ \ \bar{x}+z_{\gamma}\frac{\sigma}{\sqrt{n}}\right],\] em que \(z_{\gamma}\) é o valor da distribuição \(N(0,1)\) tal que \(P(Z<z_{\gamma})=\frac{1+\gamma}{2}\).
Exemplo: Sabe-se que o consumo mensal per capita de um determinado produto tem distribuição normal com desvio padrão de 2kg. Foi realizada uma pesquisa de mercado com 25 indivíduos dados a seguir:
9.8 8.5 7.0 10.4 8.9 7.9 5.9 7.7 8.8 7.7 7.9 5.1 9.9 7.6 8.4 3.7 6.8 7.5 8.5 5.1 7.6 5.0 6.1 6.1 10.2
Deseja-se obter um intervalo de 95% de confiança para o consumo médio mensal per capita.
Como fazer no R?
No R:
dados=c(9.8, 8.5, 7.0, 10.4, 8.9, 7.9, 5.9, 7.7, 8.8, 7.7, 7.9, 5.1, 9.9, 7.6, 8.4, 3.7, 6.8, 7.5, 8.5, 5.1, 7.6, 5.0, 6.1, 6.1, 10.2) ## Fazendo passo a passo v = 4 ## Valor assumido para variância populacional n <- length(dados) m <- mean(dados) ## Intervalo de Confiança gama = 0.95 probs = c(((1-gama)/2), (1-(1-gama)/2)) IC <- m + qnorm(probs) * sqrt(v/n) IC
## [1] 6.740014 8.307986
## Escrevendo uma funcao
ic.m <- function(x, var.pop, conf = 0.95) {
n <- length(x)
media <- mean(x)
probs = c(((1-conf)/2), (1-(1-conf)/2))
quantis <- qnorm(probs)
ic <- media + quantis * sqrt(var.pop/n)
return(ic)
}
## 95% de confiança
ic.m(dados, var.pop = 4)
## [1] 6.740014 8.307986
## 99% de confiança ic.m(dados, var.pop = 4, conf = 0.99)
## [1] 6.493668 8.554332
E necessário programar tudo isso?
## Usando pacote require(TeachingDemos) z.test(dados, mu=8, sd=2)$conf.int
## [1] 6.740014 8.307986 ## attr(,"conf.level") ## [1] 0.95
Suposições: Distribuição normal é adequada para os dados, ou a amostra é grande o suficiente, e a variância populacional \(\sigma^2\) é conhecida.
A estatística de teste é \[Z=\frac{\bar{X}-\mu_0}{\sigma/\sqrt{n}}\] Assumindo \(H_0\) verdadeira, a estatística de teste tem distribuição \(N(0,1)\), pois \(\bar{X}\sim N(\mu_0,\sigma^2/n)\).
Exemplo:
Sabe-se que o consumo mensal per capita de um determinado produto tem distribuição normal com desvio padrão de 2kg. A diretoria de uma firma que fabrica esse produto resolveu que retiraria o produto da linha de produção se a média de consumo per capita fosse menor que 8kg. Caso contrário, continuaria a fabricá-lo. Foi realizada uma pesquisa de mercado com 25 indivíduos dados a seguir:
9.8 8.5 7.0 10.4 8.9 7.9 5.9 7.7 8.8 7.7 7.9 5.1 9.9 7.6 8.4 3.7 6.8 7.5 8.5 5.1 7.6 5.0 6.1 6.1 10.2
Qual a decisão da firma a 5% de significância?
# Fazendo o teste de hipóteses z.test(dados, mu=8, alternative="less", sd=2)
## ## One Sample z-test ## ## data: dados ## z = -1.19, n = 25.0, Std. Dev. = 2.0, Std. Dev. of the sample mean = ## 0.4, p-value = 0.117 ## alternative hypothesis: true mean is less than 8 ## 95 percent confidence interval: ## -Inf 8.181941 ## sample estimates: ## mean of dados ## 7.524
Quase todas as funções possuem os seguintes resultados:
Os principais argumentos presentes nestas funções são:
Os argumentos e resultados podem variar de função para função. Dessa forma, sempre consulte o help para mais informações.
Neste caso, usamos o desvio padrão amostral e o intervalo de \(\gamma 100\%\) de confiança é dado por \[IC(\mu,\gamma)=\left[\bar{x}-t_{\gamma}\frac{s}{\sqrt{n}},\ \ \bar{x}+t_{\gamma}\frac{s}{\sqrt{n}}\right],\] onde \(t_{\gamma}\) é o valor encontrado na tabela da distribuição \(t\)-Student com \(n-1\) graus de liberdade tal que \(P(T>t_{\gamma})=\frac{1-\gamma}{2}\).
Exemplo: 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
No R:
### Com variância populacional desconhecida
## Intervalo de confiança e teste de hipoteses
## Dados
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, 3.4, 5.9, 6.3, 4.6, 5.5, 6.2)
## Obtendo o IC através da função que realiza teste de hipóteses
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
## Vendo apenas o IC res = t.test(tempo) res$conf.int
## [1] 4.278843 5.211157 ## attr(,"conf.level") ## [1] 0.95
Suposições: Os dados são normalmente distribuídos, mas não conhecemos a variância populacional \(\sigma^2\).
A estatística de teste é \[T=\frac{\bar{X}-\mu_0}{s/\sqrt{n}}.\] Assumindo \(H_0\) verdadeira, a estatística de teste tem distribuição \(t\)-Student com \(n-1\) graus de liberdade.
Exemplo: O site mercado imobiliário afirma que o preço médio de apartamentos no Barreiro é de 195 mil reais. O PROCON de Belo Horizonte diz que os preços dos imóveis naquela região tem preços maiores do que 195 mil reais. Para checar esta informação o procon colheu uma amostra dos preços de vendas de 20 imoveis vendidas no Barreiro:
198; 185; 205,2; 225,3; 206,7; 201,85; 200; 189; 192,1; 200,4
302; 195; 215,4; 235,3; 216,3; 191,85; 174,6; 199; 195,4; 202,9
Teste a afirmação do PROCON usando 10% de significância.
amostra = c(198, 185, 205.2, 225.3, 206.7, 201.85, 200, 189, 192.1, 200.4,
302, 195, 215.4, 235.3, 216.3, 191.85, 174.6, 199, 195.4, 202.9)
t.test(amostra, mu=195, alternative="greater")
## ## One Sample t-test ## ## data: amostra ## t = 1.9633, df = 19, p-value = 0.03221 ## alternative hypothesis: true mean is greater than 195 ## 95 percent confidence interval: ## 196.3792 Inf ## sample estimates: ## mean of x ## 206.565
Consideremos \(X_i\) a variável aleatória que representa a presença (ou não) de determinada característica para a \(i\)-ésima observação.
Exemplo: Numa pesquisa com 50 eleitores, o candidato José João obteve 0,34 da preferência dos eleitores. Construa, para a confiança 94%, os intervalos de confiança otimista e conservador para a proporção de votos a serem recebidos pelo candidato mencionado, supondo que a eleição fosse nesse momento.
## Preferencia de votos do José João n = 50 phat = 0.34 k = n*phat
## Podemos construir intervalos de confiança e testes de hipóteses manualmente seguindo as metodologias apresentadas nos livros de estatística básica. # Abordagem otimista SE = sqrt(phat*(1-phat)/n) E = qnorm(.975)*SE; E
## [1] 0.131303
IC1 = round(phat + c(-E, E),4) IC1
## [1] 0.2087 0.4713
# Abordagem conservadora SE = sqrt(1/(4*n)) E = qnorm(.975)*SE; E
## [1] 0.1385904
IC2 = round(phat + c(-E, E),4) IC2
## [1] 0.2014 0.4786
Também podemos utilizar a função prop.test. Ao utilizar esta função, não é efetuado o mesmo teste que aparecem nos livros de estatística básica.
#Sem correcao de continuidade: prop.test(x = k, n = n, correct = FALSE, conf.level = 0.94)
## ## 1-sample proportions test without continuity correction ## ## data: k out of n, null probability 0.5 ## X-squared = 5.12, df = 1, p-value = 0.02365 ## alternative hypothesis: true p is not equal to 0.5 ## 94 percent confidence interval: ## 0.2283482 0.4727952 ## sample estimates: ## p ## 0.34
#Com correção de continuidade: prop.test(x = k, n = n, correct = TRUE, conf.level=0.94)
## ## 1-sample proportions test with continuity correction ## ## data: k out of n, null probability 0.5 ## X-squared = 4.5, df = 1, p-value = 0.03389 ## alternative hypothesis: true p is not equal to 0.5 ## 94 percent confidence interval: ## 0.2198449 0.4829145 ## sample estimates: ## p ## 0.34
Também podemos ter a situação em que temos uma tabela com os resultados:
dados <- c(rep("sim", 120), rep("não", 80))
prop.test(x = table(dados), alternative="two.sided", conf.level = 0.94, p=0.7)
## ## 1-sample proportions test with continuity correction ## ## data: table(dados), null probability 0.7 ## X-squared = 84.292, df = 1, p-value < 2.2e-16 ## alternative hypothesis: true p is not equal to 0.7 ## 94 percent confidence interval: ## 0.3347410 0.4688671 ## sample estimates: ## p ## 0.4
Para o caso acima gostaríamos de testar se a proporção é igual a 70% ou diferente. Analisando a saída, rejeitamos a hipótese nula, ou seja, temos evidências amostrais que a proporção de votos de Joosé João é diferente de 70%.
Por último, podemos efetuar o teste binomial exato utilizando a função binom.test:
binom.test(x = k, n=n, alternative="two.sided", conf.level = 0.94, p=0.7)
## ## Exact binomial test ## ## data: k and n ## number of successes = 17, number of trials = 50, p-value = 1.791e-07 ## alternative hypothesis: true probability of success is not equal to 0.7 ## 94 percent confidence interval: ## 0.2163828 0.4820349 ## sample estimates: ## probability of success ## 0.34
# Vários outros intervalos podem ser obtidos library(binom) binom.confint(x=k, n=n, conf.level = 0.95)
## method x n mean lower upper ## 1 agresti-coull 17 50 0.3400000 0.2238941 0.4789371 ## 2 asymptotic 17 50 0.3400000 0.2086970 0.4713030 ## 3 bayes 17 50 0.3431373 0.2167061 0.4729678 ## 4 cloglog 17 50 0.3400000 0.2137049 0.4703930 ## 5 exact 17 50 0.3400000 0.2120547 0.4876525 ## 6 logit 17 50 0.3400000 0.2229732 0.4804687 ## 7 probit 17 50 0.3400000 0.2204090 0.4784223 ## 8 profile 17 50 0.3400000 0.2190717 0.4769531 ## 9 lrt 17 50 0.3400000 0.2190585 0.4769669 ## 10 prop.test 17 50 0.3400000 0.2159464 0.4885541 ## 11 wilson 17 50 0.3400000 0.2243695 0.4784617
Consideremos uma amostra aleatória \(X_1,\ldots,X_n\) de tamanho \(n\) de uma população com distribuição normal com média \(\mu\) e variância \(\sigma^2\).
Neste caso, para testar hipóteses para uma variância populacional utilizamos a estatística: \[Q=\frac{(n-1)s^2}{\sigma^2}\sim\chi_{n-1}^2.\]
Para o cálculo do intervalo de confiança, temos que:
\[IC(\sigma^2,1-\alpha)=\left(\frac{(n-1)s^2}{Q_{1-\alpha/2}},\frac{(n-1)s^2}{Q_{\alpha/2}}\right),\] em que \(Q_{\alpha/2}\) e \(Q_{1-\alpha/2}\) são tais que \(\mathbb{P}(Q < Q_{\alpha/2})=\alpha/2\) e \(\mathbb{P}(Q > Q_{\alpha/2})=\alpha/2\).
Exemplo: O peso de componentes mecânicos produzidos por uma determinada empresa é uma variável aleatória que se supõe ter distribuição normal. Pretende-se estudar a variabilidade do peso dos referidos componentes. Para isso, uma amostra de tamanho 11 foi obtida,cujos valores em grama são: 98 97 102 100 98 101 102 105 95 102 100.
Para realizar o teste e calcular o intervalo de confiança, usaremos o pacote EnviStats e a função varTest().
dados_peso <- c(98, 97, 102, 100, 98, 101, 102, 105, 95, 102, 100) varTest(dados_peso, sigma.squared = 0.5)
## ## Chi-Squared Test on Variance ## ## data: dados_peso ## Chi-Squared = 160, df = 10, p-value < 2.2e-16 ## alternative hypothesis: true variance is not equal to 0.5 ## 95 percent confidence interval: ## 3.905644 24.638334 ## sample estimates: ## variance ## 8
Acima, estamos testando se a variância é igual a 0.5. Então, pelo p-valor ao nível de significância de 5%, rejeitamos a hipótese nula, ou seja, temos evidêcia amostrais que a variância é diferente de 0.5.
Abaixo são apresentados exemplos de uso das funções para obtenção de testes e intervalos de confiança para médias, proporções e variâncias populacionais.
m1 <- c(0.9, 2.5, 9.2, 3.2, 3.7, 1.3, 1.2, 2.4, 3.6, 8.3) m2 <- c(5.3, 6.3, 5.5, 3.6, 4.1, 2.7, 2, 1.5, 5.1, 3.5) t.test(m1, m2, var.eq = TRUE, conf = 0.99) ## Assumindo variancias populacionais iguais
## ## Two Sample t-test ## ## data: m1 and m2 ## t = -0.31724, df = 18, p-value = 0.7547 ## alternative hypothesis: true difference in means is not equal to 0 ## 99 percent confidence interval: ## -3.324208 2.664208 ## sample estimates: ## mean of x mean of y ## 3.63 3.96
t.test(m1, m2, var.eq = FALSE, conf = 0.99) ## Assumindo variancias populacionais diferentes
## ## Welch Two Sample t-test ## ## data: m1 and m2 ## t = -0.31724, df = 14.028, p-value = 0.7557 ## alternative hypothesis: true difference in means is not equal to 0 ## 99 percent confidence interval: ## -3.425609 2.765609 ## sample estimates: ## mean of x mean of y ## 3.63 3.96
Observação: Para amostras pareadas é necessário utilizar o argumento paired = TRUE na função t.test
prop.test(x = c(253,196), n=c(300,300))
## ## 2-sample test for equality of proportions with continuity correction ## ## data: c(253, 196) out of c(300, 300) ## X-squared = 27.753, df = 1, p-value = 1.379e-07 ## alternative hypothesis: two.sided ## 95 percent confidence interval: ## 0.1189026 0.2610974 ## sample estimates: ## prop 1 prop 2 ## 0.8433333 0.6533333
prop.test(x = c(18,14), n=c(20,20))
## ## 2-sample test for equality of proportions with continuity correction ## ## data: c(18, 14) out of c(20, 20) ## X-squared = 1.4063, df = 1, p-value = 0.2357 ## alternative hypothesis: two.sided ## 95 percent confidence interval: ## -0.09004558 0.49004558 ## sample estimates: ## prop 1 prop 2 ## 0.9 0.7
# Amostra A ma <- c(145, 127, 136, 142, 141, 137) # Amostra B mb <- c(143, 128, 132, 138, 142, 132) # Função pronta args(var.test)
## function (x, ...) ## NULL
# Note que esta saida nao eh muito informativa. Este tipo de resultado indica que var.test()
# eh um metodo com mais de uma funcao associada. Portanto devemos pedir os argumentos da
# funcao "default".
args(getS3method("var.test", "default"))
## function (x, y, ratio = 1, alternative = c("two.sided", "less",
## "greater"), conf.level = 0.95, ...)
## NULL
var.test(ma, mb)
## ## F test to compare two variances ## ## data: ma and mb ## F = 1.0821, num df = 5, denom df = 5, p-value = 0.9331 ## alternative hypothesis: true ratio of variances is not equal to 1 ## 95 percent confidence interval: ## 0.1514131 7.7327847 ## sample estimates: ## ratio of variances ## 1.082056
Podemos utilizar o teste Qui-Quadrado para realizar:
O teste de independência Qui-Quadrado é usado para descobrir se existe uma associação entre duas variáveis em uma tabela de contingência construída à partir de dados da amostra.
Para este caso, utilizaremos o banco de dados HairEyeColor:
data(HairEyeColor) HairEyeColor
## , , Sex = Male ## ## Eye ## Hair Brown Blue Hazel Green ## Black 32 11 10 3 ## Brown 53 50 25 15 ## Red 10 10 7 7 ## Blond 3 30 5 8 ## ## , , Sex = Female ## ## Eye ## Hair Brown Blue Hazel Green ## Black 36 9 5 2 ## Brown 66 34 29 14 ## Red 16 7 7 7 ## Blond 4 64 5 8
## Juntando homens e mulheres dados = HairEyeColor[, , 1] + HairEyeColor[, , 2]
Teste com distribuição teórica da estatística de teste:
chisq.test(dados)
## ## Pearson's Chi-squared test ## ## data: dados ## X-squared = 138.29, df = 9, p-value < 2.2e-16
Teste com distribuição simulada da estatística de teste:
chisq.test(dados, sim = T)
## ## Pearson's Chi-squared test with simulated p-value (based on 2000 ## replicates) ## ## data: dados ## X-squared = 138.29, df = NA, p-value = 0.0004998
Se a variável é contínua, é necessário separá-la em faixas para executar o teste. Temos o seguindo exemplo com a distribuição exponencial:
lambda_real = 2 x = rexp(100, rate = lambda_real) quebras = c(0, 0.5, 1, 1.5, 2, +Inf) histograma = hist(x, plot=FALSE, breaks = quebras) histograma
## $breaks ## [1] 0.0 0.5 1.0 1.5 2.0 Inf ## ## $counts ## [1] 58 27 11 4 0 ## ## $density ## [1] 1.16 0.54 0.22 0.08 0.00 ## ## $mids ## [1] 0.25 0.75 1.25 1.75 Inf ## ## $xname ## [1] "x" ## ## $equidist ## [1] FALSE ## ## attr(,"class") ## [1] "histogram"
obs = histograma$counts lambda_est = 1/mean(x) esp = pexp(quebras[-1], rate = lambda_est)-pexp(quebras[-length(quebras)], lambda_est) chisq.test(obs, p=esp)
## Warning in chisq.test(obs, p = esp): Chi-squared approximation may be incorrect
## ## Chi-squared test for given probabilities ## ## data: obs ## X-squared = 3.3912, df = 4, p-value = 0.4946
chisq.test(obs,p=esp,simulate.p.value = TRUE)
## ## Chi-squared test for given probabilities with simulated p-value (based ## on 2000 replicates) ## ## data: obs ## X-squared = 3.3912, df = NA, p-value = 0.4658
Em um caso em que tenho dados do número de lentes oculares com e sem defeito, obtenho a seguinte tabela
dados<-matrix(c(36, 25, 5, 12), nrow=2,
dimnames = list(c("Lente tipo A", "Lente tipo B"),
c("Não defeituosas", "defeituosa")))
dados
## Não defeituosas defeituosa ## Lente tipo A 36 5 ## Lente tipo B 25 12
O teste de homogeneidade teste a afirmação de duas populações diferentes (com defeito, sem defeito) têm a mesma proporção de lentes defeituosas. Então:
chisq.test(dados)
## ## Pearson's Chi-squared test with Yates' continuity correction ## ## data: dados ## X-squared = 3.5612, df = 1, p-value = 0.05914
Dessa forma, vemos que o teste é marginalmente significativo. Ou seja, com 5% de significância não rejeitamos a hipótese nula, logo, temos evidências amostrais que a proporção de lenes defeituosas são iguais.
Os testes de normalidade são utilizados para verificar se a distribuição de probabilidade associada a um conjunto de dados pode ser aproximada pela distribuição normal.
No primeiro caso vamos trabalhar com dados simulados de uma distribuição normal. Para uma análise descritiva, temos que o histograma coma forma da densidade dos dados:
x.norm <- rnorm(n=80, m=10, sd=2) ## Gerando uma amostra da distribuição Normal hist(x.norm,probability = TRUE, main ="", ylab="Densidade", xlab = "Amostra")
Abaixo temos o gráfico com a distribuição acumulada dos dados, onde a linha em vermelho é a curva teórica, enquanto que os pontos são os dados simulados anteriormente.
x<-3:17 plot(x, pnorm(x, 10, 2), type="l", col="red", main="ECDF and Normal CDF", lwd=2) plot(ecdf(x.norm), add=TRUE)
No qq-plot, o ideal é que os pontos, que representam os dados simulados, estejam sob a linha teórica vermelha, ou o mais perto possivel, sem muitos desvios.
qqnorm(x.norm); qqline(x.norm, col = 2) ## Fazendo um QQ-plot
Também é possivel usar a função qqplot para outras distribuições de probabilidade. Para isso, basta trocar “rnorm” pela distribuição que queira testar.
qqplot(x.norm, rnorm(50000, m =mean(x.norm), sd=sd(x.norm))); abline(0, 1)
Vamos mostrar 2 testes de hipóteses com o objetivo de checar a normalidade dos dados:
*Shapiro-wilks:
shapiro.test(x.norm)
## ## Shapiro-Wilk normality test ## ## data: x.norm ## W = 0.99313, p-value = 0.9502
y <- rt(n=200, 2) ks.test(y, 'pnorm')
## ## One-sample Kolmogorov-Smirnov test ## ## data: y ## D = 0.11891, p-value = 0.006991 ## alternative hypothesis: two-sided