Usando o R para Inferência Estatística


26 de novembro de 2007

Inferência Estatística

Os procedimentos de inferência estatística analisam dados de uma amostra e desenham conclusões que sejam válidas para um universo maior de elementos chamado de população.

Ao analisar a amostra o interesse final é produzir um modelo probabilístico ou preditivo.

Uma boa amostra deve representar bem a população e , em geral, as amostras aleatórias conseguem fazer isto.

O principio da aleatoriedade traz consigo o princípio da independência. Em termos probabilisticos, as observações de uma amostra são realizações independentes de uma variável aleatória.

Em outras palavras, o valor registrado para uma variável no i-ésimo elemento selecionado não deve "influenciar"o valor do (i+1)-ésimo.

O Tamanho de Amostra também é importante para produzir boa inferência estatística.

Considere o seguinte exemplo abaixo que contém uma população de estudantes na faixa dos 10 aos 18 anos e os seguintes dados:

  > # Leitura dos dados de obesidade em adolescentes
  > ob<-read.csv("http://www.leg.ufpr.br/~joel/dados/obesidade.csv")

Neste caso, não seria necessário o uso de inferência estatística pois conhecemos todas as informações de interesse na população.

Logo a descrição das variáveis na população pode ser obtida pelo comando summary.

  > # resumindo as informações
  > summary(ob)

         ID       SEX           AGE             HEIGHT          WEIGHT
   Min.   :   1   F:2643   Min.   : 62.61   Min.   : 53.0   Min.   : 13.50
   1st Qu.:1251   M:2357   1st Qu.:146.62   1st Qu.:150.0   1st Qu.: 40.00
   Median :2500            Median :164.19   Median :158.0   Median : 49.00
   Mean   :2500            Mean   :167.16   Mean   :157.6   Mean   : 49.72
   3rd Qu.:3750            3rd Qu.:183.67   3rd Qu.:165.0   3rd Qu.: 57.00
   Max.   :5000            Max.   :667.77   Max.   :194.0   Max.   :128.00

Entretanto, vamos utilizar estes dados populacionais para alguns exercícios de amostragem e posteriormente inferência estatística.

Amostra Aleatória Simples

A primeira coluna, chamada de ID contém uma identificação do estudante e está numerada de 1 até 5000. Vamos retirar uma amostra de 50 estudantes e armazená-la.

  > # Selecionando os elementos da amostra
  > set.seed(1)
  > id.amo<-sample(ob$ID,50)
  > # Criando subconjunto com informações da amostra
  > ob.amo50<-ob[id.amo,]
  > # Exibição dos 10 primeiros elementos da amostra
  > head(ob.amo50,10)

         ID SEX    AGE HEIGHT WEIGHT
  1328 1328   F 138.47  146.0   57.0
  1861 1861   F 168.17  166.0   57.0
  2864 2864   F 204.01  159.0   45.0
  4539 4539   M 144.48  152.0   55.0
  1008 1008   F 155.49  147.0   37.0
  4488 4488   F 181.60  156.0   52.0
  4718 4718   F 146.12  155.0   46.0
  3300 3300   F 184.03  155.5   53.0
  3141 3141   M 140.67  145.0   33.4
  309   309   F 133.67  159.0   65.0

Construimos agora o resumo da amostra. Compare com o da população.

  > #resumo da amostra de tamanho 50
  > summary(ob.amo50)

         ID       SEX         AGE            HEIGHT          WEIGHT
   Min.   :  67   F:32   Min.   :129.4   Min.   :139.0   Min.   :33.40
   1st Qu.:1734   M:18   1st Qu.:147.9   1st Qu.:150.0   1st Qu.:43.00
   Median :2803          Median :158.9   Median :156.0   Median :48.00
   Mean   :2650          Mean   :165.6   Mean   :156.7   Mean   :51.60
   3rd Qu.:3786          3rd Qu.:180.3   3rd Qu.:162.8   3rd Qu.:57.75
   Max.   :4943          Max.   :223.4   Max.   :177.0   Max.   :85.00

