Robustez de um teste estatístico, testando o teste t

Vamos supor que eu tenho duas populações amostrais diferentes x_{1},x_{2},x_{3}, \dots ,x_{m} e y_{1},y_{2},y_{3}, \dots ,y_{n}.

A partir do teste t de student ou teste de diferença de médias eu posso perguntar se as médias são iguais.

H_{0}: \mu_x = \mu_y

Bem \mu_x é a média da distribuição normal que saíram os x_{1},x_{2},x_{3}, \dots ,x_{m}, se tivéssemos todos os infinitos valores de x, x_{1},x_{2},x_{3}, \dots ,x_{\infty}, a média seria \mu_x, mas como temos um número discreto de amostras, a gente denota que tem a média da amostra \bar{X}

Ta então a gente tem \bar{X} e \bar{Y} e seus respectivos desvios padrões  s_x e s_y

Com isso a gente pode fazer o teste t de student e verificar se as médias \bar{X} e \bar{Y} são iguais.

Realizamos a seguinte continha:

T = { { \bar{X}-\bar{Y} } \over{ S_{p}  \cdot  \sqrt{ {1 \over m }+ {1 \over n} } } }

E aquele S_{p} é o seguinte:

S_{p} = \sqrt{ (m-1) \cdot s_x^2 + (n-1) \cdot s_y^2 \over{m+n-2} }

E assim calculamos a estatística T, e sabemos que sobre a hipótese H_{0}: \mu_x = \mu_y a estatística T tem uma distribuição t com m+n-2 graus de liberdade. Claro que assumindo o seguinte:

Ambos x_{m} e y_{n} são amostras independentes e aleatorias de distribuições normais.
Os desvios padrões de x_{m} e y_{n}, os  \sigma_x e \sigma_y são iguais.

Certo então quando a gente faz um teste t, a gente ta assumindo essas duas coisas são verdades.

Ai o teste t vai consistir de uma comparação, a gente vai escolher um nível de significância \alpha que a gente quiser e vai comparar o valor que calculamos da estatística t

\vert T \vert \geq t_{n+m-2,{\alpha \over 2}}

onde t_{df,{\alpha \over 2}} é o quantile de  1 - {\alpha \over 2} de uma variável aleatória com df graus de liberdade. (df é de Degree of Freedom, grau de liberdade em inglês).

Bem se você leu até aqui parabéns, isso é meio chato né, parece muto abstrato, mas é exatamente assim que o teste t é feito. No R ou qualquer outro software. O teste t para diferenças de médias com variâncias iguais. Essa formula que mostramos ai em cima esta no wikipedia também aqui. Embaixo dela tem uma outra forma para não ter que assumir variâncias iguais, mas vamos testar esse caso aqui.

Vamos primeiro escrever na linguagem do R aquelas formulas que colocamos ali em cima.

Essa aqui:

S_{p} = \sqrt{ (m-1) \cdot s_x^2 + (n-1) \cdot s_y^2 \over{m+n-2} }

Que é o desvio padrão agrupado (pooled standard deviation) ou sei la como chama em português fica assim:

sp<-sqrt(((m-1)*sd(x)^2+(n-1)*sd(y)^2)/(m+n-2))

É so transcrevera formula para a linguagem do R, sqrt é a função que calcula raiz quadrado, ai como na formula ta tudo dentro dela, ai é uma coisa em cima e uma coisa em baixo, em cima temos m-1 onde m é o tamanho do vetor x, daria pra colocar length(x), mas salvamos o tamanho na variável m pra ficar mais curtinho e fácil de ver, ^2 eleva ao quadrado, menos pro y e o / faz a divisão. e pronto temos calculado toda aquela parte que multiplica em baixo a estatística T, agora é so fazer o mesmo esqueminha pro T inteiro.

Pegamos a formula:

T = { { \bar{X}-\bar{Y} } \over{ S_{p}  \cdot  \sqrt{ {1 \over m }+ {1 \over n} } } }

E passamos para a linguagem do R.

t<-(mean(x)-mean(y))/(sp*sqrt(1/m+1/n))

Essa foi mais fácil não? diminui as média e divide pelo sp multiplicando tudo aquilo la e pimba, temos a estatística t.

Então vamos fazer um exemplo com essa funções e ver se é isso mesmo que o R faz.

> set.seed(123456) > x<-rnorm(10,mean=50,sd=10) > y<-rnorm(10,mean=50,sd=10) > m<-length(x) > n<-length(y) > sp<-sqrt(((m-1)*sd(x)^2+(n-1)*sd(y)^2)/(m+n-2)) > t<-(mean(x)-mean(y))/(sp*sqrt(1/m+1/n)) > t [1] 1.189262

