#=============================================================================
#    Modelagem e análise de dados experimentais com o programa computacional R
#           http://www.leg.ufpr.br/doku.php/pessoais:walmes:cursoragrarias2012
#                                                               Walmes Zeviani
# Aula 8 - 01/11/2012
# Análise de experimento com respostas de proporção (binomial)
#=============================================================================

#------------------------------------------------------------------------------------------
# carregando os dados

mud <- read.table("http://www.leg.ufpr.br/~walmes/data/mudas.txt", sep=";", header=TRUE)
str(mud)
mud$bloc <- factor(mud$bloc)
levels(mud$trat)

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

require(lattice)
xyplot(viab~trat|bloc, data=mud)
xyplot(viab~trat|bloc, data=mud, type=c("p","a"))

#------------------------------------------------------------------------------------------
# qual foi o tamanho dessa binomial

xtabs(~trat+bloc, mud)     # é o número de plantas por combinação
xtabs(viab~trat+bloc, mud) # o número que sobreviveu

#------------------------------------------------------------------------------------------
# fazendo outro data.frame com o número de viáveis

require(plyr)

mud2 <- ddply(mud, .(bloc,trat), summarise, nvia=sum(viab))
str(mud2)
xyplot(nvia~trat, groups=bloc, data=mud2, type="o")

#------------------------------------------------------------------------------------------
# análise usando glm()

g0 <- glm(cbind(nvia, 20-nvia)~bloc+trat, data=mud2, family=binomial)
par(mfrow=c(2,2)); plot(g0); layout(1) # análise dos resíduos
anova(g0, test="Chisq") # quadro de análise de DEVIANCE (variância = desviancia)
                        # checar diferença entre Resir. Df e Resid. Dev

pchisq(deviance(g0), df=df.residual(g0), lower.tail=FALSE) # não tem superdispersão (ótimo)

#------------------------------------------------------------------------------------------
# obtendo o intervalo de confiança para \mu_\alpha_i, reajustar modelo com contr.sum
# remove intercepto e colocar os efeitos de bloco sob a retrição soma 0

g1 <- glm(cbind(nvia,20-nvia)~-1+trat+bloc, data=mud2, family=binomial,
          contrast=list(bloc=contr.sum))
summary(g1) # estimativas não são probabilidades, p=exp(e)/(1+exp(e))

ef <- coef(g1)[1:16] # isso são odds ratio
exp(ef)/(1+exp(ef))  # isso são probabilidades

x <- 0:2
p <- exp(x)/(1+exp(x)); p
exp(diff(x))/(1+exp(diff(x)))

#------------------------------------------------------------------------------------------
# obtendo intervalos de confiança

IC.nu <- confint(g1) # perfilhamento, IC assimétrico
IC.nu # intervalo de confiança do parâmetro, temos que passar o inverso do link
IC.nu <- cbind(fit=coef(g1), IC.nu)
IC.mu <- 1/(1+exp(-IC.nu))
IC.mu # intervalo de confiança para proporção de viáveis
IC.mu <- cbind(trat=rownames(IC.mu), as.data.frame(IC.mu))
names(IC.mu)[3:4] <- c("L","U")
IC.mu <- IC.mu[grep("trat", IC.mu$trat),]
str(IC.mu)

#------------------------------------------------------------------------------------------
# gráfico, note os intervalos de confiança assimétricos a medida que afastam de 0.5
# quando mais perto de 0.5 mais amplo também

require(latticeExtra)

segplot(reorder(trat, fit)~L+U, center=fit, data=IC.mu, draw.bands=FALSE,
        segments.fun=panel.arrows, ends="both", angle=90, length=0.1, col=1)

#------------------------------------------------------------------------------------------
# podemos comparar as médias (no caso as proporções, média=n*p) dos substratos

require(multcomp)

glht0 <- summary(glht(g1, linfct=mcp(trat="Tukey")), # contrastes de Tukey
                 test=adjusted(type="single-step"))  # tipo de ajuste de p-valor
glht0
let <- cld(glht0) # se for para associar as letras, nem sempre recomendável
let <- let$mcletters$Letters
let

#------------------------------------------------------------------------------------------
# apromirando o gráfico