Vamos analisar agora a variável altura em cm (HEIGHT). Para facilitar a manipulação, vamos criar um objeto para altura na população e outro para altura na amostra de tamanho n=50.

  > # Altura na populacao
  > ALTURA<-ob$HEIGHT
  > # Altura na amostra
  > altura<-ob.amo50$HEIGHT


  > par(mfrow=c(2,1))
  > # titulo do primeiro grafico
  > tit1<-'ALTURA NA POPULACAO'
  > # titulo do segundo grafico
  > tit2<-'altura na amostra'
  > hist(ALTURA,prob=T,xlim=c(100,200),main=tit1)
  > hist(altura,prob=T,xlim=c(100,200),main=tit2)

pict

Figura 1: Histograma da População e Amostra de Pesos (n=50)


A verdadeira média de altura é igual a 157.58 cm, enquanto para a amostra de tamanho 50, esta foi igual a 156.69. A princípio a estimativa pontual está bem próxima do valor verdadeiro na população.

Como utilizaremos as verdadeiras média e variância da variável altura mais adiante, vamos armazenar estas duas quantidades.

  > #Média da variável altura
  > mu<-mean(ALTURA)
  > mu
  [1] 157.5843
  > #Variância da variável altura
  > sigma2<-var(ALTURA)
  > sigma2
  [1] 125.9230
  > #Desvio padrão
  > sigma<-sqrt(sigma2)
Distribuição Amostral

A informação da estimativa pontual representa um único "chute"sobre o verdadeiro valor do parâmetro. Entretanto, selecionar a amostra trata-se de um experimento aleatório com outros possíveis resultados para as estatísticas de interesse. Como utilizamos a média amostral ¯x para estimar a média populacional, é importante saber como ela varia de amostra para amostra, ou seja, conhecer sua distribuição amostral.

Considere agora a seleção de 100 amostras de tamanho 50 em vez de apenas uma e o registro do variável altura em cada seleção feita.

  > # indices da populacao
  > i.p<-ob$ID
  > ###
  >
  > # 5000 amostras
  > # indices das amostras
  > i.a<-sample(i.p,50*100,rep=T)
  > ###
  >
  > # Matriz com 100 amostras de tamanho 50
  > matriz.alturas<-matrix(ALTURA[i.a],ncol=100)
  > ###
  >
  > # Médias das 100 amostras de tamanho 50
  > xbarra<-apply(matriz.alturas,2,mean)
  > xbarra

    [1] 157.88 161.75 157.79 158.01 157.64 160.90 156.21 156.23 157.33 152.99
   [11] 154.03 158.66 157.58 157.62 158.82 159.10 157.89 158.28 157.00 159.15
   [21] 158.78 156.48 158.68 156.31 158.52 157.48 158.57 155.38 156.69 153.70
   [31] 156.65 157.27 159.82 156.50 159.33 155.48 158.69 159.29 157.55 159.34
   [41] 157.66 156.50 157.90 158.78 159.30 157.38 153.46 158.83 157.55 157.00
   [51] 157.01 157.91 157.49 156.11 154.92 153.86 155.24 159.74 160.02 160.56
   [61] 159.24 157.32 156.06 158.09 158.43 158.12 157.16 159.94 159.20 158.47
   [71] 156.11 155.69 158.92 160.78 153.86 161.82 156.02 157.50 155.87 155.77
   [81] 157.87 162.88 157.19 157.30 156.13 156.94 160.36 157.19 160.73 157.45
   [91] 159.27 154.60 154.57 159.25 154.76 157.67 156.89 159.36 154.69 158.17

Para visualizar, empiricamente, a distribuição amostral da média, vamos usar novamente o histograma como ferramenta.


  > hist(xbarra,prob=T,main=' Histograma de 100 medias amostrais(n=50)')
  > rug(xbarra)
  > lines(density(xbarra),col='red')

pict

Figura 2: Ilustração da Distribuição Amostral de X¯


Vamos resumir a informações sobre esta amostra de valores da média de amostras de tamanho 50 da população de estudantes.

  > summary(xbarra)
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    153.0   156.4   157.6   157.6   158.8   162.9

