Calculando o valor p

Algo que sempre gera discussão em estatística, é o tal do valor p. Porque quase ninguém entende direito como isso funciona na verdade. E pode me incluir nesse grupo, para falar a verdade. Nessas horas que eu penso, como estatística bayesiana é mais difícil calcular mas mais simples de pensar, para alguém que não seja exatamente um estatístico.

Existem muitas concepções erradas quanto ao valor p, uma é que ele prova, quando a gente diz que existe uma diferença significativa, o que está escrito em toda publicação quase, o que isso significa? A gente está olhando uma estatística ocorrer ao acaso, mas não vamos olhar isso agora, vamos simplesmente olhar como calcula o valor p, em termos, digamos, mais práticos.

Todo mundo já ouviu falar do teste t, de diferença de médias, distribuição t, etc. Se está por fora, de uma olhada aqui, aqui e aqui.

As vezes eu mesmo leio meus posts para relembrar as coisas, mas lendo aqueles posts acima, da para ter uma noção boa de teste de diferença de média, o teste t de Student.

Mas vamos la, calcula o valor p. Primeiro passo, é ter dados, porque sem dados nada acontece.

Então geramos uns dados, e damos uma olhada neles

01

Então temos um fator de dois níveis, chamado de A e B, veja que isso é o que está no eixo x, e temos medidas para eles.

So para visualizar de outra forma, não melhor nem pior, apenas de outra forma.

02

A gente pode olhar, que temos 30 medidas em cada tratamento, contando as bolinhas e podemos ter uma noção de média e desvio padrão olhando para as figuras, e queremos fazer a pergunta, existe diferença entre os tratamentos A e B, tudo bem que o teste t não responde exatamente essa pergunta, responde se a gente repetir esse experimento, exatamente, qual a probabilidade (p) de a diferença de médias seja diferente de zero. Talvez não exatamente isso ainda, mas eu não sou exatamente um epistemólogo mesmo.

Primeiro passo é calcular o valor t, a estatística t, que é a forma de padronizarmos a diferença de médias, colocarmos a diferença de médias num valor da distribuição t, que tem desvio zero e varia em unidades de desvio padrão.
Como fazer isso, so consultar em algum livro de estatística, ou no wikipedia, todo mundo reclama do wikipedia, mas para algumas coisas simples o wikipedia para mim é ok, principalmente coisas de ciências ou matemática, já que normalmente que faz o artigo la tem algum interesse na área, e dificilmente você vê gente vandalizando artigos de matemática por exemplo, mas vamos olhar aqui a página do teste t.

Para duas amostras, existe três formulas para o valor t.

A primeira para tamanho de amostras iguais e variâncias iguais:

t=\frac{\bar{X_1}-\bar{X_2}}{S_{X_1 X_2} \cdot \sqrt{\frac{2}{n}}}

onde S_{X_1 X_2} = \sqrt{ \frac{1}{2} (S^2_{X_1} + S^2_{X_2}) }

A segunda para tamanho de amostras iguais ou diferentes mas ainda com variâncias iguais

t=\frac{\bar{X_1}-\bar{X_2}}{S_{X_1 X_2} \cdot \sqrt{\frac{1}{n_1}+\frac{1}{n_2}}}

onde S_{X_1 X_2} = \sqrt{ \frac{(n_1-1)S^2_{X_1} + (n_2-1)S^2_{X_2}}{n_1+n_2-2}  }

E a terceira para para tamanho de amostras iguais ou diferentes porem com variâncias diferentes.

t=\frac{\bar{X_1}-\bar{X_2}}{S_{\bar{X_1}-\bar{X_2}}}

onde S_{\bar{X_1}-\bar{X_2}} = \sqrt{ \frac{S^2_1}{n_1} + \frac{S^2_2}{n_2}}

Nesse caso aqui ainda temos que calcular um valor de graus de liberdade diferente para usar para a distribuição t, que é:

df= \frac{( \frac{S^2_1}{n_1} + \frac{S^2_2}{n_2})^2}  { \frac{(\frac{S^2_1}{n_1})^2}{n-1} +    \frac{(\frac{S^2_2}{n_2})^2}{n-2}   }

O df é degree of freedom, ou grau de liberdade.

