Archive

Archive for the ‘testes’ Category

Fuzzy c-means, um exemplo!

Pois bem, por sugestão do Professor Julio da Motta Singer apresentarei um script para realizar uma análise de agrupamento no escopo da lógica fuzzy. Seja um experimento fatorial (a mesma natureza dos dados aqui disponibilizados) a seguinte implementação realiza a clusterização de uma dado fator marginalmente ao outro fator.

A heurística da lógica fuzzy tem um incremento interpretativo de, ao contrário do agrupamento k-means, que é dicotômico, associar a cada cluster uma “probabilidade” de associação. Todavia, no método arbitramos uma constante (m) de fuzzyficação. Assim o seguinte script realiza o agrupamento fuzzy (por intermédio da função cmeans() da biblioteca e1071) para vários valores m assim é possível observar o comportamento (decaimento) do erro associado a clusterização e também a disposição dos fatores em relação ao cluster.

Infelizmente (ou para alegria geral) não vou me atrever a discorrer da teoria envolvida nessa análise e um estudo específico deixo a cargo dos maiores interessados.

Todo esse mise en cene, me soa um pouco genérico demais, eu sei (–Me desculpem!). Mas acredito que o script possa ser útil de alguma forma.

Aproveito o ensejo para agradecer as sugestões do Prof Julio da Motta Singer ao trabalho e o incentivo a esse post. Também deixo um alô ao colegas que participaram, de uma forma ou de outra!

Vamos ao trabalho, primeiro instalando (para quem ainda não o fez) e carregando a biblioteca e o conjunto de dados!

## install.packages('e1071')
library(e1071) # carrega pacote com implementação 'c-means'

dados <- transform(read.table("http://dl.dropbox.com/u/38195533/fac_exemplo.txt", # arquivo
                              header = TRUE,
                              sep = '\t',
                              na.string = '.',
                              dec = ','),
                   fatorA = factor(fatorA),
                   fatorB = factor(fatorB))

str(dados)

ajuste <- lm(resp ~ fatorA + fatorB, data = dados)
anova(ajuste) # ANOVA

matriz <- tapply(dados$resp, list(dados$fatorB, dados$fatorA), mean)

A primeira etapa foi realizada. Não posso passar ou próximo trecho sem deixa-los a par do quem está por vir! O agrupamento que “tentamos” é com relação ao um padrão de comportamento do fator sob análise, em relação ao outro fator, assim constituímos ideótipos e assim procedemos. O modo como constituimos esses ideótipos é aplicado unicamente a natureza desse conjunto de dados, mas outros ideótipos, ao seu critério podem ser usados, fique a vontade!

Antes de passarmos agora a análise em si manipulamos os dados e constituímos os centróides:

FinWil <- apply(matriz, 2, mean) - mean(matriz) # Finlay & Wilkinson

altos <- which(FinWil >= 0) # ambientes acima da media
baixo <- which(FinWil < 0) # ambientes abaixo da media

medias <- cbind(apply(matriz[, altos], 1, mean),
                apply(matriz[, baixo], 1, mean))

##--------------------------------------------------------------------
## IDEOTIPOS
## I   - maximo em todos 
## II  - maximo nos altos minimo nos baixos
## III - minimo nos altos maximo nos baixos
## IV  - minimo em todos 
## ADICIONAIS:
## V   - media em todos
## VI  - maximo nos altos e media nos baixos
## VII  - media nos altos maximo nos baixos

I <- c(max(matriz[, altos]), max(matriz[, baixo]))
II <- c(max(matriz[, altos]), min(matriz[, baixo]))
III <- c(min(matriz[, altos]), max(matriz[, baixo]))
IV <- c(min(matriz[, altos]), min(matriz[, baixo]))
V <- c(mean(matriz[, altos]), mean(matriz[, baixo]))
VI <- c(max(matriz[, altos]), mean(matriz[, baixo]))
VII <- c(mean(matriz[, altos]), max(matriz[, baixo]))

centroides <- rbind(I, II, III, IV, V, VI, VII) # centróides