Novamente vale ressaltar que o resultado médio encontrado acima está bem próximo do verdadeiro valor populacional igual a 157.58 cm.

Vamos medir a variabilidade da média de amostras de tamanho 50 da população de estudantes, pelo desvio padrão. O valor resultante desta operação é uma estimativa do verdadeiro erro padrão de ¯x para amostras de tamanho 50. Este valor é igual a σ
√
n= 11.22
√
50=1.59

  > sd(xbarra)
  [1] 1.917299

Todo o procedimento acima, tem interesse na estimação da média. Entretanto, em muitos casos, o interesse recai sobre a variância ou desvio-padrão.

O melhor estimador para a variância é S2 = ∑  (xi - x¯)2
--i--------
   n- 1

  > s2<-apply(matriz.alturas,2,var)
  > s2
    [1]  96.46490 152.59439 113.31724 141.54582 116.70449 137.23469 107.80704
    [8] 131.74704 116.15929 294.65806  86.12663 105.03510 149.53429 104.63837
   [15] 120.34449 109.68367 103.92133 113.79755  82.11224 190.69643 124.95061
   [22]  81.03020 122.67102  91.50908 123.71388 111.92816 124.71439 105.65878
   [29] 121.82541 102.54082 123.50255 130.13480  85.91592 230.41837 116.51643
   [36] 206.32612 102.25398 118.39888  98.77806 146.02490 116.36163  85.54082
   [43]  82.19388 110.63429  97.44898  93.76082 133.27388  92.41439 111.03316
   [50] 133.16327 129.88255 133.47643 135.27031  92.56418  93.97306 155.17388
   [57]  90.15551 103.63510 133.12204 120.45551 112.39020 125.48735 123.45551
   [64] 127.50704 115.74500 107.13837 111.82082  88.46571 135.52041 127.92255
   [71] 153.62541 220.70296 113.41184 183.74653 315.43918 146.88531 122.82612
   [78] 152.30612 107.97765  97.54296 146.51847 109.52612 137.90704 137.14286
   [85] 121.07969 169.99633 148.51061 122.10092 121.60418 162.41071 111.29806
   [92] 197.10204  78.41847 105.63520  97.46163 114.68990 124.54378 141.73510
   [99] 174.47847 119.57765

Repare com a evidência empírica que a distribuição amostral desta estatística (S2) tem formato completamente diferente da distribuição da média amostral.


  > hist(s2,prob=T,ylim=c(0,0.018),main="Histograma de 100 variancias amostrais (n=50)")
  > rug(s2)
  > lines(density(s2),col='red')

pict

Figura 3: Ilustração da Distribuição Amostral da Variância


Veja o resumo para os valores de S2 encontrados neste experimento.

  > summary(s2)
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    78.42  105.50  120.40  127.20  135.90  315.40

A variância populacional é igual a 125.92 cm2.

Intervalos de Confiança

A importante interpretação do intervalo de 100γ% de confiança é a de que o verdadeiro valor do parâmetro estará contido dentro deste intervalo em 100γ% das amostras de tamanho n que forem retiradas.1

Intervalo para a Média