So de ver esse monte de formula as coisas já parecem complicado certo? Mas não é, tudo é bem simples. Aqui o n, é o numero de amostras, que para o A é 30 e para o B é 30 (conte as bolinhas no gráfico la em cima), ou seja, como são números iguais, podemos chamar tudo de n que é 30. Mas as segunda e a terceira formula servem para o caso onde o número de amostras diferente, o caso mais genérico de uso.

Mas para a gente, no nosso exemplo:

n=n_1=n_2=30

O S^2 é variância, S^2_1 é a variância do conjunto 1, no nosso caso o nível do tratamento A e S^2_2 a variância do conjunto 2, o nível do tratamento B. Lembrando que a raiz quadrada da variância é igual ao desvio padrão.

\bar{X_1} é média de A e \bar{X_2} é media de B

Assim, para calcular no R o valor t, basta separar os valores

1
2
3
A<-resposta[tratamento=="A"]  #Valores de A
B<-resposta[tratamento=="B"]  #Valores de B
n<-nA<-nB<-length(A)          #Tamanho das amostras

E ai é só aplicar a formula diretamente, usando a função var() para variância e mean() para média.

Primeiro para tamanho de amostras iguais e variâncias iguais:

1
2
t1<-( mean(A)-mean(B) ) /
    ( sqrt(0.5*(var(A)+var(B)))*sqrt(2/n) )

Segundo para tamanho de amostras iguais ou diferentes mas ainda com variâncias iguais

1
2
t2<-( mean(A)-mean(B) ) /
  (sqrt( ( (nA-1)*var(A)+(nB-1)*var(B) ) / (nA+nB-2) ) * sqrt(1/nA + 1/nB))

E terceira para para tamanho de amostras iguais ou diferentes porem com variâncias diferentes.

1
2
t3<-( mean(A)-mean(B) )/
    sqrt( var(A)/nA+(var(B)/nB) )

Agora olha que emocionante:

> t1 [1] -3.084094 > t2 [1] -3.084094 > t3 [1] -3.084094

É tudo o mesmo valor, mas pelo que estou falando, esses são casos especiais, a primeira formula é só uma simplificação da segunda por exemplo, olhando só o divisor, ja que o numerador é sempre igual, a média de um menos a média do outro.

Para o primeiro casos tínhamos:
\sqrt{\frac{1}{2} (S^2_{X_1} + S^2_{X_2}) } \cdot \sqrt{\frac{2}{n}}

Agora guarde ali o primeiro caso, e olhamos o segundo caso:
 \sqrt{ \frac{(n_1-1)S^2_{X_1} + (n_2-1)S^2_{X_2}}{n_1+n_2-2}  } \cdot \sqrt{\frac{1}{n_1} + \frac{1}{n_2}}

Mas aqui n=n_1=n_2 logo

 \sqrt{ \frac{(n-1)S^2_{X_1} + (n-1)S^2_{X_2}}{n+n-2}  } \cdot \sqrt{\frac{1}{n} + \frac{1}{n}}

E somando o que tem dentro da segunda raiz da

 \sqrt{ \frac{(n-1)S^2_{X_1} + (n-1)S^2_{X_2}}{n+n-2}  } \cdot \sqrt{\frac{2}{n}}

Opa, uma raiz já esta igual, agora podemos mexer na outra raiz ali.

 \sqrt{ \frac{(n-1)(S^2_{X_1} + S^2_{X_2}) }{2n-2}  } \cdot \sqrt{\frac{2}{n}}

Veja que da para isolar o (n-1)

 \sqrt{ \frac{(n-1)(S^2_{X_1} + S^2_{X_2}) }{ 2 \cdot (n-1)}  } \cdot \sqrt{\frac{2}{n}}

Ai cortamos o (n-1)

 \sqrt{ \frac{(S^2_{X_1} + S^2_{X_2}) }{ 2}  } \cdot \sqrt{\frac{2}{n}}

E se você tirar aquele 2 dali de baixo para um 1/2 multiplicando

\sqrt{\frac{1}{2} (S^2_{X_1} + S^2_{X_2}) } \cdot \sqrt{\frac{2}{n}}