Resta agora realizar a análise em si, veja que o laço for() repete a análise para vários valores da constante m e de cada uma, armazena as “probabilidades” de associação a cada centróide e o erro do ajustamento!


ms <- seq(1.2, 10, by = .1)
lista <- array(NA, dim = c(nlevels(dados$fatorB), 
                           nrow(centroides), 
                           length(ms)))
erros <- matrix(NA, nrow = length(ms), ncol = 1)
dimnames(lista) <- list(levels(dados$TRAT), 
                        c('I', 'II', 'III', 'IV', 'V', 'VI', 'VII'), 
                        ms)

for(i in 1:length(ms))
  {
    j <- ms[i]
    ## função de cluster 'fuzzy c-means'
    grupamentos <- cmeans(medias, # matriz dos dados
                          centroides, # centróides dos clusters
                          verbose = FALSE, # não imprime histórico
                          dist = 'euclidean', # distância empregada
                          m = j, # constante de desfuzzyficação
                          method = 'cmeans') # método
    lista[,,i] <- grupamentos$membership
    erros[i,] <- grupamentos$withinerror
  }

Pronto… feito, quer dizer, quase! Realizamos portanto a análise para os vários valores de m resta agora “escolher” um deles e partir para o resto! Até a próxima…

PCA de três dimensões

Olá a todos. O objetivo desse post é mostrar como fazer a análise de componentes principais em três dimensões. Mas esse post tem o destaque não ter sido escrito sozinho… Sem a participação da Msc Thalita (ou “Hermínia”, ninguém vai entender!) não teria post. Agradeçam a ela, mas reclamem comigo! Pois bem, Vamos usar como exemplo dados de um estudo com fatorial duplo (4 \times 4), em que foram tomadas quatro características (dados disponíveis para teste!).

Vou aproveitar e deixar aqui também um alô para o Dsc Guilherme pela ajuda nos detalhes do script, Valeu!

Um pouco de teoria:

Experimentos fatoriais (n \times n) caracterizam-se pela apresentação dos dados em tabelas de duas entradas (matriz), sendo que cada casela da tabela contém a resposta média de cada combinação dos fatores sob análise. Todas possibilidades desses experimentos estão fora do escopo desse post, por isso vamos um pouco além…

Acrescentando-se um outro fator, ou ainda outra característica (n \times n \times n, por exemplo), com 3 interações duplas e uma tripla, isso fica um pouco mais complexo. Neste caso, os dados são organizados em “CUBO” (arranjo de três entradas).

Nesta situação não é possível aplicar a análise de componentes principais usual. Uma alternativa é a utilização da decomposição PCAn, proposta por Tuckey em 1964.

Portanto nesse post pretendemos (afinal esse foi um trabalho em equipe) mostrar como realizar a decomposição de um arranjo de três entradas (CUBO). E mais que isso, mostrar uma função que realiza a escolha entre os vários modelos possíveis que variam conforme a dimensão dos dados.

Vamos ao código… mas calma, primeiro preparamos os dados!

Veja que usamos um pacote (PTAk, disponível no CRAN). Caso você não o tenha, proceda com install.packages(‘PTAk’)!

require(PTAk) # pacote analises pscicrometricas

## le e transforma os dados em fatores
dados <- transform(read.table("http://dl.dropbox.com/u/38195533/dados.txt")
                              na.string = '.',
                              header = TRUE,
                              dec = ',',
                              sep = '\t'),
                   elemento = factor(elemento),
                   residuo = factor(residuo),
                   metodo = factor(metodo))

CUBO.inicial <- with(dados,
                     tapply(conc, list(residuo, elemento, metodo),
                            mean, na.rm = TRUE))
## padronização dos dados
CUBO.standard <- apply(CUBO.inicial, 2, scale, center = FALSE)

## reconstrucao da planilha
estrutura <- expand.grid(elemento = factor(1:4),
                         residuo = factor(1:4),
                         metodo = factor(1:4))

plan <- data.frame(estrutura, conc = c(CUBO.standard)) ## dados

ajuste <- aov(conc ~ elemento * residuo * metodo - elemento:residuo:metodo,
              data = plan)
erros <- ajuste$residuals

