[R-br] [OFF-Topic] Distribuição Normal e teste de t

[Tropidurus Torquatus]

Caros R-users,

Coloquei como off-topic pois acredito que esta seja mais uma dúvida conceitual do que de R ppd.

Dado que tenho dois conjuntos de dados:

x <- c ( 13.772 , 10.119 , 11.224 , 21.999 , 17.181 , 27.477 , 21.693 , 26.706 , 14.576 , 19.689 , 12.494 , 19.511 , 16.785 , 16.551 , 20.780 ,  0.000 , 25.867 , 28.223 , 26.059 , 21.668 , 18.209 , 20.186 , 28.320 , 26.946 , 16.218 , 16.814 , 16.859 ,  0.000 ,  0.000 ,  3.986 ,  1.625 , 17.286 , 28.149 , 19.404 , 21.669 , 22.058 , 14.947 , 14.512 , 21.177 , 29.997 , 20.958 , 15.123 , 28.932 , 23.807 , 30.589 , 10.208 , 20.076 , 18.897 , 14.666 , 13.955)

y <-  c ( 13.902 , 10.102 , 11.301 , 22.081 , 17.234 , 27.540 , 21.741 , 26.677 , 14.621 , 19.766 , 12.547 , 19.622 , 16.857 , 16.628 , 20.787 ,  0.000 , 25.943 , 28.313 , 26.166 , 21.724 , 18.255 , 20.307 , 28.387 , 27.079 , 16.202 , 16.857 , 16.875 ,  0.000 ,  0.000 ,  3.951 ,  1.681 , 17.353 , 28.273 , 19.422 , 21.724 , 22.164 , 14.902 , 14.458 , 21.174 , 30.024 , 20.949 , 15.101 , 28.926 , 23.850 , 30.597 , 10.147 , 20.035 , 18.873 , 14.691 , 13.892)

As duas amostras estão pareadas.
As duas amostras falham para o teste de normalidade:

> shapiro.test ( x )

	Shapiro-Wilk normality test

data:  testetcompounds$tR...tM.UV.ABC
W = 0.9373, p-value = 0.01057

> shapiro.test ( y )

	Shapiro-Wilk normality test

data:  testetcompounds$tR...tM.UV.ABC.Group10
W = 0.9372, p-value = 0.01043

No livro Bayesian Computation with R

Sugere que para distribuições não normais é necessário fazer o teste de Monte Carlo 

Com a seguinte função
=======================================================
 tstatistic = function ( y , x )
  {
  m = length ( y )
  n = length ( x )
  sp = sqrt( ( ( m - 1 ) * sd ( y ) ^ 2 + ( n - 1 ) * sd ( x ) ^ 2 ) / ( m + n - 2 ) )
  t = ( mean ( y ) - mean ( x ) ) / ( sp * sqrt ( 1 / m + 1 / n) )
  return ( t )
  }
========================================================
Retornando

t = -0.02242564

Em seguida o teste:
========================================================
alpha = 0.01 ; m = 50 ; n = 50 	# sets alpha, m and n ( alpha is the stated significance level, m and n are the samples sizes )
N = 10000                     			# sets the number of simulations
n.reject = 0                 			# counter of num. of rejections
for ( i in 1 : N )
  {
  x = rnorm ( m , mean = x , sd = x )     																		# simulates xs from population 1
  y = rnorm ( n , mean = y , sd = y )   																		# simulates ys from population 2
 t = tstatistic ( y , x )              	 	# computes the t statistic
  if ( abs ( t ) > qt ( 1 - alpha / 2 , n + m - 2) )
    n.reject = n.reject + 1               																			# reject if |t| exceeds critical pt
  }
true.sig.level = n.reject/N               																			# est. is proportions of rejections
========================================================
Retornando

true.sig.level = 0.0049


Problema Conceitual 

1) Como interpretar esse resultado?


Porém ainda há uma diferença entre o t calculado por essa fórmula e o t calculado com o comando t.stat (padrão do R)

================================================================
t = t.test ( y , x , mu = 0 , paired = TRUE , var.equal = TRUE , conf.level = 0.9999 )
t

	Paired t-test

data:  y and x
t = 4.6437, df = 49, p-value = 2.594e-05
alternative hypothesis: true difference in means is not equal to 0
99.99 percent confidence interval:
 0.003082235 0.067077765
sample estimates:
mean of the differences 
                0.03508 

================================================================

E finalmente… minha média da diferencia é conhecida e não é zero, assim utilizei o seguinte comando:

t = t.test ( y , x , mu = 0.87 , paired = TRUE , var.equal = TRUE , conf.level = 0.9999, alternative = "less” )

Paired t-test

data:  y and x
t = -110.5231, df = 49, p-value < 2.2e-16
alternative hypothesis: true difference in means is less than 0.87
99.99 percent confidence interval:
       -Inf 0.06545182
sample estimates:
mean of the differences 
                0.03508 

2) Problema de Comando… como fazer a simulação de Monte Carlo na situação acima?

* Obs.: aceito sugestões de leitura sobre o assunto e comando alternativos para comparar a médias


Muito obrigado e me desculpe pelo longo e-mail.


Thales