Notebook aula bioestatística

Teste de Hipótese

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introdução aos testes de hipóteses\n",
    "\n",
    "### Reginaldo K Fukuchi\n",
    "\n",
    "Com frequência lemos em artigos científicos descrições de natureza similar a esta.\n",
    "\n",
    "“With a two-tailed t-test, the null hypothesis was rejected at the 5% alpha-level. The results were statistically significant with a p-value of 0.018”\n",
    "\n",
    "\n",
    "Hoje iremos apresentar conceitos sobre o teste de hipótese. O teste de hipótese é um teste estatístico que usa dados obtidos de uma amostra para tomar decisões a respeito de aceitar ou não a **hipótese nula** ($H_0$). \n",
    "\n",
    "Embora estamos tratanto de estatística, podemos considerar esta tomada de decisão como a mesma decisão que um juiz deve tomar ao condenar ou absolver um réu. Esta decisão é baseada nas evidências apresentadas. Tipicamente, pressupomos que todos somos inocentes até que provem o contrário. Na estatística, também tipicamente assumimos, pelo menos temporariamente, que a $H_0$ é verdadeira. Então com base nas evidências decidimos rejeitá-la ou não em favor da **hipótese alternativa** ($H_1$). \n",
    "\n",
    "Portanto,\n",
    "* $H_0$ consiste declarar estatisticamente que não existe diferença, resposta ou mudança.\n",
    "* $H_1$ consiste em contradizer a $H_0$\n",
    "\n",
    "Por exemplo, uma suposição é feita sobre uma determinada característica da população que pode ou não ser verdadeira. Então, tudo se inicia com uma boa pergunta. Exemplos:\n",
    "* A renda média do brasileiro diminiu em relação ao ano passado?\n",
    "* A média de temperatura caiu em relação ao inverno do ano passado?\n",
    "* A proporção de mulheres matriculadas aumentou?\n",
    "\n",
    "Para responder estas e outras questões recorremos então aos testes estatísticos.\n",
    "\n",
    "A melhor forma de estudar os testes de hipóteses é por meio de exemplos. Então vamos usar alguns exemplos de (Bluman, 2012)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Teste Z para uma amostra\n",
    "Example. A researcher claims that the average cost of men’s athletic shoes is less than 80 dollars. He selects a random sample of 36 pairs of shoes from a catalog and finds the following costs (in dollars). (The costs have been rounded to the nearest dollar.) Is there enough evidence to support the researcher’s claim at a $\\alpha=$ 0.10? Assume $\\sigma=$ 19.2.\n",
    "\n",
    "Baseado na descrição do problema acima, podemos notar que o desvio padrão da população ($\\sigma=$) é conhecido e que o tamanho da amostra (*n* > 30). Então podemos usar o teste Z para responder esta questão.\n",
    "\n",
    "Ainda que existam particularidades, a formulação geral de um teste estatístico segue a seguinte convenção.\n",
    "\n",
    "$$\\mathbf{Teste} = \\frac{\\mathbf{Observado} - \\mathbf{Esperado}}{\\mathbf{Erro Padrão}} $$\n",
    "\n",
    "Para o teste Z, podemos formular assim.\n",
    "\n",
    "$$\\mathbf{Z} = \\frac{\\mathbf{\\bar{X}} - \\mathbf{\\mu}}{\\frac{\\mathbf{\\sigma}}{\\mathbf{\\sqrt{n}}}} $$\n",
    "\n",
    "Onde,\n",
    "\n",
    "$\\mathbf{\\bar{X}}$: média da amostra\n",
    "\n",
    "$\\mathbf{\\mu}$: média da população\n",
    "\n",
    "$\\mathbf{\\sigma}$: desvio padrão da população\n",
    "\n",
    "*n*: tamanho da amostra\n",
    "\n",
    "Antes de resolver este problema, precisamos formular as $H_0$ e $H_1$. Vejamos que a palavra \"*less*\" foi mencionada. Então, $H_1$ é unicaudal e portanto.\n",
    "\n",
    "$$ \\mathbf{H_0} : \\mathbf{\\mu} = 80$$ $$\\mathbf{H_1} : \\mathbf{\\mu} < 80$$ \n",
    "\n",
    "Perceba que se a palavra *different* fosse usada no lugar de *less*, as hipóteses seriam.\n",
    "\n",
    "$$ \\mathbf{H_0} : \\mathbf{\\mu} = 80$$ $$\\mathbf{H_1} : \\mathbf{\\mu} \\neq 80$$ \n",
    "\n",
    "Portanto, existem diferentes frases que são usadas para expressar as hipóteses do estudo como mostra a Figura abaixo.\n",
    "\n",
    "<figure><img src=\"http://4.bp.blogspot.com/-INU3DejOWOc/UHmQYBWA7eI/AAAAAAAAAQw/orlLDkXnL7k/s1600/jj.jpg\" alt=\"Torque\" width=\"400\"/><figcaption><center><i>Figura. Frases comumente empregadas em testes de hipóteses </i></center></figcaption></figure>\n",
    "\n",
    "Outra informação importante descrita no problema foi $\\alpha$ = 0.10 que é considerado o **nível de significância** do teste estatístico. O nível de significância é a probabilidade de rejeitar $H_0$ quando ela é verdadeira que também é conhecido como **erro tipo I**. O valor $\\alpha$ = 0.05 tem sido adotado pela comunidade científica e sua origem se deveu ao fato que Ronald Fisher considerava conveniente uma probabilidade de 1 em 20 de rejeitar erroneamente a $H_0$.\n",
    "\n",
    "Assim, como um juiz pode culpar um inocente também corremos o risco de rejeitar uma $H_0$ verdadeira (**erro tipo I**). Ao mesmo tempo, podemos também aceitar (ou não rejeitar) uma $H_0$ quando ela não é verdadeira. Neste útimo caso estamos cometendo o **erro tipo II**. De fato existem quatro possibilidades (duas decisões corretas e duas incorretas) quando tomamos uma decisão baseado no teste estatístico como ilustra a Figura abaixo.\n",
    "\n",
    "<figure><img src=\"http://allpsych.com/wp-content/uploads/2014/08/type1and2error.gif\" alt=\"Torque\" width=\"400\"/><figcaption><center><i>Figura. Tabela da verdade do teste de hipóteses </i></center></figcaption></figure>\n",
    "\n",
    "Agora que entedemos algumas questões fundamentais podemos resolver o problema original. Vamos usar o R para isto."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Import the necessary libraries\n",
    "library(\"psych\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Raw data\n",
    "shoePrices <- c(60,70,75,55,80,55,50,40,80,70,50,95,120,90,75,85,80,60,110,65,80,85,85,45,75,60,90,90,60,95,110,85,45,90,70,70)\n",
    "\n",
    "\n",
    "mu <- 80 # population mean\n",
    "sPricesM <- mean(shoePrices) # Mean shoe prices\n",
    "sigma <- 19.2 # population SD\n",
    "alpha <- 0.10 # level of significance\n",
    "n <- 36 # sample size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "-1.5625"
      ],
      "text/latex": [
       "-1.5625"
      ],
      "text/markdown": [
       "-1.5625"
      ],
      "text/plain": [
       "[1] -1.5625"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "-1.2815515655446"
      ],
      "text/latex": [
       "-1.2815515655446"
      ],
      "text/markdown": [
       "-1.2815515655446"
      ],
      "text/plain": [
       "[1] -1.281552"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Z test\n",
    "ztest <- (sPricesM - mu)/(sigma/sqrt(n)) \n",
    "ztest\n",
    "# Finding critical value for alpha = 0.05 two-tailed test\n",
    "cv <- qnorm(alpha, mean = 0, sd = 1, lower.tail = TRUE)\n",
    "cv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Como a hipótese foi direcional, vamos usar o teste unicaudal. Como o valor do teste é menor que o valor crítico, a hipótese nula pode ser rejeitada. Portanto, podemos concluir que existe evidência suficiente para indicar que o custo médio de um calçado é **menor** que U$80.\n",
    "\n",
    "Dependendo da formulação da hipótese estatística, o teste estatístico pode ser **unicaudal** ou **bicaudal**. O teste unicaudal é aquele onde a região de rejeição só pode estar a direita (maior) ou a esquerda da média (menor). Por outro lado o bicaudal a região de rejeição pode se encontrar em qualquer lado ($\\neq$) da média. Portanto, a hipótese deve ser formulada **antes** de ter acesso aos dados do estudo. A figura abaixo ilustra as diferentes condições de teste.\n",
    "\n",
    "<figure><img src=\"https://www.cliffsnotes.com/assets/267184.png\" alt=\"Torque\" width=\"350\"/><figcaption><center><i>Figura. Testes estatítico unicaudal e bicaudal </i></center></figcaption></figure>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "### Teste t para uma amostra\n",
    "\n",
    "Veja a descrição de um outro problema sobre infecções hospitalares (Bluman, 2012).\n",
    "\n",
    "\"*A medical investigation claims that the average number of infections per week at a\n",
    "hospital in southwestern Pennsylvania is 16.3. A random sample of 10 weeks had a\n",
    "mean number of 17.7 infections. The sample standard deviation is 1.8. Is there enough\n",
    "evidence to reject the investigator’s claim at $\\alpha$ = 0.05?*\"\n",
    "\n",
    "De acordo com a descrição podemos notar que:\n",
    "* o tamanho da amostra é pequena (*n*=10)\n",
    "* o desvio padrão da população ($\\sigma$) é desconhecido\n",
    "\n",
    "Vimos anteriormente que para as condições acima, não é recomendado usar a distribuição Z mas sim uma outra família chamada distribuição t de Student. Portanto, para resolver este problema vamos usar o teste t com a seguinte formulação.\n",
    "\n",
    "\n",
    "$$\\mathbf{t} = \\frac{\\mathbf{\\bar{X}} - \\mathbf{\\mu}}{\\frac{\\mathbf{s}}{\\mathbf{\\sqrt{n}}}} $$\n",
    "\n",
    "Onde,\n",
    "\n",
    "$\\mathbf{\\bar{X}}$: média da amostra\n",
    "\n",
    "$\\mathbf{\\mu}$: média da população\n",
    "\n",
    "$\\mathbf{s}$: desvio padrão da amostra\n",
    "\n",
    "*n*: tamanho da amostra\n",
    "\n",
    "g.l.: graus de liberdade (n-1)\n",
    "\n",
    "De acordo com o problema as seguintes hipóteses são formuladas:\n",
    "$$ \\mathbf{H_0} : \\mathbf{\\mu} = 16.3$$ $$\\mathbf{H_1} : \\mathbf{\\mu} \\neq 16.3$$ \n",
    "\n",
    "Portanto, o teste é **bicaudal**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# presented data\n",
    "mu <- 16.3 # avg population\n",
    "x  <- 17.7 # avg sample\n",
    "s  <- 1.8 # sd sample\n",
    "n  <- 10 # sample size\n",
    "df <- n-1 # dof\n",
    "alpha <- 0.05"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "2.46"
      ],
      "text/latex": [
       "2.46"
      ],
      "text/markdown": [
       "2.46"
      ],
      "text/plain": [
       "[1] 2.46"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>-2.262</li>\n",
       "\t<li>2.262</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item -2.262\n",
       "\\item 2.262\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. -2.262\n",
       "2. 2.262\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] -2.262  2.262"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# t test\n",
    "ttest <- (x - mu)/(s/sqrt(n))\n",
    "round(ttest,3)\n",
    "\n",
    "cv <- qt(alpha/2, df)\n",
    "round(cv*c(1,-1),3) # critical values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Como o valor do teste é maior que o valor crítico, rejeita-se $H_0$. Portanto, conclui-se que o número médio de infecções hospitalares por semana é diferente da suspeitada. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparação entre duas médias com teste Z\n",
    "\n",
    "Em alguns casos (para não dizer a maioria) estamos interessados em comparar duas amostras como por exemplo um grupo que fez tratamento e o outro controle.\n",
    "\n",
    "De novo vamos usar um exemplo do livro.\n",
    "\n",
    "\"*A survey found that the average hotel room rate in New Orleans is \\$ 88.42 and the\n",
    "average room rate in Phoenix is \\$ 80.61. Assume that the data were obtained from two\n",
    "samples of 50 hotels each and that the standard deviations of the populations are \\$ 5.62\n",
    "and \\$ 4.83, respectively. At $\\alpha$ =0.05, can it be concluded that there is a significant\n",
    "difference in the rates?*\"\n",
    "\n",
    "Considerando que ambas as amostras sejam **aleatórias**, **independentes** e que o **desvio padrão da população**, das quais elas pertencem, seja conhecido podemos usar o teste Z. Note também que o tamanho da amostra é maior que 30. \n",
    "\n",
    "Mas agora temos duas amostras e precisamos formular as hipóteses um pouco diferente das que foram feitas anteriormente como segue.\n",
    "\n",
    "$$ \\mathbf{H_0} : \\mathbf{\\mu_1}  = \\mathbf{\\mu_2}$$ $$\\mathbf{H_1} : \\mathbf{\\mu_1}  \\neq \\mathbf{\\mu_2}$$\n",
    "\n",
    "Ou, alternativamente\n",
    "\n",
    "$$ \\mathbf{H_0} : \\mathbf{\\mu_1}  - \\mathbf{\\mu_2} = 0$$ $$\\mathbf{H_1} : \\mathbf{\\mu_1} - \\mathbf{\\mu_2} \\neq 0$$\n",
    "\n",
    "Então, a distribuição agora será de diferenças entre as médias de pares de amostras retiradas da população. Se as populações tiverem médias iguais, esta média das médias das diferenças também será zero. A variância da diferença $\\bar{X_1} - \\bar{X_2}$ é igual a soma das variâncias individuais,\n",
    "\n",
    "$$ \\mathbf{\\sigma^2}_{\\bar{X_1} - \\bar{X_2}} = \\sigma^2_{\\bar{X}_1} + \\sigma^2_{\\bar{X}_2} $$\n",
    "\n",
    "E então,\n",
    "\n",
    "$$ \\mathbf{\\sigma^2}_{\\bar{X_1} - \\bar{X_2}} = \\frac{\\sigma^2_1}{n_1} + \\frac{\\sigma^2_2}{n_2}$$\n",
    "\n",
    "portanto, o erro padrão da diferença $\\bar{X_1} - \\bar{X_2}$ é\n",
    "\n",
    "$$ \\mathbf{\\sqrt{\\frac{\\sigma^2_1}{n_1} + \\frac{\\sigma^2_2}{n_2}}} $$\n",
    "\n",
    "Então, o teste Z para comparação entre amostras quando o desvio padrão da população é conhecido é calculado por\n",
    "\n",
    "$$ \\mathbf{z} = \\frac{(\\mathbf{\\bar{X}_1} - \\mathbf{\\bar{X}_2}) - (\\mathbf{\\mu_1} - \\mathbf{\\mu}_2)}{\\mathbf{\\sqrt{\\frac{\\sigma^2_1}{n_1} + \\frac{\\sigma^2_2}{n_2}}}} $$\n",
    "\n",
    "Resolvemos então o problema apresentado anteriormente no R."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Presented data\n",
    "x1 <- 88.42 # avg sample 1\n",
    "x2 <- 80.61 # avg sample 2\n",
    "sigma1 <- 5.62 # sd population 1\n",
    "sigma2 <- 4.83 # sd population 2\n",
    "n1 <- 50 # sample size 1\n",
    "n2 <- n1 # sample size 2\n",
    "alpha <- 0.05 # sig. level"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "7.452"
      ],
      "text/latex": [
       "7.452"
      ],
      "text/markdown": [
       "7.452"
      ],
      "text/plain": [
       "[1] 7.452"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>-1.96</li>\n",
       "\t<li>1.96</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item -1.96\n",
       "\\item 1.96\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. -1.96\n",
       "2. 1.96\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] -1.96  1.96"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Calculating Z test\n",
    "ztest <- ((x1-x2) - 0) / sqrt((sigma1^2/n1) + (sigma2^2/n2))\n",
    "round(ztest,3)\n",
    "\n",
    "# Finding critical value for two-tailed test\n",
    "cv <- qnorm(alpha/2, mean=0, sd = 1)\n",
    "round(cv*c(1,-1),3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Como o teste estatítico é maior que o valor crítico podemos rejeitar $H_0$ em favor da $H_1$ e concluir que as diárias dos hotéis nas duas cidades são diferentes. Note que não foi dado uma hipótese direcional."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparação entre duas médias com teste **t**\n",
    "#### Amostras independentes\n",
    "\n",
    "Na maioria dos casos vamos nos deparar com situações onde o desvio padrão das populações estudadas é desconhecido e apenas o das amostras pode ser obtido. Em adição, as amostras obtidas também são inferiores a 30. Nestes casos não podemos usar o teste Z pois violaremos as suas pressuposições. Vamos analisar o exemplo abaixo do livro.\n",
    "\n",
    "*\"The average size of a farm in Indiana County, Pennsylvania, is 191 acres. The average size\n",
    "of a farm in Greene County, Pennsylvania, is 199 acres. Assume the data were obtained\n",
    "from two samples with standard deviations of 38 and 12 acres, respectively, and sample\n",
    "sizes of 8 and 10, respectively. Can it be concluded at $\\alpha$ = 0.05 that the average size of the\n",
    "farms in the two counties is different? Assume the populations are normally distributed.\"*\n",
    "\n",
    "Como pode ser observado, os desvios padrões das populações não foram informados e o tamanho da amostra é pequeno. Vimos anteriormente que a família de distribuição *t* pode ser usado nestes casos pois o formato da distribuição varia em função do tamanho da amostra.\n",
    "\n",
    "A formulação do teste *t* para comparação de duas médias amostrais é similar ao do teste Z como segue.\n",
    "\n",
    "$$ \\mathit{t} = \\frac{(\\mathbf{\\bar{X}_1} - \\mathbf{\\bar{X}_2}) - (\\mathbf{\\mu_1} - \\mathbf{\\mu}_2)}{\\mathbf{\\sqrt{\\frac{s^2_1}{n_1} + \\frac{s^2_2}{n_2}}}} $$\n",
    "\n",
    "Note que a diferença foi o uso do desvio padrão das amostras. Se você recordar, para a distribuição *t* é necessário levar em conta o grau de liberdade (g.l. = n-1). No caso, como temos um tamanho de amostra diferente usamos o menor para ser conservador.\n",
    "\n",
    "É importante lembrar sempre de forumlar as hipóteses antes de realizar o teste estatítico. Para este caso, não houve informação de direção e consideramos hipótese bidirecional.\n",
    "\n",
    "$$ \\mathbf{H_0} : \\mathbf{\\mu_1}  - \\mathbf{\\mu_2} = 0$$ $$\\mathbf{H_1} : \\mathbf{\\mu_1} - \\mathbf{\\mu_2} \\neq 0$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Presented data\n",
    "x1 <- 191\n",
    "x2 <- 199\n",
    "n1 <- 8\n",
    "n2 <- 10\n",
    "s1 <- 38\n",
    "s2 <- 12\n",
    "alpha <- 0.05\n",
    "df1 <- n1-1\n",
    "df2 <- n2-1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "-0.573"
      ],
      "text/latex": [
       "-0.573"
      ],
      "text/markdown": [
       "-0.573"
      ],
      "text/plain": [
       "[1] -0.573"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>-2.365</li>\n",
       "\t<li>2.365</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item -2.365\n",
       "\\item 2.365\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. -2.365\n",
       "2. 2.365\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] -2.365  2.365"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Performing t test\n",
    "ttest <- ((x1-x2) - 0) / sqrt((s1^2/n1) + (s2^2/n2))\n",
    "round(ttest,3)\n",
    "\n",
    "# Calculating critical value\n",
    "if (df1 < df2) {df <- df1} else {df <- df2}\n",
    "\n",
    "cv <- qt(alpha/2,df)\n",
    "round(cv*c(1,-1),3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Como o valor do teste estatístico está dentro da região de aceitação, não existe evidência suficiente para rejeitar a $H_0$. Portanto, podemos concluir que não existe diferença na dimensão das fazendas dos dois condados."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Amostras dependentes\n",
    "Em algumas situações, é de interesse que a mesma amostra seja comparada em dois instantes (antes e depois) como por exemplo para testar o efeito de uma intervenção. Nestes casos, não podemos considerar valores independentes. Considere o seguinte exemplo abaixo retirado do livro.\n",
    "\n",
    "*\"A dietitian wishes to see if a person’s cholesterol level will change if the diet is\n",
    "supplemented by a certain mineral. Six subjects were pretested, and then they\n",
    "took the mineral supplement for a 6-week period. The results are shown in the table.\n",
    "(Cholesterol level is measured in milligrams per deciliter.) Can it be concluded that the\n",
    "cholesterol level has been changed at $\\alpha$ = 0.10? Assume the variable is approximately\n",
    "normally distributed.\"*\n",
    "\n",
    "Perceba no exemplo que as medidas foram obtidas da mesma amostra de sujeito e, portanto, são dependentes. A formulação das hipóteses segue o mesmo princípio que anteriormente.\n",
    "\n",
    "$$ \\mathbf{H_0} : \\mathbf{X_{pre}}  - \\mathbf{X_{pos}} = 0$$ $$\\mathbf{H_1} : \\mathbf{X_{pre}}  - \\mathbf{X_{pos}} \\neq 0$$\n",
    "\n",
    "Se considerarmos que $\\bar{D} = \\mathbf{X_{pre}}  - \\mathbf{X_{pos}}$, então a formulação do teste estatístico para condições dependentes é,\n",
    "\n",
    "$$ \\mathbf{t} = \\frac{\\mathbf{\\bar{D}} - \\mathbf{\\mu_D}}{\\frac{\\mathbf{s_D}}{\\sqrt{n}}}   $$\n",
    "\n",
    "Onde, $S_D$ é o desvio padrão das diferenças e o d.f. = n - 1.\n",
    "\n",
    "Vamos resolver o problema no R."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Raw data\n",
    "prior <- c(210,235,208,190,172,244)\n",
    "after <- c(190,170,210,188,173,228)\n",
    "\n",
    "n <- 6 # sample size\n",
    "df <- n-1 # dof\n",
    "alpha <- 0.1 # sig level\n",
    "# Calculated data\n",
    "dM <- mean(prior-after) # avg of differences\n",
    "sD <- sd(prior-after) # sd of differences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "1.608"
      ],
      "text/latex": [
       "1.608"
      ],
      "text/markdown": [
       "1.608"
      ],
      "text/plain": [
       "[1] 1.608"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Statistical test\n",
    "ttest <- (dM - 0)/(sD/sqrt(n))\n",
    "round(ttest,3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>-2.015</li>\n",
       "\t<li>2.015</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item -2.015\n",
       "\\item 2.015\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. -2.015\n",
       "2. 2.015\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] -2.015  2.015"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Critical value\n",
    "cv <- qt(alpha/2,df)\n",
    "round(cv*c(1,-1),3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Como o valor do teste está dentro da região de aceitação, não existe evidência para refutar a $H_0$ em favor de $H_1$ e, portanto, conclui-se que os níveis de colesterol não mudaram em função da dieta. Veja se o desenho do estudo permite tirar esta conclusão."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Poder estatístico\n",
    "\n",
    "<figure><img src=\"http://www.nature.com/nmeth/journal/v10/n12/images/nmeth.2738-F2.jpg\" alt=\"Torque\" width=\"640\"/><figcaption><center><i>Figura. (a) Observations are assumed to be from the null distribution (H0) with mean μ0. We reject H0 for values larger than x* with an error rate α (red area). (b) The alternative hypothesis (HA) is the competing scenario with a different mean μA. Values sampled from HA smaller than x* do not trigger rejection of H0 and occur at a rate β. Power (sensitivity) is 1 − β (blue area). (c) Relationship of inference errors to x*. </i></center></figcaption></figure>\n",
    "\n",
    "\n",
    "<figure><img src=\"http://www.nature.com/nmeth/journal/v10/n12/images_article/nmeth.2738-F4.jpg\" alt=\"Torque\" width=\"640\"/><figcaption><center><i>Figura. H0 and HA are assumed normal with σ = 1. (a) Increasing n decreases the spread of the distribution of sample averages in proportion to 1/√n. Shown are scenarios at n = 1, 3 and 7 for d = 1 and α = 0.05. Right, power as function of n at four different α values for d = 1. The circles correspond to the three scenarios. (b) Power increases with d, making it easier to detect larger effects. The distributions show effect sizes d = 1, 1.5 and 2 for n = 3 and α = 0.05. Right, power as function of d at four different a values for n = 3. </i></center></figcaption></figure>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Referências\n",
    "* Bluman, Allan G. Elementary statistics : a step by step approach / Allan Bluman. — 8th ed.\n",
    "* Krzywinski & Altman (2013). Points of significance: Power and sample size. Nature Methods 10, 1139–1140."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.3.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}