## CUBO
Z <- tapply(erros, list(plan$residuo, plan$elemento, plan$metodo), mean) 

Nesse trecho são lidos os dados, na forma usual de uma planilha de dados (uma coluna para cada variável). Pela natureza dos dados, estes são padronizados (N(\mu = 0, \sigma^2 = 1)). Depois a panilha de dados é “reconstruída” e feito o ajuste do modelo linear respectivo (sem a interação tripla — objeto do nosso estudo!). Do resultado da análise extraímos os resíduos e construímos o CUBO.

Quase tudo pronto! O próximo trecho é a função que efetivamente faz os vários ajustes possíveis (isso é função da dimensão dos dados).

otimim.PCAn <- function(cubo) {
  dimensoes <- dim(Z) - 1
  diagnostico <- vector('list', length = prod(dimensoes - 1))
  ajustes <- vector('list', length = prod(dimensoes - 1))
  contador <- 0
  for(i in 2:dimensoes[ 1 ])
    {
      for(j in 2:dimensoes[ 2 ])
        {
          for(k in 2:dimensoes[ 3 ])
            {
              modelo <- c(i, j, k)
              tamanho <- sum(modelo)
              ## decomposicao Tucker -- PTAk
              tucker <- PCAn(Z, # array
                             dim = modelo, # dimensoes
                             test = 1E-12,
                             Maxiter = 1000,
                             smoothing = FALSE,
                             verbose = FALSE,
                             file = NULL)
              ## diagnostico
              porcentagem <- tucker[[ 3 ]]$pct # % explicada
              contador <- contador + 1
              ajustes[[ contador ]] <- tucker
              diagnostico[[ contador ]] <- c(modelo, tamanho, porcentagem)
            }
        }
    }
  ## diagnosticos dos ajustes
  tab.ajustes <- data.frame(do.call(rbind, diagnostico))
  tab.ajustes <- tab.ajustes[order(-tab.ajustes[, 4], -tab.ajustes[, 5]), ]
  por.total <- split(tab.ajustes, tab.ajustes[, 4])
  melhores <- vector('list', length = length(por.total))
  for(i in 1:length(melhores))
    {
      melhores[[i]] <- por.total[[i]][1, ]
    }
  ## melhores ajustes por dimensao
  tab.melhores <- do.call(rbind, melhores)
  tab.melhores <- tab.melhores[order(-tab.melhores[, 4]), ]
  ## melhor modelo
  melhor <- tab.melhores[(which((tab.melhores[,5]-
                          c(tab.melhores[,5][-1],0))>=5)[1]+1),]
  ## resposta -- melhor ajuste
  return(list(escolhido = melhor,
              resposta = ajustes[[ as.numeric(rownames(melhor)) ]]))
}

UFA!

Feita a função agora é só usa-lá (simples não?). Um uso possível dessa decomposição é representar e interpretar os dados graficamente em um biplot. Nessa condição de arranjos de três entradas são possíveis diferentes formas de gráficos biplots. Mas o objetivo do nosso post não é apresentá-los.

O trecho final aplica a função recém-construída (e retorna um controle de tempo de cada ajuste). Mostramos também como extrair as matrizes (A, B e C) e o arranjo G do melhor ajuste.

As matrizes de resposta e o arranjo G é que são úteis na construção dos gráficos. Aproveitem!

ajuste.PCA <- otimim.PCAn(Z)

A <- ajuste.PCA$resposta[[1]]$v
B <- ajuste.PCA$resposta[[2]]$v
C <- ajuste.PCA$resposta[[3]]$v

G <- ajuste.PCA$resposta[[3]]$coremat

summary(ajuste) # ANOVA do modelo
summary.PTAk(ajuste.PCA$resposta) # quanto cada componente explica

Voi lá…, Até mais pessoal!

Paralelização de Processos II

Dando continuidade ao que foi apresentado no post sobre paralelização de processos e motivado pelo comentário do colega J. Franky e fundamentalmente com a ajuda de Benilton Carvalho (aos quais deixo meu agradecimento), nesse “fast-post” vou mostrar uma outra implementação em paralelo.

