Noções de notação matemática para o teste t

Opa, após alguma semanas parado apanhando pra conseguir voltar o blog, surge um novo post.

Vamos ver aqui um pouquinho de como é a notação desses modelos que a gente faz, como escrever eles com aquelas letrinhas bonitinhas da matemática que normalmente confunde as pessoas.
A parte boa de tentar entender essa notação, é que ela pode ajudar nosso entendimento dos modelos, assim como melhorar nossa intuição do que deve estar acontecendo.

Bem teste-t é uma análise bem básica, que a gente ja viu milhares de vezes, como aqui e aqui, para ver diferença de médias.
Mas vemos ver de novo ele para os seguintes dados.

massa<-c(6,8,5,7,9,11)
regiao<-factor(c("Campo Grande","Campo Grande","Campo Grande","Campo Grande","Dourados","Dourados"))

Podemos fazer um gráfico rápido para facilitar o entendimento.

01

Então temos um fator, de dois níveis, que são as cidades, campo grande e dourados, e temos a massa de indivíduos de alguma espécie qualquer.

No R, fazer o teste t significa digitar

t.test(massa~regiao,var.equal=T)

E temos como resposta…

Two Sample t-test data: massa by regiao t = -3.0551, df = 4, p-value = 0.03784 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -6.6808175 -0.3191825 sample estimates: mean in group Campo Grande mean in group Dourados 6.5 10.0

Certo, mas isso na verdade a gente está fazendo um modelo linear, é algo como um caso especial, não um caso especial, simplesmente um caso de modelo linear com 1 fator de 2 níveis, falar um caso significa que queremos ver o teste-t dentro da teoria que usamos para fazer todos os modelo lineares ou modelos gerais linearizados.

Bem um modelo linear no R é feito simplesmente usando o comando lm.

lm(massa~regiao,var.equal=T)
Call: lm(formula = massa ~ regiao) Coefficients: (Intercept) regiaoDourados 6.5 3.5

Bem se a gente calcular a média das massas para as duas cidades a gente vai ver que:

> aggregate(massa,list(regiao),mean) Group.1 x 1 Campo Grande 6.5 2 Dourados 10.0

Agora esse 6.5, que está indicado como intercept da para saber o que é, é a media de campo grande, agora e esse 3.5 de dourados, e aonde estão os valores p?

Primeiro, temos que perceber que antes de tudo temos que estimar parâmetros, para construir um modelo, mas o que diabo é o modelo? Agora quando a gente tentar escrever a notação matemática, vamos conseguir talvez entender melhor, visualizar o quadro maior.

Quando a gente fala modelo, a gente ta falando em fazer uma função matemática, que representa a situação que a gente observa, qual a melhor função matemática para descrever o que a gente ta vendo?
Uma função exata não da para fazer quase sempre, um função probabilística sempre da
O porque não conseguimos fazer uma função matemática perfeita para os dados? Porque existe quase sempre um componente estocástico, um componente chamado acaso, que a gente pode esperar algo, mas não temos como saber exatamente o que. Mesmo que você tenha todos os dados dos mundo, a posição exata de cada atomo, sua carga, toda a informação, ainda sim não da para prever o tempo perfeitamente. Perfeitamente não, mas razoavelmente bem sim.

Ta ficando confuso, mas muita calma nessa hora.

Matematicamente, o que estamos fazendo aqui é uma função assim:

massa_i = \alpha + \beta \cdot regiao_i + \epsilon_i

Poxa, a resposta então é a massa, ali a gente tem a equação de reta, a+bx e um erro, esse sendo representado por esse simbolo epsilon (\epsilon_i).
Os i nas são referentes aos índices, índices do que? Dos casos, das amostras que temos, vamos olhar a massa, no R a gente tem:

> massa [1] 6 8 5 7 9 11 > massa[1] [1] 6 > massa[3] [1] 5

Então a gente tem 6 amostras, o i é o índice, quando o índice = 1, o valor de massa que observamos é 6, quando o índice=3, o valor observado é 5, isso é o i, índice de posição, o mesmo vale para a região.