Então geramos 10 valores ao acaso, com média 50 e desvio 10, contamos o tamanho deles e salvamos em m e n e então usamos nossas formulinhas, e pimba, temos o valor t para esse caso, que é de 1.189262.

Mas para confirmar que estamos no caminho certo, o r tem um comando para fazer o teste t, o comando é t.test(), se usarmos ele e declararmos que as variações são iguais, que é o teste que estamos tentado realizar podemos conferir o valor t que ele vai calcular.

> t.test(x,y,var.equal =T,alternative = c("two.sided")) Two Sample t-test data: x and y t = 1.1893, df = 18, p-value = 0.2498 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -4.132755 14.915122 sample estimates: mean of x mean of y 57.93401 52.54283

E olha la em cima, t=1.1893, ele apenas arredondou o valor para 4 casas depois da virgula, mas é exatamente o mesmo valor, então implementamos corretamente a formula no R.

Por curiosidade, quando a gente faz um modelo linear, usando lm por exemplo:

> summary(lm(c(x,y)~rep(c("a","b"),c(m,n)))) Call: lm(formula = c(x, y) ~ rep(c("a", "b"), c(m, n))) Residuals: Min 1Q Median 3Q Max -13.682 -10.094 0.407 7.069 17.092 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 57.934 3.205 18.073 5.49e-13 *** rep(c("a", "b"), c(m, n))b -5.391 4.533 -1.189 0.25 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 10.14 on 18 degrees of freedom Multiple R-squared: 0.07285, Adjusted R-squared: 0.02134 F-statistic: 1.414 on 1 and 18 DF, p-value: 0.2498

Notem que calculamos o valor t também, é a mesma conta, so que ele saiu negativo, porque nao esta apresentado em modulo, mas o valor é o mesmo, somente arredondado para três casas decimais.

Então a gente pode colocar todos aqueles passos dentro de uma função, assim:

estatistica.t<-function(x,y) {
    m<-length(x)
    n<-length(y)
    sp<-sqrt(((m-1)*sd(x)^2+(n-1)*sd(y)^2)/(m+n-2))
    t<-(mean(x)-mean(y))/(sp*sqrt(1/m+1/n))
    return(t)
    }

E temos nossa própria função de calcular a estatística t, podemos até dar uma testada, é so entrar com dois vetores cheio de números e temos a estatística t.

> dados.x<-c(1,4,3,6,5) > dados.y<-c(5,4,7,6,10) > estatistica.t(dados.x,dados.y) [1] -1.937926

Hehe, muito divertido. Certo agora a gente pode fazer uma pequena simulação para entender como o teste funciona. Basicamente, a gente vai estabelecer m como 10 e n como 10 também, ai vamos gerar 10 valores ao acaso para um vetor x e um y, preenchendo então 10 amostras em x_{1},x_{2},x_{3}, \dots ,x_{m} e y_{1},y_{2},y_{3}, \dots ,y_{n}
vamos pegar, gerar 10 valores de uma distribuição normal com média 0 e desvio padrão de um, depois vamos calcular o valor da estatistica t com nossa função para esse exemplo e fazer essa comparação:

\vert T \vert \geq t_{n+m-2,{\alpha \over 2}}

Ai vamos registrar quantas vezes essa comparação vai ser verdadeira. Na comparação com a distribuição t, para estabelecer o valor para comparar, vamos usar um valor de alpha de 0.05, ja que é uma convenção usar esse nível de significância, ou essa chance para assumir que H_{0}: \mu_x = \mu_y esta errada,é falsa. Aqui ta meio confuso, mas o que a gente faz é basicamente isso.

#Definir um alpha > alpha<-0.05 #Tamano das amostras > m<-10 > n<-10 #Número de simulações > N<-10000 #Uma variável para contar se há diferença ou não entre as médias > n.rejeitar<-0 #Simulando N vezes > for (i in 1:N) { #Criando as duas populações com média 0, desvio 1 e m ou n amostras + x<-rnorm(m,mean=0,sd=1) + y<-rnorm(n,mean=0,sd=1) #Calculando estatistica t + t<-estatistica.t(x,y) #Se T calculado é maior que o quantile da distribuição t, conto 1 rejeição + if (abs(t)>qt(1-alpha/2,n+m-2)) { + n.rejeitar<-n.rejeitar+1 + } #E começamos denovo, até 10000 vezes + } #Depois é so ver quantas vezes isso aconteceu dado o número de simulações > nivel.sig.verdadeiro<-n.rejeitar/N > nivel.sig.verdadeiro [1] 0.0484