O cenário é exatamente o mesmo do apresentado no post anterior (veja para entender). Agora, entretanto o processamento paralelo se dará em dois “níveis”, a saber: As nsmc simulações correm como um fork, ou seja, cópias do mesmo processo. Elas são implementadas como uma função e “correm” com a função mclapply(), o argumento mc.set.seed = TRUE garante que cada uma das nsmc terá uma semente geradora de números diferente.

Dentro de cada simualação temos ainda dois processos independentes dentro de cada ciclo (as duas estratégias de amostragem). Para um caso em que cada uma dessas estratégias demore um tempo grande, ao invés de esperar que a primeira termine para então começar a segunda, usamos uma versão de thread, com as funções parallel() e collect() e paralelizamos as duas estratégias.

O argumento mc.set.seed = TRUE tem o mesmo propósito que anteriormente, o argmento name = ‘foo’ é um grande facilitador… Vocês vão entender! Ao passar uma tarefa (a estratégia de amostragem) pela função parallel() é mesmo que dar um comando no terminal com um ‘nohup &’, ou seja, o computador executa a tarefa mas não “bloqueia” tal terminal para outro comando.

Depois de paralelizados, os processos são coletados pela função collect() na estrutura de uma lista, temos como argumentos da função (que não estão sendo usados, veja a documentação) o wait, ele é muito útil quando queremos que se espera ‘X’ tempo pelo fim das tarefas, após esse tempo, o processo continua.

A vantagem do uso do argumento name = ‘foo’ é que por default os nomes na estrutura da lista do collect recebem o número PID do processo, ao atribuir um nome fica muito mais fácil distribuir o que foi paralelizado.

Espero que tenham gostado… Lembrando que este é um exemplo puramente didático, que provavelmente não se aplica ao uso desses procedimentos. Segue portando código da terceira implementação, agora com o uso de paralelização em dois níveis.

#-----------------------------------------------------------------------------
## PARAMETROS DA SIMULAÇÃO -- usar os mesmo para as duas situações
n <- 5000 # numero de individuos
ciclos <- 30 # numero de ciclos
nsmc <- 100 # numero de simulacoes

## SEGUNDO CASO -- paralelizada
#-----------------------------------------------------------------------------
## IDEM (realiza um simulacao de 'nsmc')
simula.ii <- function(x) { # inicio da funcao -- sem argumentos
  ## mesmos comentarios do caso acima
  resultados <- matrix(NA, ciclos, 4) 
  p0 <- sample(c(0, 1, 2), n, replace = TRUE, prob = c(.25, .5, .25))
  resultados[1, ] <- rep(c(mean(p0), var(p0)), times = 2)
  estrA <- sample(p0, n/5, replace = FALSE)
  estrB <- p0[seq(1, 5000, 5)]
  resultados[2, ] <- c(mean(estrA), var(estrA), mean(estrB), var(estrB))
  for(k in 3:ciclos) {
    ## pA e pB paralelizados em 'thread'
    pA <- parallel(tipoA(estrA), name = 'A', mc.set.seed = TRUE)
    pB <- parallel(tipoB(estrB), name = 'B', mc.set.seed = TRUE)
    estrG <- collect(list(pA, pB)) # coleta dos processos paralelizados
    estrA <- estrG$A; estrB <- estrG$B # redistribui os processos
    resultados[k, ] <- c(mean(estrA), var(estrA), mean(estrB), var(estrB))
  }
  return(resultados) # retorna uma simulacao
}

tempoC <- system.time({ # armazena o tempo de processamento
saida3 <- mclapply(1:nsmc, # numero de nsmc
# aplica a funcao -- faz a simulacao
                   simula.ii,
                   mc.preschedule = FALSE,
                   # 'expande' quantos processos forem possiveis
                   mc.set.seed = TRUE,
                   # uma semente para cada processo
                   mc.cores = getOption('cores'))
                   # usa quantos processadores forem possiveis
})

## TEMPO
tempoC[3]

Até a próxima!

Paralelização de processos