Muito legal não, por isso que da o mesmo valor, da para fazer a mesma coisa ali com a outra formula, mas não vou fazer porque isso gasta muito espaço e tempo.
E veja que é exatamente essa conta que está na função t.test

> t.test(resposta~tratamento) Welch Two Sample t-test data: resposta by tratamento t = -3.0841, df = 56.559, p-value = 0.003156 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -1.1965426 -0.2543416 sample estimates: mean in group A mean in group B 0.4528962 1.1783383

Então temos o valor t que calculamos, e sabemos como é a distribuição t, vamos dar uma olhada nela aqui.

03

Agora sabendo que integrar, calcular a integral é calcular a área embaixo de uma curva, de uma olhada aqui e aqui.

Toda distribuição estatística, tem uma área embaixo dela que vale 1, 100%.

> integrate(function(x) dt(x, df=2*n-2),-Inf,Inf) 1 with absolute error < 2.4e-05

Mas suponha, que ao invés de ver tudo, os 100%, eu queira ver somente 95%
Se eu olhar do 2.5% até o 97.5%, eu consigo olhar somente 95% do meio dessa curva

> 0.975-0.025 [1] 0.95

A gente pode olhar em quais quantiles esses valores estão, quantiles são esses valores em termos de t

> qt(p=c(0.025,0.975), df=2*n-1) [1] -2.000995 2.000995

E daqui que acho que vem a confusão de depois desvios para cima e dois desvios para baixo, mas integrando a função de densidade da distribuição t nesses valores temos

> integrate(function(x) dt(x, df=2*n-2),-2.000995,2.000995) 0.94992 with absolute error < 1.2e-10

Veja que o integrate usa um método numérico, não e exato, mas isso ai é 0.95, ou 95% de probabilidade.

Dai a gente pergunta, o valor t calculado, está no meio do 95%, da área de 0.95, ou esta pra fora, é mais extremo?

04

Está pra fora, temos que p<=0.05, podemos fazer essa afirmação, a essa diferença, padronizada em valor t, é rara para esses dados, é mais fácil, parcimonioso assumir que algo acontece, algum processo gera essa diferença, e não o acaso, dai que falamos o tal do existe diferença significativa, que não é exatamente essa a afirmação, mas tudo bem. Podemos perguntar além, esse valor t, está no bolo, dentro da área de 99%?

> 0.995-0.005 [1] 0.99 > qt(p=c(0.005,0.995), df=2*n-2) [1] -2.663287 2.663287

Veja que do 0.005 ao 0.995 da uma área de 99%, calculamos quais são os quantiles que representam isso na distribuição e olhamos, se o nosso t está dentro ou fora.

05

E está fora, a ai falamos que a probabilidade é menor que 1% de esse valor ter aparecido ao acaso, e blabla p<=0.01. Com a figura fica fácil né, basicamente estamos perguntando se nosso valor calculado está dentro ou fora da área do que ocorre comumente, seja ela 95% ou 99%. Salvo definições e discussões que duram para sempre, é isso que ta rolando no teste, está dentro, está fora. Claro que no passado, essa pergunta era feita usando aquelas tabelas do fim do livro, nas tabelas as gente tinha os valores dos quantiles, e usávamos eles, comparávamos ele com nosso valor calculado para saber se observamos um t dentro ou fora dessa área. Mas era sempre somente um valor positivo, mas já já chego nesse ponto, continuando... O valor p que o teste de t nos retorna, mais do que essa resposta, é ou não é, ta dentro ou não, é o valor exato de probabilidade, probabilidade de ser menor ou mais extremo, basicamente isto. 06

Ou seja, a área que está para frente do que calculamos para t, veja que é a área fora do bolo do mais comum, a área de caras iguais ou mais incomuns que o valor t que calculamos.

E isso é a integral do menos infinito até o valor t que calculamos

> integrate(function(x) dt(x, df=2*n-2),-Inf,t1) 0.001562451 with absolute error < 2.7e-05

Veja que essa hipótese, no teste t, é perguntar se a media é menor, “less” e falar que as variâncias são iguais, senão o grau de liberdade vai cair no terceiro caso, vamos considerar o caso 1 aqui, que o grau de liberdade, que é outro papo, mas é 58, número de amostras menos 2