Certo, então rejeitamos a hipótese que as populações são iguais 0.0484 vezes, ou 04.84% das vezes, praticamente o alpha de 0.05, ou de 5% ou 1/20. Então so existe 5% de chance de a gente falar que duas populações são diferentes no caso de elas serem iguais. Quando a gente ve isso nos nossos dados, o valor T calculado é muito grande, é mais fácil procurar um motivo para a diferença nas médias que achar que elas foram diferentes ao acaso, a gente so tem 5% de chance de ta fazendo cagada.

Agora vamos fazer uma pequena modificação, ao invés de descartar o valor T que calculamos a cada teste, vamos guardar eles e fazer um gráfico para ver como é a distribuição deles para compararmos com a distribuição t, que estamos usando para tomar essa decisão.

alpha<-0.05
m<-10
n<-10
N<-10000
n.rejeitar<-0
t.caso01<-vector()
 
for (i in 1:N) {
    x<-rnorm(m,mean=0,sd=1)
    y<-rnorm(n,mean=0,sd=1)
    t.caso01[i]<-estatistica.t(x,y)
    if (abs(t.caso01[i])>qt(1-alpha/2,n+m-2)) {
        n.rejeitar<-n.rejeitar+1
        }
    }

Veja que agora eu to criando um vetor vazio t.caso01, para guardar os resultados, e toda hora que eu calculo a estatística t, eu vou guardar numa posição desse vetor.
Ai fazemos a simulação, agora guardando as valores da estatística t e colocamos isso num gráfico.

01

E olha que legal, ficam praticamente idênticos, a estatística t que calculamos e a distribuição t. Isso porque essa é exatamente a distribuição t.

Agora até aqui, a gente esta respeitando todas as premissas, mas e se a gente quebrar alguma premissa, não respeitar elas.

Vamos colocar desvios muito diferentes e ver o que acontece:

> #Duas populações normais com média zero e desvios muito diferentes (σx=1, σy=10). > nivel.sig.verdadeiro<-n.rejeitar/N > nivel.sig.verdadeiro [1] 0.065

Humm, aumentamos um pouco o número de rejeições, mas bem pouco, quer dizer que não estaremos assumindo que existe diferença onde não existe de verdade com muito maior frequência, mas se fosse o juros da minha conta do banco, eu gostaria de 1% a mais, depende então do prejuízo que esse 1% pode causar, mas a gente começar a procurar algo a mais onde não existe. Mas variâncias diferentes pode ser contornado, so ver a próxima formula no wikipedia, na pagina de teste t.
Vamos fazer o gráfico também para ver o que esse 1% a mais dessa simulação significou em termos de diferença na distribuição da estatística t.

02

Bem, ainda está bem parecido, um pouquinho desviado se você desejar ser muito detalhista, mas pra mim parece ok.
Vamos testar outro caso, duas populações vinda de uma distribuição t com 4 graus de liberdade, ela vai ser um pouco mais aberta que a distribui normal.

> #Populações T com 4 graus de liberdade e desvios iguais > nivel.sig.verdadeiro<-n.rejeitar/N > nivel.sig.verdadeiro [1] 0.046

Agora rejeitamos menos vezes? É claro que cada vez que simulamos, esse número vai mudar um pouquinho, mas são dez mil simulações, a tendencia é sempre sair algo perto disso dae. Mas estamos na verdade considerando populações iguais mais que deveríamos. Olhando o gráfico vemos…

03

O pico ficou encolhido, sei la, mas ainda ficou bem bom, e acho que menos de 1% de diferença no nível de significância é um erro aceitável, podia ser pior.

Vamos testar agora se os dados viessem de outra distribuição, uma distribuição exponencial, o que aconteceria.

> #Populações exponenciais com μx = μy = 1. > nivel.sig.verdadeiro<-n.rejeitar/N > nivel.sig.verdadeiro [1] 0.0442

Um nível de significância um pouco menor denovo, menor até que o caso anterior, olhando amostras de uma distribuição t, mas 1% até que da para levar, se isso fosse o teste para ver a eficiência de um remédio, acho que daria para usar e levar a frente as pesquisas. Se vimos uma diferença quer dizer até que o remédio pode ser um pouquinho melhor, ja que a rejeição é menor que o alpha básico de 0.05 que achamos que estamos comparando, mas vamos olhar o gráfico nesse caso e ver como fica.

04

Novamente, a gente come um pedaço do cume da montainha que é a distribuição aqui, mas ainda ta parecido, mas não é a mesma coisa, mas esta parecido.

Agora vamos testar caso de cada população vem de uma distribuição diferente.

> #Uma população normal (μx = 10, σx = 2) e uma população exponencial (μy = 10). > nivel.sig.verdadeiro<-n.rejeitar/N > nivel.sig.verdadeiro [1] 0.1042