Hoje em dia os computadores estão cada vez mais rápidos e via de regra, nossos processsadores tem mais de um núcleo de processamento. Isso é um “horizonte” grande que nos permite avançar num campo “baixo nível” da informática de gerenciar nossos processos para que eles sejam feitos em paralelo entre os vários núcleos do processador. Esse procedimento tem recomendações bem específicas e portanto não tem reflexos gerais em qualquer aplicação.

Estou tentando aprender sobre isso e ultimamente andei lendo e brincando com esses procedimentos que estão disponíveis no R. Basicamente existem três tipos de paralilazação. As de baixo nível, que envolvem um profundo conhecimento do hardware e assim como li, não são tão recomendadas para usuários finais. Outras são as paralizações explícitas em que o usuário é quem “manda” e por fim, já está disponível no R (exceto plataforma Windows) pacotes em que é possível implementar o processo de paralelização implícita, em que o usuário interage com ferramentas internas de “baixo nível” e o computador faz o resto… (brincadeira).

No primeiro tipo, se esse for seu propósito veja a documentação dos pacotes multicore (funções fork() e children()). No segundo caso procure sobre o pacote snow. E no último caso, ao qual vou dar mais atenção, veja sobre os pacotes doMC e multiocore especificamente nas funções mclapply(), parallel() e collect().

Como exemplo vou apresentar um caso bem tosco por sinal, mas acho que bastante ilustrativo. A primeira parte do código é uma implementação bem naive de uma simulação em que sobre uma população são retiradas dois tipos de amostra, para esses dois casos são em sequencia em um certo número de ciclos retiradas amostras. Esse processo por sua vez é repetido nsmc vezes (abreviatura de numero de simulações de Monte Carlo). Para as duas implementações (a naive e a paralela) é armazenado o tempo de processamento. Veja a diferença… Se você quiser ver a coisa acontecendo, abra o capô do seu computador, no terminal digite top, vai abrir o gerenciador de processos, e você vai conseguir ver no caso da implementação em paralelo, as várias cópias do seu processo correndo, diferentemente da implementação naive em que cada simulação só começa quando a anterios acaba.

Segue o código

#-----------------------------------------------------------------------------
## PARAMETROS DA SIMULAÇÃO -- usar os mesmo para as duas situações
n <- 5000 # numero de individuos
ciclos <- 30 # numero de ciclos
nsmc <- 100 # numero de simulacoes

## PRIMEIRO CASO -- "naive"
#-----------------------------------------------------------------------------
saida1 <- vector('list', length = nsmc) # saida da implementacao SIMPLES
tempoA <- system.time({ # armazena o tempo de processamento
   for(i in 1:nsmc) { # laço -- para cada simulacao
   resultados <- matrix(NA, ciclos, 4) # matriz dos resultados individuais
   # populacao inicial
   p0 <- sample(c(0, 1, 2), n, replace = TRUE, prob = c(.25, .5, .25))
   resultados[1, ] <- rep(c(mean(p0), var(p0)), times = 2) # resultado ciclo 0
   estrA <- sample(p0, n/5, replace = FALSE) # estrategia A
   estrB <- p0[seq(1, n, 5)] # estrategia B
   # resultados ciclo 1
   resultados[2, ] <- c(mean(estrA), var(estrA), mean(estrB), var(estrB))
   for(k in 3:ciclos) { # laço -- ciclos subsequentes
      ## populacao estrategia A ciclo n
      pA <- sample(estrA, n, replace = TRUE)
      estrA <- sample(p0, n/5, replace = FALSE)
      ## populacao estrategia B ciclo n
      pB <- sample(estrB, n, replace = TRUE)
      estrB <- pB[seq(1, n, 5)]
      # resultados ciclo n
      resultados[k, ] <- c(mean(estrA), var(estrA), mean(estrB), var(estrB))
   }
   saida1[[i]] <- resultados # armazena na saida
}
})