> t.test(resposta~tratamento,alternative = "less",var.equal = TRUE) Two Sample t-test data: resposta by tratamento t = -3.0841, df = 58, p-value = 0.001562 alternative hypothesis: true difference in means is less than 0 95 percent confidence interval: -Inf -0.3322587 sample estimates: mean in group A mean in group B 0.4528962 1.1783383

Certo, o valor da integral é o valor que vemos no p-value. Mas veja que existe uma formula fechada para calcular isso, com calculo, e ela ta implementada na função pt

> pt(t1,df=2*n-2) [1] 0.001562452

O mesmo valor, no caso o test.t vai usar o pt, e não usa o integrate que aproxima a integral, usa o pt que faz um calculo exato e é mais eficiente computacionalmente, já que faz uma conta e pronto..

Mas não acaba ae, veja que a gente poderia ter a diferença da média de B para A, o R organiza os fatores em ordem alfabetica, mas as vezes não é isso que queremos, queremos por exemplo que tudo seja comparado a partir de um controle, mas as vezes, se você usa adudo como nome do seu nivel do fator, adubo vem antes de controle na ordem lexografica, logo você quer mudar a ordem, vamos mudar aqui o B pelo A, para calcular o t

> t1<-( mean(B)-mean(A) ) / + ( sqrt(0.5*(var(B)+var(A)))*sqrt(2/n) ) > round(t1,digits=4) [1] 3.0841

E a única coisa que mudou foi o sinal. E queremos essa área agora

07

Agora se você usar a função t.test, como antes, mudando quem é o nível base, o intercepto

> t.test(resposta~tratamento,alternative = "less",var.equal = TRUE) Two Sample t-test data: resposta by tratamento t = 3.0841, df = 58, p-value = 0.9984 alternative hypothesis: true difference in means is less than 0 95 percent confidence interval: -Inf 1.118625 sample estimates: mean in group B mean in group A 1.1783383 0.4528962

Veja que para a mesma pergunta, não veremos mais diferença, porque estamos integrando pro lado errado. Do ponto de vista de quem queremos como “controle”.

Agora a gente não quer mais integrar do menos infinito até o valor de t calculado

> integrate(function(x) dt(x, df=2*n-2),-Inf,t1) 0.9984375 with absolute error < 3.1e-05

Mas do valor calculado até o valor infinito

> integrate(function(x) dt(x, df=2*n-2),t1,Inf) 0.001562451 with absolute error < 2.7e-05

Veja que aqui temos o mesmo valor de p de antes, ou seja, temos 2 hipóteses, da diferença de médias entre B e A ser negativa (“less” como argumento do t.test), B menor que A, que é so olhar o gráfico para ver que não é isso, ou B é maior (“greater” como argumento do t.test), que é o que observamos no gráfico, a hipótese de que B é menor, não é verdade, olhe o gráfico, a outra hipotese que diz que B é maior, que é a segunda integral, com mesma probabilidade do caso anterior é a que nos interessaria aqui.

Mas aqui começa a ficar confuso certo? Dai alguém pensou, meu, so tenta ver se tem diferença, abstrai essa confusão, não interessa quem é maior, quem é menor, para que lado está o t, menor ou maior, primeiro eu respondo se existe diferença, depois eu olho no gráfico e vejo quem é menor e tal. Eu tento ver essas probabilidades aqui

09

Além da propriedade de sempre integrar para 1, a distribuição t, assim como a normal, ela é simétrica, ou seja o que você ve de um lado dela, é igual do outro, como um espelho, sendo que o espelho está no 0.

10

Então a gente poderia integrar usando integrate

> integrate(function(x) dt(x, df=2*n-2),-Inf,-t1) 0.001562451 with absolute error < 2.7e-05 > integrate(function(x) dt(x, df=2*n-2),t1,Inf) 0.001562451 with absolute error < 2.7e-05

E somar as duas probabilidades para saber se existe diferença (ou a chance de se repetir o experimento blablabla…).

> integrate(function(x) dt(x, df=2*n-2),-Inf,-t1)$value+integrate(function(x) dt(x, df=2*n-2),t1,Inf)$value [1] 0.003124903