Rapaz, agora aqui ferrou tudo, aqui esculhambou geral a estatística, o nível de significância que deveria ta pertinho do 0.05, ta em 0.1042 na verdade. Isso quer dizer que coisas que são iguais, a gente vai ta falando que são diferente, a gente vai fica procurando diferença em um monte de coisa que não tem diferença, em que as médias são iguais, se a gente olhar o gráfico desse caso…

05

Aqui não tem nem o que discutir muito, as curvas estão muito diferentes. Mas muito mesmo, nem estão no mesmo lugar. Se a gente usar o teste t quebrando as premissas desse jeito vai dar problema, um grande problema.

Mas fazendo uma revisão então:

06

Da para ver que quando as pessoas falam que um teste é robusto, eu acho que é algo no sentido que dependendo de como você quebra as premissas, ele ainda tem uma boa confiabilidade, não vamos ter conclusões muito ruins. Mas dependendo de como quebramos, como visto no ultimo exemplo, nossas conclusões vão ser péssimas, os valores que estamos lidando não tem quase nada a ver com a distribuição que estamos usando para comparar, e ficaremos tentando ver efeito onde não existe.
É claro que existe muitas outras situações ficar pensando em todas as possibilidades para simular e ver que tipo de erro vai estar acontecendo não vale muito a pena, basta ter sempre em mente que, se você esta quebrando premissas, ou não esta levando elas em conta, pelo menos olhe com cuidado as estatísticas que esta calculado e pense bem nas conclusões que esta tomando baseado nessas estatísticas.

Bem, como sempre o script fica guardado no repositório Recologia além de aqui em baixo.

set.seed(123456)
x<-rnorm(10,mean=50,sd=10)
y<-rnorm(10,mean=50,sd=10)
 
m<-length(x)
n<-length(y)
 
sp<-sqrt(((m-1)*sd(x)^2+(n-1)*sd(y)^2)/(m+n-2))
sp
 
t<-(mean(x)-mean(y))/(sp*sqrt(1/m+1/n))
t
 
t.test(x,y,var.equal =T,alternative = c("two.sided"))
summary(lm(c(x,y)~rep(c("a","b"),c(m,n))))
 
estatistica.t<-function(x,y) {
    m<-length(x)
    n<-length(y)
    sp<-sqrt(((m-1)*sd(x)^2+(n-1)*sd(y)^2)/(m+n-2))
    t<-(mean(x)-mean(y))/(sp*sqrt(1/m+1/n))
    return(t)
    }
 
dados.x<-c(1,4,3,6,5)
dados.y<-c(5,4,7,6,10)
estatistica.t(dados.x,dados.y)
 
alpha<-0.05
m<-10
n<-10
N<-10000
n.rejeitar<-0
 
for (i in 1:N) {
    x<-rnorm(m,mean=0,sd=1)
    y<-rnorm(n,mean=0,sd=1)
    t<-estatistica.t(x,y)
    if (abs(t)>qt(1-alpha/2,n+m-2)) {
        n.rejeitar<-n.rejeitar+1
        }
    }
 
nivel.sig.verdadeiro<-n.rejeitar/N
nivel.sig.verdadeiro
 
 
 
#Duas populações normais com médias zero e desvios iguais (σx = σy = 1).
alpha<-0.05
m<-10
n<-10
N<-10000
n.rejeitar<-0
t.caso01<-vector()
 
for (i in 1:N) {
    x<-rnorm(m,mean=0,sd=1)
    y<-rnorm(n,mean=0,sd=1)
    t.caso01[i]<-estatistica.t(x,y)
    if (abs(t.caso01[i])>qt(1-alpha/2,n+m-2)) {
        n.rejeitar<-n.rejeitar+1
        }
    }
 
nivel.sig.verdadeiro<-n.rejeitar/N
nivel.sig.verdadeiro
 
plot(density(t.caso01),xlim=c(min(t.caso01),max(t.caso01)),ylim=c(0,0.5),frame=F,main="",lwd=2)
lines(seq(min(t.caso01),max(t.caso01),by=0.1),dt(seq(min(t.caso01),max(t.caso01),by=0.1),df=18),lty=3,col="red",lwd=3)
legend("topright",c("Simulação","t(18)"),lwd=c(2,3),lty=c(2,3),col=c("black","red"),bty="n")
 
#Duas populações normais com média zero e desvios muito diferentes (σx=1, σy=10).
alpha<-0.05
m<-10
n<-10
N<-10000
n.rejeitar<-0
t.caso02<-vector()
 
for (i in 1:N) {
    x=rnorm(m,mean=0,sd=1)
    y=rnorm(n,mean=0,sd=10)
    t.caso02[i]<-estatistica.t(x,y)
    if (abs(t.caso02[i])>qt(1-alpha/2,n+m-2)) {
        n.rejeitar<-n.rejeitar+1
        }
    }
 
nivel.sig.verdadeiro<-n.rejeitar/N
nivel.sig.verdadeiro
 