## SEGUNDO CASO -- paralelizada
#-----------------------------------------------------------------------------
## IDEM (realiza um simulacao de 'nsmc')
simula.i <- function(x) { # inicio da funcao -- sem argumentos
   ## mesmos comentarios do caso acima
   resultados <- matrix(NA, ciclos, 4)
   p0 <- sample(c(0, 1, 2), n, replace = TRUE, prob = c(.25, .5, .25))
   resultados[1, ] <- rep(c(mean(p0), var(p0)), times = 2)
   estrA <- sample(p0, n/5, replace = FALSE)
   estrB <- p0[seq(1, n, 5)]
   resultados[2, ] <- c(mean(estrA), var(estrA), mean(estrB), var(estrB))
   for(k in 3:ciclos) {
      pA <- sample(estrA, n, replace = TRUE)
      estrA <- sample(p0, n/5, replace = FALSE)
      pB <- sample(estrB, n, replace = TRUE)
      estrB <- pB[seq(1, n, 5)]
      resultados[k, ] <- c(mean(estrA), var(estrA), mean(estrB), var(estrB))
   }
   return(resultados) # retorna uma simulacao
}

tempoB <- system.time({ # armazena o tempo de processamento
saida2 <- mclapply(1:nsmc, # numero de nsmc
# aplica a funcao -- faz a simulacao
                   simula.i,
                   mc.preschedule = FALSE,
                   # 'expande' quantos processos forem possiveis
                   mc.set.seed = TRUE,
                   # uma semente para cada processo
                   mc.cores = getOption('cores'))
                   # usa quantos processadores forem possiveis
})

## TEMPOS
tempoA[3]
tempoB[3]

Vale lembrar que é preciso estar usando ou Linux, ou Mac e é preciso antes fazer o download e chamar o pacote “multicore”. Espero que tenham gostado e que com esse exemplo consigam realizar suas próprias paraleilazações.

Categorias:linux, notícias, testes

Os testes chi-quadrado

Gráfico de barras representando as frequências absolutas das classes de aspecto de agregado em função da profundidade de coleta. Valores dentro do gráfico são resultados do teste de homogeneidade.

Os testes \chi^2 (qui-quadrado) a que me refiro são:

  • Teste de aderência: testa a hipótese da amostra ser proveniente de uma distribuição de probabilidade definida em H_0. Com essa distribuição definida em H_0 são obtidos as frequências esperadas (E);
  • Teste de homogeneidade: testa a hipótese H_0 de duas ou mais amostras serem provenientes de uma mesma distribuição de probabilidades. Os valores esperados são obtidos pelo produto da linha marginal e tamanho das amostras;
  • Teste de independência: testa a hipótese H_0 de que a distribuição conjunta é o produto das distribuições marginais, o que só ocorre quando existe independência entre as variáveis aleatórias. No caso de duas variáveis aleatórias organizadas numa tabela de dupla entrada, os valores esperados são obtidos como produto dos valores marginais.

Nos testes chi-quadrado o que muda é só a hipótese envolvida no cálculo dos valores esperados. Para os três tipos de hipótese, a estatística do teste é

X^2 = \sum_{i=1}^{n} \displaystyle \frac{(O_i-E_i)^2}{E_i}

sendo que sob H_0 a variável aleatória X^2 \sim \chi^2_\nu em que \nu são os graus de liberdade.

Nesse post vou apresentar cada uma dos três tipos de teste de hipótese. Para uma melhor abordagem teórica do teste chi-quadrado consulte os livros de estatística básica como Estatítica Básica do Bussab e Morettin.

Embora alguns autores façam distinção entre o teste de homogeneidade e de independência, você vai perceber que são o mesmo teste para a mesma hipótese escrita de duas formas. Veja, se há independência entre as classificações então é esperado que os valores para a combinações sejam o produto das probabilidades marginais, pois P(A_i \cap B_i) = P(A_i)\cdot P(B_i) sob independência. Logo, se as probabilidades na linha marginal representam os as probabilidades nas linhas de cada amostra, então há homogeneidade. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# teste de aderência 1: frequencia de acidentes nos dias da semana
# hipótese H_0 é de as frequências são dadas por uma distribuição
# uniforme discreta com n=5, ou seja, p_i=5 para todo i={seg,ter,qua,qui,sex}
# dados do Bussab & Morettin - Estatística Básica - 6 edição, pg 404