Agora para usar o pt, a gente precisa daquela propriedade da simetria, porque ele calcula a integral do menos infinito até o valor de t fornecido, assim, se o fornecido t é positivo, a gente calcula a área do lado errado que queremos

> pt(t1,df=2*n-2) [1] 0.9984375

Mas como sabemos que queremos o outro lado, e toda a área vale 1, então se pegarmos 1 menos esse valor

> 1-pt(t1,df=2*n-2) [1] 0.001562452

Temos a probabilidade que queremos, ou seja 1 menos o que calculamos é igual ao menos o valor retornado por pt de -t

> pt(-t1,df=2*n-2) [1] 0.001562452 > 1-pt(t1,df=2*n-2) [1] 0.001562452

Veja que é so comparar os dois e pimba…

> pt(-t1,df=2*n-2)==(1-pt(t1,df=2*n-2)) [1] FALSE

Ohh meu deus, o R ta falando que dois número iguais são diferentes, bem vindo ao inferno da precisão no computador. Os números tem um tamanho máximo de zeros depois da vírgula, são ponto flutuantes, e segundo o faq do R, direto o povo manda e-mail falando que encontrou um erro no R, mas isso é uma questão da memória apenas, de como a informação é guardada, e a gente vê muito isso em estatística, porque os números comumente ficam muito, muito pequenos, tanto que a verosimilhança é calculado como log, log(verosimilhança), loglikelihood, pra tentar escapar um pouco disso, além de facilitar a matemática. Mas existe uma função no R para fazer comparações desse tipo, e não cair nesse problema.

> all.equal(pt(-t1,df=2*n-2),1-pt(t1,df=2*n-2)) [1] TRUE

E sim, são a mesma coisa, ufa, quase achei que tudo que eu sabia estava errado a primeira vez que vi isso, mas muita calma nessa hora. Voltando ao ponto, que pensando desse jeito, normalmente nos exemplos, probabilidades seram calculadas assim nos programas, e exemplos.

> 2*pt(-abs(t1),df=2*n-2) [1] 0.003124903

Duas vezes menos o valor absoluto de t, valor absoluto para garantir que ele ta positivo, e menos porque queremos o cantinho da esquerda, porque é assim que o pt calcula, e vai estar assim para muitas outras distribuições simétricas, veja que toda distribuição que tem no R, tem uma função p_nome_da_distribuição, para calcular probabilidades, valores p, senão, não podemos fazer afirmações.

E o teste t é por padrão desse jeito, bi-caudal, two-sided, não interessa quem é maior que quem, so qual a probabilidade de existir diferença, mais facil de não fazer cagada, que pode acontecer dependendo de escolher entre qual cauda olhar, less ou greater

> t.test(resposta~tratamento,var.equal=TRUE) Two Sample t-test data: resposta by tratamento t = 3.0841, df = 58, p-value = 0.003125 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: 0.2545972 1.1962870 sample estimates: mean in group B mean in group A 1.1783383 0.4528962

Mas além disso, ele usa a aproximação de Welch para o grau de liberdade, ou seja, ele usa aquela formula 3 para os graus de liberdade, e não simplismente 58 graus de liberdade, veja o grau de liberdade, o df aqui em baixo é até um numero quebrado, que vem daquela fração la da fórmula 3.

> t.test(resposta~tratamento) Welch Two Sample t-test data: resposta by tratamento t = 3.0841, df = 56.559, p-value = 0.003156 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: 0.2543416 1.1965426 sample estimates: mean in group B mean in group A 1.1783383 0.4528962

Bem é isso ai, o script vai estar la no repositório recologia, eu escrevi esse post, porque esses dias fui conferir um resultado de teste t num trabalho, e tive que parar para pensar no valor p, como calcular ele. Se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail, hoje não tem referência, e até mais.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
#http://recologia.com.br/2015/05/calculando-o-valor-p/
##################################
## Gerando dados
##################################
set.seed(123)
tratamento<-factor(rep(c("A","B"),each=30))
resposta<-rnorm(60,as.numeric(tratamento)*0.5,1)
 
#jpeg("01.jpg")
boxplot(resposta~tratamento,xlab="tratamento")
#dev.off()
 