plot(density(t.caso02),xlim=c(min(t.caso02),max(t.caso02)),ylim=c(0,0.5),frame=F,main="",lwd=2)
lines(seq(min(t.caso02),max(t.caso02),by=0.1),dt(seq(min(t.caso02),max(t.caso02),by=0.1),df=18),lty=3,col="red",lwd=3)
legend("topright",c("Simulação","t(18)"),lwd=c(2,3),lty=c(2,3),col=c("black","red"),bty="n")
 
#Populações T com 4 graus de liberdade e desvios iguais
alpha<-0.05
m<-10
n<-10
N<-10000
n.rejeitar<-0
t.caso03<-vector()
 
for (i in 1:N) {
    x=rt(m,df=4)
    y=rt(n,df=4)
    t.caso03[i]<-estatistica.t(x,y)
    if (abs(t.caso03[i])>qt(1-alpha/2,n+m-2)) {
        n.rejeitar<-n.rejeitar+1
        }
    }
 
nivel.sig.verdadeiro<-n.rejeitar/N
nivel.sig.verdadeiro
 
plot(density(t.caso03),xlim=c(min(t.caso03),max(t.caso03)),ylim=c(0,0.5),frame=F,main="",lwd=2)
lines(seq(min(t.caso03),max(t.caso03),by=0.1),dt(seq(min(t.caso03),max(t.caso03),by=0.1),df=18),lty=3,col="red",lwd=3)
legend("topright",c("Simulação","t(18)"),lwd=c(2,3),lty=c(2,3),col=c("black","red"),bty="n")
 
#Populações exponenciais com μx = μy = 1.
alpha<-0.05
m<-10
n<-10
N<-10000
n.rejeitar<-0
t.caso04<-vector()
 
for (i in 1:N) {
    x=rexp(m,rate=1)
    y=rexp(n,rate=1)
    t.caso04[i]<-estatistica.t(x,y)
    if (abs(t.caso04[i])>qt(1-alpha/2,n+m-2)) {
        n.rejeitar<-n.rejeitar+1
        }
    }
 
nivel.sig.verdadeiro<-n.rejeitar/N
nivel.sig.verdadeiro
 
plot(density(t.caso04),xlim=c(min(t.caso04),max(t.caso04)),ylim=c(0,0.5),frame=F,main="",lwd=2)
lines(seq(min(t.caso04),max(t.caso04),by=0.1),dt(seq(min(t.caso04),max(t.caso04),by=0.1),df=18),lty=3,col="red",lwd=3)
legend("topright",c("Simulação","t(18)"),lwd=c(2,3),lty=c(2,3),col=c("black","red"),bty="n")
 
#Uma população normal (μx = 10, σx = 2) e uma população exponencial (μy = 10).
alpha<-0.05
m<-10
n<-10
N<-10000
n.rejeitar<-0
t.caso05<-vector()
 
for (i in 1:N) {
    x=rnorm(m,mean=10,sd=2)
    y=rexp(n,rate=1/10)
    t.caso05[i]<-estatistica.t(x,y)
    if (abs(t.caso05[i])>qt(1-alpha/2,n+m-2)) {
        n.rejeitar<-n.rejeitar+1
        }
    }
 
nivel.sig.verdadeiro<-n.rejeitar/N
nivel.sig.verdadeiro
 
plot(density(t.caso05),xlim=c(min(t.caso05),max(t.caso05)),ylim=c(0,0.5),frame=F,main="",lwd=2)
lines(seq(min(t.caso05),max(t.caso05),by=0.1),dt(seq(min(t.caso05),max(t.caso05),by=0.1),df=18),lty=3,col="red",lwd=3)
legend("topright",c("Simulação","t(18)"),lwd=c(2,3),lty=c(2,3),col=c("black","red"),bty="n")
 
#Graficão
par(mfrow=c(3,2))
 
plot(density(t.caso01),xlim=c(min(t.caso01),max(t.caso01)),ylim=c(0,0.5),frame=F,
     main="Duas populações normais com médias zero e desvios iguais",lwd=2)
lines(seq(min(t.caso01),max(t.caso01),by=0.1),dt(seq(min(t.caso01),max(t.caso01),by=0.1),df=18),lty=3,col="red",lwd=3)
legend("topright",c("Simulação","t(18)"),lwd=c(2,3),lty=c(2,3),col=c("black","red"),bty="n")
 
plot(density(t.caso02),xlim=c(min(t.caso02),max(t.caso02)),ylim=c(0,0.5),frame=F,
     main="Duas populações normais com média zero e desvios muito diferentes",lwd=2)