Beleza, mas a notação completa do modelo é assim;

massa_i = \alpha + \beta \cdot regiao_i + \epsilon_i
\epsilon_i \sim Normal(0,\sigma^2)

Ou as vezes assim:

massa_i = Normal(\alpha + \beta \cdot regiao_i,\sigma^2)

Então a primeira linha é a equação determinística, mas nela ta somado um número que vem de uma distribuição normal com média zero, a segunda linha, no segundo caso a gente so juntou tudo isso numa linha, mas o primeiro jeito eu acho que é mais simples de entender.

Mas como a gente ve isso em relação aos nossos dados? Assim.

  6  = \alpha \cdot 1 + \beta \cdot 0 + \epsilon_1 \\  8  = \alpha \cdot 1 + \beta \cdot 0 + \epsilon_2 \\  5  = \alpha \cdot 1 + \beta \cdot 0 + \epsilon_3 \\  7  = \alpha \cdot 1 + \beta \cdot 0 + \epsilon_4 \\  9  = \alpha \cdot 1 + \beta \cdot 1 + \epsilon_5 \\  11 = \alpha \cdot 1 + \beta \cdot 1 + \epsilon_6 \\

Essencialmente esse é o meu modelo, eu to tentando ver ele em relação a cada amostra de dados que tenho.
Quando a gente estima parâmetros, o que eu quero saber é qual valor eu posso colocar para alpha(\alpha) e para beta(\beta) para solucionar esse monte de equação ai. So que lembre-se que o epsilon (\epsilon) é um valor que vem de uma distribuição normal, uma distribuição de eventos estocásticos, dados ao acaso, então ele pode ser qualquer valor, eu posso colocar qualquer número pra ele, qualquer número mesmo, ou seja, eu sempre vou conseguir solucionar todas essa equações ae em cima, para qualquer valores de alpha e beta. Por exemplo , eu posso colocar o valor de alpha como 100, beta como um milhão e epsilon como -94 e ta resolvido.

 6  = 100 \cdot 1 + 1000000 \cdot 0 + (-94)

Espero não ter errado, mas 100 vezes 1 da 100, dai um milhão vezes zero da zero, e tirando 94 da 6. So que essa solução é ruim, porque esses valores de alpha e beta fazem a gente ter que usar valores muito grandes para epsilon, mais precisamente, muito diferentes para cada equação. Imagina la em baixo quando o um milhão é multiplicado por 1, a gente vai precisar de um epsilon de mais de um milhão la, enquanto usamos 94 na primeira equação, ou seja para conseguir valores assim precisamos de um desvio muitooooo grande para o epsilon, e a pergunta é, quais valores eu ponho em alpha e beta para usar o menor desvio possível no epsilon. Pra variar o mínimo possível os valores que vou ter que colocar no epsilon para satisfazer todas essas equações ai de cima.

Certo, até aqui tudo bem, mas como isso está no R?

Bem o calculo não é feito com aquele monte equação, tem um jeito mais simples de fazer isso no R, que é usando álgebra linear(ou não, dependendo do tamanho do seu ódio para com álgebra linear). Mas a gente usa operações com matrizes.

Fica assim a conta para o nosso exemplo.

 \left(\begin{array}{c} 6 \\ 8 \\ 5 \\ 7 \\ 9 \\ 11 \end{array} \right) = \left(\begin{array}{cc} 1&0 \\ 1&0 \\ 1&0 \\ 1&0 \\ 1&1 \\ 1&1 \end{array} \right) \ast  \left(\begin{array}{c} \alpha \\ \beta \end{array} \right) +  \left(\begin{array}{c} \epsilon_1 \\ \epsilon_2 \\ \epsilon_3 \\ \epsilon_4 \\ \epsilon_5 \\ \epsilon_6 \end{array} \right)

Se a gente resolver essa multiplicação e a soma de matrizes e vetores, a gente vai resolver exatamente aquele conjunto de equaçãozinha la em cima.

Agora é so pensar, antes do igual, é nosso vetor chamado massa:

> massa [1] 6 8 5 7 9 11

Essa é a primeira matriz, antes do sinal de igual, a matriz de 1 coluna e 6 linhas(o número de linhas é o número de amostras) de resposta.

A segunda matriz, é a matriz do modelo. Pra ver ela é so ver como é guardado os dados da região.
Ele é um fator, que tem 2 componentes, os níveis e como eles estão distribuídos.

Com o comando levels a gente ve os níveis no R

> levels(regiao) [1] "Campo Grande" "Dourados" > as.integer(regiao) [1] 1 1 1 1 2 2

E temos números que representam onde eles estão, é assim porque a gente faz conta com número, para nos os nomes tem significado, mas para a álgebra das matrizes não.

Mas é so unir essas duas informações que temos nossa informação devolta

> levels(regiao)[as.integer(regiao)] [1] "Campo Grande" "Campo Grande" "Campo Grande" "Campo Grande" "Dourados" "Dourados"

Bem com essa informação, o R constrói aquela matriz que de 0 e 1 do modelo

> model.matrix(~regiao) (Intercept) regiaoDourados 1 1 0 2 1 0 3 1 0 4 1 0 5 1 1 6 1 1 attr(,"assign") [1] 0 1 attr(,"contrasts") attr(,"contrasts")$regiao [1] "contr.treatment"

E assim ja temos a nossa função, agora é só minimizar os quadrados para obter a resposta. Isso é os outros dois componente la de cima.

O comando lm ta achando aquele vetor que tem o alpha e o beta:

> lm(massa~regiao) Call: lm(formula = massa ~ regiao) Coefficients: (Intercept) regiaoDourados 6.5 3.5

Se você não se ligou ainda, essas são as médias, do que todas são dadas a partir do intercepto, ou seja a média de campo grande é 6.5, e a média de dourados é 10, que é o 6.5+3.5, aqui é quando o beta é multiplicado por 1 na matriz la em cima.

E com os menores resíduos possíveis, resido é os números que temos que atribuir a epsilon, esses aqui:

> resid(lm(massa~regiao)) 1 2 3 4 5 6 -0.5 1.5 -1.5 0.5 -1.0 1.0

Se a gente substituir esses valores la, temos a solução de todas as equações, usando o menor desvio possível.

Bem, então o nosso exemplo é assim:

 \left(\begin{array}{c} 6 \\ 8 \\ 5 \\ 7 \\ 9 \\ 11 \end{array} \right) = \left(\begin{array}{cc} 1&0 \\ 1&0 \\ 1&0 \\ 1&0 \\ 1&1 \\ 1&1 \end{array} \right) \ast  \left(\begin{array}{c} 6.5 \\ 3.5 \end{array} \right) +  \left(\begin{array}{c} -0.5 \\ 1.5 \\ -1.5 \\ 0.5 \\ -1 \\ 1 \end{array} \right)

Essa operação ai, é exatamente esse conjunto de equações.

  6  = 6.5 \cdot 1 + 3.5 \cdot 0 + -0.5 \\  8  = 6.5 \cdot 1 + 3.5 \cdot 0 + 1.5 \\  5  = 6.5 \cdot 1 + 3.5 \cdot 0 + -1.5 \\  7  = 6.5 \cdot 1 + 3.5 \cdot 0 + 0.5 \\  9  = 6.5 \cdot 1 + 3.5 \cdot 1 + -1 \\  11 = 6.5 \cdot 1 + 3.5 \cdot 1 + 1 \\

Por exemplo, o primeiro caso, 6.5 vezes 1 da 6.5, mais 3.5 vezes zero da zero, e continuamos com 6.5, tirando -0.5 da 6, que é o nosso valor observado.