#jpeg("02.jpg")
stripchart(resposta~tratamento,vertical=T,method="jitter",pch=19,
           at=c(1.2,1.8),xlab="tratamento")
#dev.off()
 
##################################
## Calculando estatística t
##################################
A<-resposta[tratamento=="A"]  #Valores de A
B<-resposta[tratamento=="B"]  #Valores de B
n<-nA<-nB<-length(A)          #Tamanho das amostras
 
#Equal sample sizes, equal variance
t1<-( mean(A)-mean(B) ) /
    ( sqrt(0.5*(var(A)+var(B)))*sqrt(2/n) )
round(t1,digits=4)
 
#Equal or unequal sample sizes, equal variance
t2<-( mean(A)-mean(B) ) /
  (sqrt( ( (nA-1)*var(A)+(nB-1)*var(B) ) / (nA+nB-2) ) * sqrt(1/nA + 1/nB))
round(t2,digits=4)
 
#Equal or unequal sample sizes, unequal variances
t3<-( mean(A)-mean(B) )/
    sqrt( var(A)/nA+(var(B)/nB) )
round(t3,digits=4)
 
#Distribuição t
#jpeg("03.jpg")
plot(dt(x=seq(-4,4,by=0.01), df=2*n-2)~seq(-4,4,by=0.01),
     type="l",ylab="Probabilidade",xlab="Valor t")
 
arrows(x0=t1,x1=t1 ,y0=0.05, y1 = 0,lwd=2,length = 0.10)
text(x=t1,y=0.06,"Valor t calculado")
#dev.off()
 
integrate(function(x) dt(x, df=2*n-2),-Inf,Inf)
 
0.975-0.025
 
qt(p=c(0.025,0.975), df=2*n-1)
integrate(function(x) dt(x, df=2*n-2),-2.000995,2.000995)
 
#95%
#jpeg("04.jpg")
plot(dt(x=seq(-4,4,by=0.01), df=2*n-2)~seq(-4,4,by=0.01),
     type="l",ylab="Probabilidade",xlab="Valor t")
 
polygon(x= c(seq(-2.000995,2.000995,length=1000) ),
        y= c(0, dt(x=seq(-2.000995,2.000995,length=998), df=2*n-1),0 ),col="red"  )
 
arrows(x0=t1,x1=t1 ,y0=0.05, y1 = 0,lwd=2,length = 0.10)
text(x=t1,y=0.06,"Valor t calculado")
#dev.off()
 
#99%
0.995-0.005
qt(p=c(0.005,0.995), df=2*n-2)
 
#jpeg("05.jpg")
plot(dt(x=seq(-4,4,by=0.01), df=2*n-2)~seq(-4,4,by=0.01),
     type="l",ylab="Probabilidade",xlab="Valor t")
 
polygon(x= c(seq(-2.661759,2.661759,length=1000) ),
        y= c(0, dt(x=seq(-2.661759,2.661759,length=998), df=2*n-2),0 ),col="red"  )
arrows(x0=t1,x1=t1 ,y0=0.05, y1 = 0,lwd=2,length = 0.10)
text(x=t1,y=0.06,"Valor t calculado")
#dev.off()
 
#jpeg("06.jpg")
plot(dt(x=seq(-4,4,by=0.01), df=2*n-2)~seq(-4,4,by=0.01),
     type="l",ylab="Probabilidade",xlab="Valor t")
polygon(x= c(seq(t1,-4,length=1000) ),
        y= c(0, dt(x=seq(t1,-4,length=998), df=2*n-2),0 ),col="red")
arrows(x0=t1,x1=t1 ,y0=0.05, y1 = 0,lwd=2,length = 0.10)
text(x=t1,y=0.06,"Valor t calculado")
#dev.off()
 
#valor p
integrate(function(x) dt(x, df=2*n-2),-Inf,t1)
t.test(resposta~tratamento,alternative = "less",var.equal = TRUE)
 
pt(t1,df=2*n-2)
 
#invertendo a ordem da pergunta
t1<-( mean(B)-mean(A) ) /
    ( sqrt(0.5*(var(B)+var(A)))*sqrt(2/n) )
round(t1,digits=4)
 
