#==============================================================================
#  Comparação de Médias - Teste F - ANOVA
#==============================================================================
#
#  Descrição do Estudo: Volume Expiratório Forçado VEF - Três Grupos
#                        Livro do Marcelo Pagano Cap. 12
#
#===============================================================================
#
#  Leitura dos Dados
#
 ex01 <- read.table("dados_aula2.txt", head=T)
 ex01
# Caso o arquivo esteja em outro diretório deve-se colocar o caminho completo 
# deste diretório no argumento de read.table acima. 
#
# Inspecionando os dados e suas componentes. 
#
 is.data.frame(ex01)
 names(ex01)
 ex01$volume
 ex01$grupo
#
 is.factor(ex01$grupo)
 ex01$grupo<-factor(ex01$grupo)
 is.numeric(ex01$volume)
#
# Concluímos que o objeto ex01 é um data-frame com duas variáveis, sendo uma delas um 
# fator (a variável grupo) e a outra uma variável numérica. 
#
# Análise descritiva 
#
  summary(ex01)
  tapply(ex01$volume, ex01$grupo, length)
  tapply(ex01$volume, ex01$grupo, mean)
  tapply(ex01$volume, ex01$grupo, median)
  tapply(ex01$volume, ex01$grupo, sd)
  tapply(ex01$volume, ex01$grupo, max)-tapply(ex01$volume, ex01$grupo, min)
  tapply(ex01$volume, ex01$grupo, summary)
#
# Existe uma indicação inicial que o grupo 2 tem média/mediana superior aos
# outros dois.
# Media ~= mediana indicação de simetria de distribuições e mesmo dp/range é uma
# indicação de homocedasticidade.
#  
# Há um mecanismo no R de "anexar" objetos ao caminho de procura que permite 
# economizar um pouco de digitação. Veja os comandos abaixo e compare com o 
# comando anterior. 
#
  attach(ex01)
#
 tapply(volume, grupo, mean)
#
# Interessante não? Quando "anexamos" um objeto do tipo list ou data.frame 
# no caminho de procura com o comando attach() fazemos com que os componentes 
# deste objeto se tornem imediatamente disponíveis e portanto podemos, por 
# exemplo, digitar somente grupo ao invés de ex01$grupo. 
#
# No entanto, este comando deve ser utilizado com cuidado pois o uso abusivo
# vai colocando objeto em cima de objeto.
#
# Vamos prosseguir com a análise exploratória, obtendo algumas medidas e 
# gráficos. 
#
  plot(ex01)
#  plot(ex01, pch="x", col=2, cex=1.5) # muda a marca, cor e intensidade
  plot(ex01$grupo,ex01$volume)
#
  boxplot(volume ~ grupo) # mesmo gráfico anterior
#
#  Agora vamos fazer a análise de variância. Vamos "desanexar" o objeto com 
#  os dados (embora isto não seja necessário). 
  detach(ex01)
  ex01.av <- aov(volume ~ factor(grupo), data = ex01)
  ex01.av
#
  summary(ex01.av)
  anova(ex01.av)
#
# Portanto o objeto ex01.av guarda os resultados da análise. Vamos 
# inspecionar este  objeto mais cuidadosamente e fazer também uma análise 
# dos resultados e resíduos. 
#
  names(ex01.av)
  class(ex01.av)
  ex01.av$coef
#
  ex01.av$res
  residuals(ex01.av)
#
#  O valor-p=0,052 somente é válido se as suposições forem confirmadas.
#
#  Vamos verificar as suposições do modelo/teste F
#   1- homocedasticidade e 2- normalidade
#
#  A homogeneidade dos desvios-padrão pode ser verificada pelo 
#  Teste de Bartlett. Um teste alternativo,o de Levene,
#  (usualmente melhor que o de Bartlett) pode ser realizado no pacote
#  lawstat
#
  bartlett.test(volume, grupo)
  require(lawstat) 
  levene.test(volume, grupo)
#
#  Ambos testes mostram não haver evidência contra a hipótese de
#  homocedasticidade. Os resíduos também são utilizados para verificar
#  as suposições.
#
  plot(ex01.av) # pressione a tecla enter para mudar o gráfico
#
 par(mfrow=c(2,2))
 plot(ex01.av)
 par(mfrow=c(1,1))
#
 plot(ex01.av$fit, ex01.av$res, xlab="valores ajustados", ylab="resíduos")
 title("resíduos vs Preditos")
#
# o primeiro (e o terceiro) mostrando faixas de mesma dispersão atestam
# a adequação da homocedasticidade. E o segundo, mostrando os pontos em torno
# da reta, mostram a adequação da normalidade.
#
# Mais gráficos para verificar a distribuição dos resíduos
#
 names(anova(ex01.av))
 s2sq <- anova(ex01.av)$Mean[2]   # estimativa da variância
 res <- ex01.av$res        # extraindo resíduos
 respad <- (res/sqrt(s2))  # resíduos padronizados 
 boxplot(respad)
 title("Resíduos Padronizados" )
#
 hist(respad, main=NULL)
 title("Histograma dos resíduos padronizados")
#
 stem(respad)
 qqnorm(res,ylab="Residuos", main=NULL)
 qqline(res)
 title("Gráfico Normal de Probabilidade dos Resíduos")
#
# O teste de Shapiro-Wilks verifica a normalidade.
#
 shapiro.test(res)