E qual diabos é a função ai? Podemos escrever ela assim.

   f(x) = \left\{       \begin{array}{lr}         6.5 + \epsilon & : x = 1\\         6.5 + 3.5 + \epsilon & : x = 2       \end{array}     \right.

Onde x é nada mais que a região. E lembre-se que a gente tem a região representada assim

> regiao [1] Campo Grande Campo Grande Campo Grande Campo Grande Dourados Dourados Levels: Campo Grande Dourados

Mas o computador le assim:

> levels(regiao) [1] "Campo Grande" "Dourados" > as.integer(regiao) [1] 1 1 1 1 2 2

E dai vem o 1 e 2 da equação la em cima, x so pode assumir esses dois valores para essa equação. Agora imagine que para generalizar aquela equação para 3 níveis, como numa anova, é so adicionar uma linha, simples certo?

Agora da para visualizar como o teste-t é um caso especial de anova com apenas 2 níveis.

Bem vamos olhar mais uma vez o resultado de um modelo linear no R

> summary(lm(massa~regiao)) Call: lm(formula = massa ~ regiao) Residuals: 1 2 3 4 5 6 -0.5 1.5 -1.5 0.5 -1.0 1.0 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 6.5000 0.6614 9.827 0.000601 *** regiaoDourados 3.5000 1.1456 3.055 0.037841 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.323 on 4 degrees of freedom Multiple R-squared: 0.7, Adjusted R-squared: 0.625 F-statistic: 9.333 on 1 and 4 DF, p-value: 0.03784

Bem agora a gente sabe o que são aqueles resíduos, os coeficientes ou parâmetros estimados, e aquele itemzinho chamado “Residual standard error” é o desvio usado no Epsilon.

> sqrt(deviance(lm(massa~regiao))/df.residual(lm(massa~regiao))) [1] 1.322876

A menor quantidade de erro para satisfazer o modelo com esses parâmetros ae.

Assim se não houve-se diferença na massa entre dourados e campo grande, esperaríamos um valor de beta próximo de zero. A nossa equação ficou.

   f(x) = \left\{       \begin{array}{lr}         6.5 + \epsilon & : x = 1\\         6.5 + 3.5 + \epsilon & : x = 2       \end{array}     \right.

Deveria ter a mesma média para x=1 e x=2, o termo de diferença, que é o beta, deveria ser zero.
Como não é, é porque tem essa diferença de 3.5, lembrando que é uma diferença de mais ou menos 3.5, não exatamente, porque o epsilon pode assumir qualquer valor.

Bem, ficamos por aqui, por enquanto. Mas acho que agora talvez esses modelos façam mais sentido para todo mundo. Eu me pergunto porque eu nunca vi uma aula as pessoas mostrando esses modelos desse jeito, acho que eu entenderia melhor, mas ta bom, ta ruim mas ta bom.

Bem apesar do pau no site, o repositório recologia sempre esteve la. E é isso, até o próximo post.

massa<-c(6,8,5,7,9,11)
regiao<-factor(c("Campo Grande","Campo Grande","Campo Grande","Campo Grande","Dourados","Dourados"))
 
plot(massa~regiao,frame=F)
 
regiao
 
levels(regiao)
 
as.integer(regiao)
 
levels(regiao)[as.integer(regiao)]
 
t.test(massa~regiao,var.equal=T)
 
lm(massa~regiao)
 
resid(lm(massa~regiao))
 
model.matrix(~regiao)
 
summary(lm(massa~regiao))
 
model.matrix(~regiao-1)
 
summary(lm(massa~regiao-1))
 
aggregate(massa,list(regiao),mean)
 
sqrt(deviance(lm(massa~regiao))/df.residual(lm(massa~regiao)))

O blog está de volta!

Ae, após algumas semanas apanhando para banco de dados, SQL e tudo mais, consegui recuperar tudo sem erros de acentos.

Acredito que todos os links e posts antigos estejam funcionando, se algo ainda estiver errado por favor me avisem.

Agora estamos no 000webhost.com ja que o host antigo nos deu um pé na bunda e eu não tenho grana nem para pagar um host ultimamente. Essa semana espero voltar a postar, principalmente porque estou aprendendo bastante coisas interessantes nos cursos de bioinformática e Machine Learning no Coursera.
E nas férias, o Udacity tem um bocado de cursos novos muito interessantes, inclusive em análise de dados.
Bem até o post dessa semana.