#jpeg("07.jpg")
plot(dt(x=seq(-4,4,by=0.01), df=2*n-2)~seq(-4,4,by=0.01),
     type="l",ylab="Probabilidade",xlab="Valor t")
polygon(x= c(seq(t1,4,length=1000) ),
        y= c(0, dt(x=seq(t1,4,length=998), df=2*n-2),0 ),col="red")
arrows(x0=t1,x1=t1 ,y0=0.05, y1 = 0,lwd=2,length = 0.10)
text(x=t1,y=0.06,"Valor t calculado")
#dev.off()
 
#mudando quem é o intercepto
tratamento <- relevel(tratamento, "B")
t.test(resposta~tratamento,alternative = "less",var.equal = TRUE)
 
#jpeg("08.jpg")
plot(dt(x=seq(-4,4,by=0.01), df=2*n-2)~seq(-4,4,by=0.01),
     type="l",ylab="Probabilidade",xlab="Valor t")
polygon(x= c(seq(-4,t1,length=1000) ),
        y= c(0, dt(x=seq(-4,t1,length=998), df=2*n-2),0 ),col="red")
arrows(x0=t1,x1=t1 ,y0=0.05, y1 = 0,lwd=2,length = 0.10)
text(x=t1,y=0.06,"Valor t calculado")
#dev.off()
 
integrate(function(x) dt(x, df=2*n-2),-Inf,t1)
integrate(function(x) dt(x, df=2*n-2),t1,Inf)
 
#
#jpeg("09.jpg")
plot(dt(x=seq(-4,4,by=0.01), df=2*n-2)~seq(-4,4,by=0.01),
     type="l",ylab="Probabilidade",xlab="Valor t")
polygon(x= c(seq(t1,4,length=1000) ),
        y= c(0, dt(x=seq(t1,4,length=998), df=2*n-2),0 ),col="red")
polygon(x= c(seq(-4,-t1,length=1000) ),
        y= c(0, dt(x=seq(-4,-t1,length=998), df=2*n-2),0 ),col="red")
arrows(x0=t1,x1=t1 ,y0=0.05, y1 = 0,lwd=2,length = 0.10)
arrows(x0=-t1,x1=-t1 ,y0=0.05, y1 = 0,lwd=2,length = 0.10)
text(x=t1,y=0.06,"Valor t calculado")
text(x=-t1,y=0.06,"Menos Valor t calculado")
#dev.off()
 
#jpeg("10.jpg")
plot(dt(x=seq(-4,4,by=0.01), df=2*n-2)~seq(-4,4,by=0.01),
     type="l",ylab="Probabilidade",xlab="Valor t")
polygon(x= c(seq(t1,4,length=1000) ),
        y= c(0, dt(x=seq(t1,4,length=998), df=2*n-2),0 ),col="red")
polygon(x= c(seq(-4,-t1,length=1000) ),
        y= c(0, dt(x=seq(-4,-t1,length=998), df=2*n-2),0 ),col="red")
abline(v=0,lwd=3,lty=3)
points(c(0,t1),c(0.05,0.05),type="l",lwd=3,lty=3)
points(c(0,-t1),c(0.05,0.05),type="l",lwd=3,lty=3)
 
arrows(x0=t1,x1=t1 ,y0=0.05, y1 = 0,lwd=3,lty=3,length = 0.10)
arrows(x0=-t1,x1=-t1 ,y0=0.05, y1 = 0,lwd=3,lty=3,length = 0.10)
#dev.off()
 
integrate(function(x) dt(x, df=2*n-2),-Inf,-t1)
integrate(function(x) dt(x, df=2*n-2),t1,Inf)
integrate(function(x) dt(x, df=2*n-2),-Inf,-t1)$value+integrate(function(x) dt(x, df=2*n-2),t1,Inf)$value
 
#
pt(-t1,df=2*n-2)
1-pt(t1,df=2*n-2)
 
#
pt(-t1,df=2*n-2)==(1-pt(t1,df=2*n-2))
all.equal(pt(-t1,df=2*n-2),1-pt(t1,df=2*n-2))
 
#
2*pt(-abs(t1),df=2*n-2)
t.test(resposta~tratamento,var.equal=TRUE)
t.test(resposta~tratamento)