Vamos utilizar o próprio exemplo das alturas para ilustrar a idéia por trás do intervalo de confiança. Nas 100 amostras de tamanho 50 da variável altura, armazenadas no objeto matriz.alturas, vamos calcular em cada uma destas amostras o intervalo de 95% de confiança para a verdadeira média(admitindo a variância conhecida).

  > #Limites inferiores
  > lmi<-function(x){mean(x)+qnorm(0.025)*sigma/sqrt(length(x))}
  > LI<-apply(matriz.alturas,2,lmi)
  > LI

    [1] 154.7696 158.6396 154.6796 154.8996 154.5296 157.7896 153.0996 153.1196
    [9] 154.2196 149.8796 150.9196 155.5496 154.4696 154.5096 155.7096 155.9896
   [17] 154.7796 155.1696 153.8896 156.0396 155.6696 153.3696 155.5696 153.1996
   [25] 155.4096 154.3696 155.4596 152.2696 153.5796 150.5896 153.5396 154.1596
   [33] 156.7096 153.3896 156.2196 152.3696 155.5796 156.1796 154.4396 156.2296
   [41] 154.5496 153.3896 154.7896 155.6696 156.1896 154.2696 150.3496 155.7196
   [49] 154.4396 153.8896 153.8996 154.7996 154.3796 152.9996 151.8096 150.7496
   [57] 152.1296 156.6296 156.9096 157.4496 156.1296 154.2096 152.9496 154.9796
   [65] 155.3196 155.0096 154.0496 156.8296 156.0896 155.3596 152.9996 152.5796
   [73] 155.8096 157.6696 150.7496 158.7096 152.9096 154.3896 152.7596 152.6596
   [81] 154.7596 159.7696 154.0796 154.1896 153.0196 153.8296 157.2496 154.0796
   [89] 157.6196 154.3396 156.1596 151.4896 151.4596 156.1396 151.6496 154.5596
   [97] 153.7796 156.2496 151.5796 155.0596

  > #Limites superiores
  > lms<-function(x){mean(x)+qnorm(0.975)*sigma/sqrt(length(x))}
  > LS<-apply(matriz.alturas,2,lms)
  > LS

    [1] 160.9904 164.8604 160.9004 161.1204 160.7504 164.0104 159.3204 159.3404
    [9] 160.4404 156.1004 157.1404 161.7704 160.6904 160.7304 161.9304 162.2104
   [17] 161.0004 161.3904 160.1104 162.2604 161.8904 159.5904 161.7904 159.4204
   [25] 161.6304 160.5904 161.6804 158.4904 159.8004 156.8104 159.7604 160.3804
   [33] 162.9304 159.6104 162.4404 158.5904 161.8004 162.4004 160.6604 162.4504
   [41] 160.7704 159.6104 161.0104 161.8904 162.4104 160.4904 156.5704 161.9404
   [49] 160.6604 160.1104 160.1204 161.0204 160.6004 159.2204 158.0304 156.9704
   [57] 158.3504 162.8504 163.1304 163.6704 162.3504 160.4304 159.1704 161.2004
   [65] 161.5404 161.2304 160.2704 163.0504 162.3104 161.5804 159.2204 158.8004
   [73] 162.0304 163.8904 156.9704 164.9304 159.1304 160.6104 158.9804 158.8804
   [81] 160.9804 165.9904 160.3004 160.4104 159.2404 160.0504 163.4704 160.3004
   [89] 163.8404 160.5604 162.3804 157.7104 157.6804 162.3604 157.8704 160.7804
   [97] 160.0004 162.4704 157.8004 161.2804

Armazenando os limites de confiança, verificamos se o verdadeiro valor da média está contido neles.

  > ic.media<-data.frame(cbind(LI,LS))
  > head(ic.media)

          LI       LS
  1 154.7696 160.9904
  2 158.6396 164.8604
  3 154.6796 160.9004
  4 154.8996 161.1204
  5 154.5296 160.7504
  6 157.7896 164.0104

  > tail(ic.media)

            LI       LS
  95  151.6496 157.8704
  96  154.5596 160.7804
  97  153.7796 160.0004
  98  156.2496 162.4704
  99  151.5796 157.8004
  100 155.0596 161.2804

Caso a variância não seja conhecida, o intervalo de confiança acima é calculado substituindo quantis da normal(-1,96; 1,96) por quantis da distribuição de t-student com 49 graus de liberdade (-2.01,2.01)

Abaixo, contamos a freqüência de vezes que o verdadeiro valor da média (157.58) está no interior do intervalo.

  > table(ic.media$LI<mu &  ic.media$LS>mu)

  FALSE  TRUE
     12    88

Intervalo para a Variância