lines(seq(min(t.caso02),max(t.caso02),by=0.1),dt(seq(min(t.caso02),max(t.caso02),by=0.1),df=18),lty=3,col="red",lwd=3)
legend("topright",c("Simulação","t(18)"),lwd=c(2,3),lty=c(2,3),col=c("black","red"),bty="n")
 
plot(density(t.caso03),xlim=c(min(t.caso03),max(t.caso03)),ylim=c(0,0.5),frame=F,
     main="Populações T com 4 graus de liberdade e desvios iguais",lwd=2)
lines(seq(min(t.caso03),max(t.caso03),by=0.1),dt(seq(min(t.caso03),max(t.caso03),by=0.1),df=18),lty=3,col="red",lwd=3)
legend("topright",c("Simulação","t(18)"),lwd=c(2,3),lty=c(2,3),col=c("black","red"),bty="n")
 
plot(density(t.caso04),xlim=c(min(t.caso04),max(t.caso04)),ylim=c(0,0.5),frame=F,
     main="Populações exponenciais",lwd=2)
lines(seq(min(t.caso04),max(t.caso04),by=0.1),dt(seq(min(t.caso04),max(t.caso04),by=0.1),df=18),lty=3,col="red",lwd=3)
legend("topright",c("Simulação","t(18)"),lwd=c(2,3),lty=c(2,3),col=c("black","red"),bty="n")
 
plot(density(t.caso05),xlim=c(min(t.caso05),max(t.caso05)),ylim=c(0,0.5),frame=F,
     main="Uma população normal e uma população exponencial",lwd=2)
lines(seq(min(t.caso05),max(t.caso05),by=0.1),dt(seq(min(t.caso05),max(t.caso05),by=0.1),df=18),lty=3,col="red",lwd=3)
legend("topright",c("Simulação","t(18)"),lwd=c(2,3),lty=c(2,3),col=c("black","red"),bty="n")

O que acontece se o modelo ou formula da análise que você usa esta errado? O exemplo do paradoxo de Simpson.

Essa é uma pergunta que eu sempre pensei, o que acontece quando a gente usa a análise errada? Lembrando que Anova, Regressão, geral ou não, linear ou não. Você esta informando não so os seus dados, mas também o que você acha que acontece, como o mundo funciona.
E o que acontece quando a gente tem a informação correta, mas informa a formula errada?

Um exemplo interessante é o paradoxo de Simpson. Vamos tentar inventar um exemplo aqui para entender melhor esse caso.

Vamos supor que estamos interessados em saber como a altitude (variável preditora) influencia no DAP (diâmetro a altura do peito) de uma espécie de árvore. Talvez alguém queria saber onde plantar árvores para produzir mais madeira. Vamos inventar alguns dados apenas para ilustrar a situação. Então temos uma montanha, coletamos em três lugares, base, lateral e topo da montanha, e para cada planta medimos o DPA e registramos com o GPS a sua altura dessa árvore. No final temos várias plantas com suas respectivas alturas e DAP.

01

Legal, parece que quanto maior a altitude, maior o DAP, aplicamos uma regressão linear e como resultado vemos:

> summary(modelo) Call: lm(formula = resposta ~ altitude, data = dados) Residuals: Min 1Q Median 3Q Max -2.4560 -0.7646 0.2089 0.7864 3.5954 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.4484 0.5035 0.890 0.38084 altitude 0.8883 0.3024 2.937 0.00656 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.375 on 28 degrees of freedom Multiple R-squared: 0.2355, Adjusted R-squared: 0.2082 F-statistic: 8.626 on 1 and 28 DF, p-value: 0.00656

Legal, vemos uma estimativa positiva para a altitude (0.8883), p significativo (0.00656 **) com muitas estrelinha. Ai analisamos os resíduos para ver se existe muito problemas no modelo…

02

E fora alguns outliers, parece que os resíduos tem uma distribuição normal, e não apresentam uma tendência. Talvez esse não seja o melhor exemplo, mas o problema aqui é ao simplificar a forma como olhamos os dados, não considerando onde coletamos os dados, somente aquilo que era o interesse, altura e DAP, há uma inversão no padrão que é o mais importante da nossa questão, a inclinação da reta, a relação entre altitude e DAP.
Vamos olhar novamente a figura, agora separando por cores cada população coletada, a da base, entorno e pico da montanha.

03

Agora vemos que para cada população a tendência é que quanto maior a altitude, menor o DAP.

> summary(modelo2) Call: lm(formula = resposta ~ altitude + pop, data = dados) Residuals: Min 1Q Median 3Q Max -1.6268 -0.6942 -0.1737 0.6092 1.8863 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.7965 0.4470 4.019 0.000445 *** altitude -2.2413 0.5928 -3.781 0.000826 *** poppop2 2.7382 0.6269 4.368 0.000178 *** poppop3 6.7670 1.1999 5.640 6.27e-06 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.9565 on 26 degrees of freedom Multiple R-squared: 0.6566, Adjusted R-squared: 0.617 F-statistic: 16.57 on 3 and 26 DF, p-value: 3.184e-06