IC.mu <- IC.mu[order(IC.mu$fit),]       # ordena pelas médias
IC.mu$let <- sort(let, decreasing=TRUE) # poe coluna com as letras
IC.mu$tratlet <- paste(gsub("^trat","", as.character(IC.mu$trat)), # rótulo com letras
                       " (", IC.mu$let,")", sep="")

segplot(reorder(tratlet, fit)~L+U, center=fit, data=IC.mu, draw.bands=FALSE,
        segments.fun=panel.arrows, ends="both", angle=90, length=0.1, col=1)

#------------------------------------------------------------------------------------------
# como ficaria se eu assumisse normalidade?

m0 <- lm(nvia~bloc+trat, data=mud2)
par(mfrow=c(2,2)); plot(m0); layout(1)

# vamos aos fatos:
# * o size é grande (20), quanto maior, mais p=x/n fica contínuo;
# * o p é próximo de 0.5, ficou entre 0.2 e 0.8;
# * essas coisas combinadas fazem com que a normal aproxime bem a binomial.

#------------------------------------------------------------------------------------------
# dados de número de sementes viáveis de soja

soja <- read.table("http://www.leg.ufpr.br/~walmes/data/soja.txt", header=TRUE, dec=",")
soja <- transform(soja, k=factor(potassio), a=factor(agua))
names(soja)[1:3] <- c("K","A","bloc")
str(soja)

#------------------------------------------------------------------------------------------
# ver

xyplot(nv/(nv+nvi)~K|a, soja, type=c("p","a"))

#------------------------------------------------------------------------------------------
# ajuste modelo de caselas aos dados assumindo distribuição binomial (link=logit)

g0 <- glm(cbind(nv, nvi)~bloc+k*a, data=soja, family=binomial)
par(mfrow=c(2,2)); plot(g0); layout(1) # ok

#------------------------------------------------------------------------------------------
# quadro de análise de deviance, faz a vez da anova

anova(g0, test="Chisq")
pchisq(deviance(g0), df=df.residual(g0), lower.tail=FALSE) # tem superdispersão

#------------------------------------------------------------------------------------------
# usando polinômios para efeito do K

g0 <- glm(cbind(nv, nvi)~bloc+k*a, data=soja, family=quasibinomial)
g1 <- glm(cbind(nv, nvi)~bloc+a*poly(K,2), data=soja, family=quasibinomial)
g1 <- glm(cbind(nv, nvi)~bloc+a*(K+I(K^2)), data=soja, family=quasibinomial)
par(mfrow=c(2,2)); plot(g1); layout(1) # ok
anova(g0, g1, test="F") # abandono de termos de maior ordem não significativo
anova(g1, test="F")
summary(g1)

#------------------------------------------------------------------------------------------
# declarando modelo reduzido

g2 <- glm(cbind(nv, nvi)~bloc+a*K+I(K^2), data=soja, family=quasibinomial)
par(mfrow=c(2,2)); plot(g2); layout(1) # ok
anova(g0, g2, test="F") # abandono de termos de maior ordem não significativo
anova(g2, test="F")
summary(g2)

#------------------------------------------------------------------------------------------
# não perca seu tempo interpretando coeficientes, vamos ver o gráfico

pred <- expand.grid(bloc="I", a=levels(soja$a), K=seq(0,180,by=1))
pred$rengrao <- predict(g2, pred, type="response")

xyplot(rengrao~K, groups=a, data=pred, type="l", auto.key=TRUE,
       ylab="Probabilidade de viabilidade da vagem")

#------------------------------------------------------------------------------------------
# se eu fizesse o ajuste assumindo normal?

soja <- transform(soja, proporc=nv/(nv+nvi))

m0 <- lm(proporc~bloc+a*k, data=soja)
par(mfrow=c(2,2)); plot(m0); layout(1) # fuga!

require(MASS)

boxcox(m0, seq(0,10,l=30)) # que? elevar à 7ª?

#------------------------------------------------------------------------------------------
# e se eu ignorar os pressupostos?

m1 <- lm(proporc~bloc+a*K+I(K^2), data=soja)

pred <- expand.grid(bloc="I", a=levels(soja$a), K=seq(0,180,by=1))
pred$rengrao <- predict(m1, pred)

xyplot(rengrao~K, groups=a, data=pred, type="l", auto.key=TRUE,
       ylab="Probabilidade de viabilidade da vagem")

# epa! probabilidade maior que 1?

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