A situação mais real é aquela em que além da média, também não conhecemos o valor populacional da variância. Os passos abaixo estimarão intervalos para a variância, baseados na distribuição Qui-quadrado.

  > #Limites inferiores
  > lmi<-function(x){(length(x)-1)*var(x)/qchisq(0.975,length(x)-1)}
  > LI<-apply(matriz.alturas,2,lmi)
  > LI

    [1]  67.31156 106.47776  79.07084  98.76825  81.43440  95.76002  75.22591
    [8]  91.93083  81.05396 205.60736  60.09769  73.29170 104.34247  73.01486
   [15]  83.97433  76.53539  72.51453  79.40599  57.29652 133.06471  87.18840
   [22]  56.54149  85.59774  63.85347  86.32543  78.10156  87.02357  73.72689
   [29]  85.00769  71.55123  86.17797  90.80584  59.95066 160.78200  81.30317
   [36] 143.97084  71.35108  82.61671  68.92564 101.89368  81.19516  59.68892
   [43]  57.35348  77.19871  67.99823  65.42470  92.99623  64.48518  77.47704
   [50]  92.91905  90.62982  93.13757  94.38931  64.58971  65.57280 108.27768
   [57]  62.90897  72.31480  92.89028  84.05180  78.42396  87.56293  86.14515
   [64]  88.97223  80.76488  74.75932  78.02665  61.72986  94.56382  89.26217
   [71] 107.19718 154.00275  79.13684 128.21519 220.10807 102.49406  85.70597
   [78] 106.27661  75.34496  68.06381 102.23808  76.42546  96.22918  95.69594
   [85]  84.48734 118.62053 103.62817  85.19993  84.85332 113.32742  77.66188
   [92] 137.53444  54.71907  73.71044  68.00706  80.02865  86.90452  98.90033
   [99] 121.74809  83.43924

  > #Limites superiores
  > lms<-function(x){(length(x)-1)*var(x)/qchisq(0.025,length(x)-1)}
  > LS<-apply(matriz.alturas,2,lms)
  > LS

    [1] 149.7954 236.9559 175.9645 219.7992 181.2244 213.1047 167.4080 204.5832
    [9] 180.3778 457.5593 133.7416 163.1036 232.2041 162.4875 186.8767 170.3221
   [17] 161.3741 176.7103 127.5079 296.1226 194.0294 125.8276 190.4895 142.0997
   [25] 192.1089 173.8075 193.6625 164.0721 189.1764 159.2303 191.7807 202.0796
   [33] 133.4144 357.8048 180.9323 320.3932 158.7849 183.8555 153.3873 226.7545
   [41] 180.6920 132.8319 127.6346 171.7983 151.3235 145.5963 206.9541 143.5055
   [49] 172.4177 206.7824 201.6879 207.2687 210.0543 143.7381 145.9259 240.9615
   [57] 139.9978 160.9296 206.7183 187.0491 174.5249 194.8628 191.7077 197.9991
   [65] 179.7344 166.3696 173.6408 137.3738 210.4426 198.6443 238.5570 342.7182
   [73] 176.1114 285.3305 489.8292 228.0906 190.7303 236.5083 167.6729 151.4694
   [81] 227.5210 170.0775 214.1487 212.9621 188.0184 263.9785 230.6145 189.6042
   [89] 188.8329 252.1992 172.8290 306.0696 121.7720 164.0355 151.3431 178.0960
   [97] 193.3976 220.0931 270.9386 185.6860

Veja agora os limites de confiança e verifique se o verdadeiro valor da variância está contido nele.

  > ic.var<-data.frame(cbind(LI,LS))
  > head(ic.var)

           LI       LS
  1  67.31156 149.7954
  2 106.47776 236.9559
  3  79.07084 175.9645
  4  98.76825 219.7992
  5  81.43440 181.2244
  6  95.76002 213.1047

  > tail(ic.var)

             LI       LS
  95   68.00706 151.3431
  96   80.02865 178.0960
  97   86.90452 193.3976
  98   98.90033 220.0931
  99  121.74809 270.9386
  100  83.43924 185.6860

Tal como feito para a média, contamos a freqüência de vezes que o verdadeiro valor da variância (125.92) está no interior do intervalo.

  > table(ic.var$LI<sigma2 &  ic.var$LS>sigma2)

  FALSE  TRUE
     10    90