O que acontece é que cada população tem um intercepto diferente, ele começa em um ponto diferente. Podemos pensar por exemplo, que com o aumento da altitude, espécies vão sumindo, e que isso deixa a espécie de nosso interesse em melhor situação, crescendo mais, mas o padrão é que quanto maior a altitude, mais o DPA diminui. Outra explicação é que altitude também altera os nutrientes disponíveis, o que pode alterar o DAP, mas ainda sim quanto maior a altitude, menor o DPA, essa relação altitude altitude vs DPA é negativa. O ponto é que na tentantiva de simplificar o mundo com menos parâmetros, invertemos a relação que observamos. A conclusão ficou totalmente diferente. Apesar que os dois casos tudo é significativo, a conclusão é totalmente diferente. Pense nos conselhos que daríamos para “onde plantar árvores” para cada caso.
Esse caso particular é conhecido como paradoxo de Simpson, que diz respeito a como um padrão se altera ao não considerarmos uma diferença entre interceptos. Eu tentei exemplificar um caso que acho mais comum a biologia, mas no artigo do Wikipédia tem excelentes exemplo.

Agora uma coisa que podemos fazer é olhar os likelihood. Nesse caso se olharmos os likelihoods, vemos que:

> logLik(modelo) 'log Lik.' -51.09479 (df=3) > logLik(modelo2) 'log Lik.' -39.0885 (df=5) > logLik(modelo)/logLik(modelo2) 'log Lik.' 1.307157 (df=3)

O “modelo” que não considera a diferença dos interceptos tem um loglikelihood maior (mais negativo) que o do “modelo2”. Podemos dividir um pelo outro e vemos que somos mais inclinados a acreditar no modelo2, se o resultado fosse 1, seria porque os likelehoods são iguais, então os modelos seriam igualmente verosimilhantes, teriam uma capacidade igual de previsão, mas desvios dos 1 dizem respeito a probabilidade de previsão melhor de um em relação ao outro. So pensar que sempre o que for valores mais da “ponta do gradiente” sempre se pareceram com outliers no “modelo” mais simples, se você não vê esses valores, cada vez vai ficar mais parecido que é inútil olhar a diferença entre as populações. Ou seja quanto menos a amostra conter todo o gradiente, pior vai ser para descobrir que caimos dentro desse paradoxo. Alias aqui se a amostra não tiver uma distribuição uniforme de altitudes, e ficar concentrado numa altitude, maior o perigo de ficar enroscado nesse tipo de paradoxo, qualidade da informação acaba sendo mais importante que quantidade.

Mas podemos ainda testar esse desvio da razão dos likelihood usando a distribuição F. Assim temos.

> anova(modelo,modelo2) Analysis of Variance Table Model 1: resposta ~ altitude Model 2: resposta ~ altitude + pop Res.Df RSS Df Sum of Sq F Pr(>F) 1 28 52.965 2 26 23.789 2 29.177 15.944 3.027e-05 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Uma diferença entre os modelos, ou seja eles são diferentes, então deveríamos considerar nesse caso aquele que tem o maior likelihood. (lembrando que o maior likelihood é o valor do modelo 2 -39.09).
Esse paradoxo eu achei legal principalmente para ter idéia de que nem sempre o modelo mais simples é o melhor. Apesar de tendermos sempre a simplificar o máximo possível os modelos do mundo seguinte a teoria da Navalha de Occan, como disse Einstein “A scientific theory should be as simple as possible, but no simpler.”

Uma teoria cientifica tem que ser o mais simples possível, mas não mais simples.

#simpson paradox
set.seed(51)
pop1<-runif(10,0,1)
pop2<-runif(10,1,2)
pop3<-runif(10,2,3)
 
resp1<-rnorm(10,2-2*pop1)
resp2<-rnorm(10,4-2*pop2)
resp3<-rnorm(10,8-2*pop3)
 
dados<-data.frame(resposta=c(resp1,resp2,resp3),altitude=c(pop1,pop2,pop3),
pop=rep(c("pop1","pop2","pop3"),each=10))
head(dados)
 
plot(resposta~altitude,data=dados,pch=19,xlab="Altitude",ylab="DAP",frame.plot=F)
abline(lm(dados$resposta~dados$altitude))
 
modelo<-lm(resposta~altitude,data=dados)
summary(modelo)
 
par(mfrow=c(2,2))
plot(modelo)
 