#
# Como as suposições foram atendidads, nós dizemos que existem
# diferença entre os grupos ao nível de 10%. Portanto, necessitamos
# de comparações múltiplas para encontrar quais grupos são diferentes
# e quais são iguais. Vamos utilizar os métodos de Bonferroni e Tukey.
#
 pairwise.t.test(ex01$volume, ex01$grupo, p.adj = "bonf")
#
 ex01.tu <- TukeyHSD(ex01.av)
 plot(ex01.tu)
 ex01.tu
#
# Concluimos que existe diferença entre os grupos 1 (John Hopkins) e o 2(Rancho
# Los Amigos).
#======================================================================================
#  TÓPICOS ESPECIAIS: Outros Métodos de Análise
#============================================================================
# 1- Teste Não-Paramétrico de Kruskal-Wallis
#
  kruskal.test(volume ~ grupo, data = ex01) 
#=============================================================================
# 2- Métodos de Comparação Múltiplas
#
  all.pairs <- function(r)
  list(first = rep(1:r,rep(r,r))[lower.tri(diag(r))],
       second = rep(1:r, r)[lower.tri(diag(r))])
  fitted<- tapply(ex01$volume, ex01$grupo, mean)
  nis<- tapply(ex01$volume,ex01$grupo,length)
  df<-ex01.av$df.residual
  MSE=sum(ex01.av$residuals^2)/df
#
#  2.1 - Tukey
 tukeyCI <- function(fitted, nis, df, MSE, conf.level=.95){
  ## fitted is a sequence of means
  ## nis is a corresponding sequence of sample sizes for each mean
  ## df is the residual df from the ANOVA table
  ## MSE = mean squared error from the ANOVA table
  ## conf.level is the family-wise confidence level, defaults to .95
  r <- length(fitted)
  pairs <- all.pairs(r)
  diffs <- fitted[pairs$first] - fitted[pairs$second]
  df <- sum(nis) - r
  T <- qtukey(conf.level, r, df)/sqrt(2)
  hwidths <-  T*sqrt(MSE*(1/nis[pairs$first] + 1/nis[pairs$second]))
  val <- cbind(diffs - hwidths, diffs, diffs + hwidths)
  dimnames(val) <- list(paste("mu",pairs$first," - mu", pairs$second,
  sep=""), c("Lower", "Diff","Upper"))
  val
}
 tukeyCI(fitted, nis, df, MSE, conf=.90)
#
# 2.2 Scheffé
 scheffeCI <- function(fitted, nis, df, MSE, conf.level=.95){
  ## fitted is a sequence of means
  ## nis is a corresponding sequence of sample sizes for each mean
  ## df is the residual df from the ANOVA table
  ## MSE = mean squared error from the ANOVA table
  ## conf.level is the family-wise confidence level, defaults to .95
  r <- length(fitted)
  pairs <- all.pairs(r)
  diffs <- fitted[pairs$first] - fitted[pairs$second]
  T <- sqrt((r-1)*qf(conf.level,r-1,df))
  hwidths <-  T*sqrt(MSE*(1/nis[pairs$first] + 1/nis[pairs$second]))
  val <- cbind(diffs - hwidths, diffs, diffs + hwidths)
  dimnames(val) <- list(paste("mu",pairs$first," - mu", pairs$second,
  sep=""), c("Lower", "Diff","Upper"))
  val
}
 scheffeCI(fitted, nis, df, MSE, conf=.90)
#
# 2.3 Bonferroni
 bonferroniCI <- function(fitted, nis, df, MSE, conf.level=.90){
  ## fitted is a sequence of means
  ## nis is a corresponding sequence of sample sizes for each mean
  ## df is the residual df from the ANOVA table
  ## MSE = mean squared error from the ANOVA table
  ## conf.level is the family-wise confidence level, defaults to .95
  r <- length(fitted)
  pairs <- all.pairs(r)
  diffs <- fitted[pairs$first] - fitted[pairs$second]
  T <- qt(1-(1-conf.level)/(r*2),df)
  hwidths <-  T*sqrt(MSE*(1/nis[pairs$first] + 1/nis[pairs$second]))
  val <- cbind(diffs - hwidths, diffs, diffs + hwidths)
  dimnames(val) <- list(paste("mu",pairs$first," - mu", pairs$second,
  sep=""), c("Lower", "Diff","Upper"))
  val
}
  bonferroniCI(fitted, nis, df, MSE, conf=.90)
#==========================================================================
#
#   3-   Teste via bootstrap
#       Teste de Permutação (Sem Reposição)
#
# amostra observada: 21 + 16 + 23
 b <- 20000 # número de amostras bootstrap
 n<-21+16+23
 chi.boot<-NULL
#
 for (j in 1:b) { # começo do for do bootstrap
        x.boot<-NULL
	i. <- sample(ex01$grupo,replace=F)
	       chi.boot[j] <- anova( aov(ex01$volume ~ factor(i.)))$F[1] 
 } # fim do for do bootstrap
# hist(chi.boot)
 boot.p.value<-length(chi.boot[chi.boot>=anova(aov(ex01$volume ~ factor(
 ex01$grupo)))$F[1] ])/b # 
 boot.p.value
#
# fim do for do bootstrap
  x<-seq(0,8,.05)
 hist(chi.boot,freq=F)
 lines(x,df(x,2,57),lty=1)
#===========================================================================
#  4- Modelando a dispersão 
#
  mod0<-lm(volume~grupo,data=ex01)
  summary(mod0)
#
  require(gamlss)
  mod1<-gamlss(volume ~ grupo,sigma.formula= ~grupo)
  summary(mod1)             
#
# Confirmando que não é necessário modelar a dispersão
#   
  wp(mod1)
  plot(mod1)
#============================================================================
#