#------------------------------
# Página 355: Danos cerebrais - 
# Exemplo usado nos slides    -
#------------------------------
rm(list=ls())

# setwd('C:/Aulas/Ano2013/Patologia/R/Fatorial/')
              
#-------------------
# Entrada de dados -
#-------------------

(fis = factor(rep(c('I','II','III','IV','V','VI'), each=12)))
(psiq = factor(rep(LETTERS[1:4], 18)))
(Trat = factor(interaction(fis,psiq)))
							
(meses = c(11.0,  9.4, 12.5, 13.2,
            9.6,  9.6, 11.5, 13.2,
	         10.8,  9.6, 10.5, 13.5,
	         10.5, 10.8, 10.5, 15.0,
    	     11.5, 10.5, 11.8, 14.6,
	         12.0, 10.5, 11.5, 14.0,
	         12.0, 11.5, 11.8, 12.8,
	         11.5, 11.5, 11.8, 13.7,
	         11.8, 12.3, 12.3, 13.1,
	         11.5,  9.4, 13.7, 14.0,
	         11.8,  9.1, 13.5, 15.0,
	         10.5, 10.8, 12.5, 14.0,
	         11.0, 11.2, 14.4, 13.0,
	         11.2, 11.8, 14.2, 14.2,
	         10.0, 10.2, 13.5, 13.7,
	         11.2, 10.8, 11.5, 11.8,
	         10.8, 11.5, 10.2, 12.8,
	         11.8, 10.2, 11.5, 12.0))
							
(danos = data.frame(Trat, fis, psiq, meses))
attach(danos)
names(danos)
						
str(danos)

#----------------------------------------------
# Criando a tabela de dados para uso no Latex -
#----------------------------------------------
# d = matrix(meses, 18, 4, byrow=T)
# require(xtable)
# xtable(d)							
# addmargins(d)

#---------------------------
# Estatísticas descritivas -
#---------------------------
summary(meses)

(n = tapply(meses, list(fis, psiq), length))
(medias = tapply(meses, list(fis, psiq), mean))
	min(medias) ; max(medias)
(variâncias = tapply(meses, list(fis, psiq), var))
	min(variâncias) ; max(variâncias)
(desvios = tapply(meses, list(fis, psiq), sd))
(cv = desvios / medias * 100)
	min(cv) ; max(cv)

(medias.fis = tapply(meses, fis, mean))
(medias.psiq = tapply(meses, psiq, mean))
(medias.Trat = tapply(meses, Trat, mean))


(soma = tapply(meses, list(fis, psiq), sum))
addmargins(soma


#-----------------------
# Gráficos descritivos -
#-----------------------
par(mai=c(1, 1, .2, .2))
plot.default(Trat, meses, las=1, xlab='Tratamentos', 
                 ylab='', col='blue', bty='l', xaxt='n')
mtext('Tempos médios  (meses)', side=2, line=3.5)
axis(1, at=c(1:24))
points(medias.Trat, pch="*", col='red', cex=1.5)

par(mai=c(1,1,.2,.2))
plot.default(fis, meses, las=1, xlab='Terapia Física', 
                 ylab='', col='blue', bty='l')
mtext('Tempos médios  (meses)', side=2, line=3.5)
points(medias.fis, pch="*", col='red', cex=1.5)

par(mai=c(1,1,.2,.2))
plot.default(psiq, meses, las=1, xlab='Tratamento Psiquiátrico', 
                 ylab='', col='blue', bty='l')
mtext('Tempos médios  (meses)', side=2, line=3.5)
points(medias.psiq, pch="*", col='red', cex=1.5)

par(mai=c(1,1,.2,.2))
boxplot(meses ~ Trat, names=1:24,
        ylab="Tempos médios  (meses)", xlab="Tratamentos",
        las=1, col='LightYellow')
points(medias.Trat, pch="+", col=2, cex=1.5)

par(mai=c(1,1,.2,.2))
boxplot(meses ~ psiq, names=1:4,
        ylab="Tempos médios  (meses)", xlab="Tratamento Psiquiátrico",
        las=1, col='LightYellow')
points(medias.psiq, pch="+", col=2, cex=1.5)

par(mai=c(1,1,.2,.2))
boxplot(meses ~ fis, names=1:6,
        ylab="Tempos médios  (meses)", xlab="Terapia Física",
        las=1, col='LightYellow')
points(medias.fis, pch="+", col=2, cex=1.5)


#-------------------------------------------------
# Visualização de todas as médias do experimento -
#           em relação à média geral             -
#-------------------------------------------------
par(mai=c(1,1,.2,.2))
plot.design(danos, xlab='Fatores', ylab="", las=1,
    	      bty='l', col='blue')
mtext('Tempos médios  (meses)', side=2, line=3.5)


#--------------------------------------------------------------
# Em experimentos fatoriais é importante verificar se existe  -
# interação entre os fatores. Inicialmente vamos fazer isto   -
# graficamente e mais a frente faremos um teste formal para   -
# presença de interação. Os comandos a seguir são usados para -
# produzir os gráficos.                                       -
#--------------------------------------------------------------
par(mai=c(1, 1, .2, .2))
interaction.plot(fis, psiq, meses, las=1, xlab='Terapia Física',
         ylab='Tempos médios (meses)', col=c('red','blue','green','brown'), bty='l',
         trace.label = deparse(substitute(psqi)), lwd=2.5)

interaction.plot(psiq, fis, meses, las=1, xlab='Tratamento Psiquiátrico',
         ylab='Tempos médios (meses)', col=c('red','blue','green','brown','magenta','black'),
         bty='l', trace.label = deparse(substitute(fis)), lwd=2.5)

#----------------------------------------------------
# Para verificarmos se as variâncias são homogêneas -
#----------------------------------------------------
anava.danos = aov(meses~Trat)
summary(anava.danos)
							
shapiro.test(anava.danos$res)
bartlett.test(meses~Trat)
							
anava.d.fat = aov(meses ~ fis + psiq + fis*psiq, data=danos)
summary(anava.d.fat)
							
require(ExpDes.pt)
fat2.dic(fis, psiq, meses, quali=c(TRUE,TRUE), fac.names=c('fis','psiq'))