plot(resposta~altitude,data=dados,pch=19,col=as.numeric(dados$pop)+1,xlab="Altitude",ylab="DAP",frame.plot=F)
abline(modelo)
abline(lm(dados[dados$pop=="pop1","resposta"]~
dados[dados$pop=="pop1","altitude"]),lty=2,col=2)
abline(lm(dados[dados$pop=="pop2","resposta"]~
dados[dados$pop=="pop2","altitude"]),lty=2,col=3)
abline(lm(dados[dados$pop=="pop3","resposta"]~
dados[dados$pop=="pop3","altitude"]),lty=2,col=4)
 
modelo2<-lm(resposta~altitude+pop,data=dados)
summary(modelo2)
 
extractAIC(modelo)
extractAIC(modelo2)
 
logLik(modelo)
logLik(modelo2)
 
anova(modelo,modelo2)

Parâmetros são mais importantes que valores p

Parâmetros são mais importantes que valores p, ou pelo menos deviam ser observados com carinho.

Os valores p muitas vezes são cruéis. E vou ilustrar com uma análise que vi durante uma disciplina para a graduação em biologia, em que alunos fazem projetos rápidos como treinamento para pesquisa.

Mudanças na diversidade, riqueza ou composição de espécies são sempre um projeto batata, diferenças entre borda e centro de matas por exemplo sempre ocorrem. Essas relações entre um gradiente abiótico e mudanças na riqueza, diversidade ou composição de espécies acabam sendo comum como projetos de treinamento desse tipo.

Pois nesse caso o projeto era o seguinte, avaliar a relação da profundidade de um local de uma lagoa com a diversidade de plantas aquáticas ali. Salvo outras considerações, tudo era simples e direto. Coletar uma área, um quadrado de 1 metro quadrado, pegando todas as plantas, quantificando a abundância de cada espécie nesse quadrado e medir a profundidade da agua nesse local.

Tudo bem, coletas feitas, muitas amostras, dados foram tabulados e análises produzidas.

E la saiu uma figura muito bonita.

Feito a análise, tudo significativo e com muitas estrelinhas:

 
> summary(modelo) Call: lm(formula = diversidade ~ profundidade) Residuals: Min 1Q Median 3Q Max -0.018415 -0.006439 -0.001135 0.004694 0.017791 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.997717 0.002693 1113.041 < 2e-16 *** profundidade 0.020549 0.002307 8.907 9.72e-12 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.009179 on 48 degrees of freedom Multiple R-squared: 0.623, Adjusted R-squared: 0.6152 F-statistic: 79.34 on 1 and 48 DF, p-value: 9.717e-12

Tudo bem e todo mundo feliz. Mas havia um problema. A inclinação do modelo era um valor muito baixo. Na verdade eu não lembro os valores reais. Mas era algo como isso, de acordo com o modelo, precisaria de 5 metros a mais de profundidade pra aumentar uma espécie, ou mesmo alterações na equitabilidade que alterassem a diversidade, seria impossível com essa inclinação realmente observar algo dentro do limite físico de profundidade de lagoas que existem.

Esse gráfico refletiria melhor o que acontece, sempre tudo igual. Existe um limite físico para profundidade das lagoas e nunca veríamos uma mudança na diversidade nesse ritmo de alteração.
Mas como o gráfico anterior ficou bonito e “significativo”, pouca atenção foi dada aos parâmetros estimados. Mas o ponto é que não havia la grandes problemas na análise , e sim no fato que mesmo os parâmetros sendo significativos eles não faziam muito sentido biológico, não naquela magnitude de mudança.
Mas as pessoas acabam por odiar quem tira valores p significativos delas.

set.seed(13)
profundidade<-runif(50,0,2)
diversidade<-rnorm(50,3+0.02*profundidade,0.01)
 
#Figura 1
plot(diversidade~profundidade,xlab="Profundidade local",
ylab="Diversidade de biológica",pch=19)
abline(lm(diversidade~profundidade))
 
#Regressão linear
modelo<-lm(diversidade~profundidade)
summary(modelo)
 
#Figura 2
#preparando para desenhar o intervalo de confiança, fazendo predições
novosdados<-data.frame(profundidade=seq(0,5,by=0.01))
predição<-predict(modelo,newdata=novosdados,interval = c("prediction"))
 
plot(diversidade~profundidade,ylim=c(2.5,3.5),xlim=c(0,5),
frame.plot=F,type="n",xlab="Profundidade local",
ylab="Diversidade biológica")
 
polygon(x=c(novosdados[,1],rev(novosdados[,1])),
y=c(predição[,3],rev(predição[,2])),
col="gray90",border="gray90")
 
points(predição[,1]~novosdados[,1],type="l",lwd=2)
points(diversidade~profundidade,type="p",pch=19)
 
legend("topleft",legend=c("Modelo","Intervalo de predição"),lty=c(1,1),
lwd=c(2,20),col=c("black","gray90"))