Oi <- c(seg=32, ter=40, qua=20, qui=25, sex=33) # observados
Ei <- sum(Oi)*1/length(Oi)                      # esperados sob H_0
X2 <- sum((Oi-Ei)^2/Ei)                         # estatística do teste
nu <- length(Oi)-1                              # graus de liberdade
pchisq(X2, df=nu, lower.tail=FALSE)             # valor-p do teste

#-----------------------------------------------------------------------------
# usando a função chis.test()

chisq.test(Oi)

#-----------------------------------------------------------------------------
# teste de aderência 2: número de plantas por m² em uma floresta
# hipótese H_0 é de o número de plantas observados é Poisson(lambda)
# lambda precisa ser estimado e isso diminui um grau de liberdade

n <- 0:10                                       # número de plantas
Oi <- c(574,922,1172,917,609,324,150,64,19,4,0) # número observado
names(Oi) <- c(n[-11],">9")
lambda <- sum(Oi*n)/sum(Oi)                     # estimativa de lambda
pi <- dpois(n[-11], lambda=lambda)              # frequência sob H_0
pi <- c(pi, 1-sum(pi))
Ei <- sum(Oi)*pi                                # número esperado sob H_0
X2 <- sum((Oi-Ei)^2/Ei)                         # estatística do teste
nu <- length(Oi)-1-1                            # grau de liberdade
pchisq(X2, df=nu, lower.tail=FALSE)             # valor p do teste

#-----------------------------------------------------------------------------
# usando a função chis.test()

chisq.test(Oi, p=pi)

# os graus de liberdade não consideram a estimação de lambda
# prestar atenção quando usar a chisq.test() nestes casos

#-----------------------------------------------------------------------------
# teste de homogeneidade: testar se a distribuição do aspecto de agregados
# muda com a profundidade de amostragem
# os dados são contínuos e serão colocados em classes para aplicação do teste
# foram usadas as classes obtidas para construir um histograma

ag <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/agreg.txt",
                 header=TRUE, sep="\t")
str(ag)
ht <- hist(ag$aspecto)
classes <- ht$breaks                                # classes de aspecto
cla <- cut(ag$asp, classes)                         # atribuição às classes
Oi <- table(ag$prof, cla)                           # observados
Ei <- outer(rowSums(Oi), colSums(Oi), "*")/sum(Oi)  # esperados sob H_0
X2 <- sum((Oi-Ei)^2/Ei)                             # estatística do teste
nu <- prod(dim(Ei)-1)                               # graus de liberdade
pchisq(X2, df=nu, lower.tail=FALSE) -> P; P         # valor p do teste

#-----------------------------------------------------------------------------
# usando a função chis.test()

chisq.test(Oi)

#-----------------------------------------------------------------------------
# gráfico

#png("f016.png", w=500, h=300)
par(mar=c(4.1,6.1,2.1,2.1))
col <- c("green3","green4")
barplot(Oi, beside=TRUE, horiz=TRUE, las=1, col=col,
        xlab="Frequência absoluta")
mtext(side=2, text="Classe de aspecto do agregado", line=5)
legend("bottomright", legend=c("0 - 5 cm","5 - 20 cm"), fill=col, bty="n")
text(30, 3, substitute(italic(X)^2==x~~~~~~italic(valor-p)==P,
                       list(x=round(X2,4), P=round(P,4))))
#dev.off()

#-----------------------------------------------------------------------------
# teste de independência: testar se há independência na classificação
# quanto ao grau de ingestão de alcool (0, <7, >7 copos por semana) e o
# desenvolvimento de doença cardíaca (sim, não)

Oi <- matrix(c(146,106,29,750,590,292), byrow=TRUE, # observados
             2, 3, dimnames=list(c("sim","não"), c(0,"<7",">7")))
Ei <- outer(rowSums(Oi), colSums(Oi), "*")/sum(Oi)  # esperados sob H_0
X2 <- sum((Oi-Ei)^2/Ei)                             # estatística do teste
nu <- prod(dim(Ei)-1)                               # graus de liberdade
pchisq(X2, df=nu, lower.tail=FALSE)                 # valor p do teste

#-----------------------------------------------------------------------------
# usando a função chis.test()

chisq.test(Oi)

#-----------------------------------------------------------------------------