### R code from vignette source 'Examen.Rnw'
### Encoding: UTF-8

###################################################
### code chunk number 1: Examen.Rnw:37-38
###################################################
  if(exists(".orig.enc")) options(encoding = .orig.enc)


###################################################
### code chunk number 2: Examen.Rnw:212-242
###################################################
library(nlme) 
library(bbmle)

#Generando la base de datos: duraznos
set.seed(1234) 

nobs <- 96 
ntrat <- nobs/3 
fert <- gl(3, ntrat, labels=c("nf", "org", "inorg"))
pnf <- rnorm(ntrat, 37, sd= 15) 
porg <- rnorm(ntrat, 41, sd=10) 
pinorg <- rnorm(ntrat, 44, sd= 7)
prod <- c(pnf, porg, pinorg)

duraznos <- data.frame(prod, fert)

boxplot(prod~fert, data=duraznos, ylab= "Biomasa de los frutos (kg)", 
xlab= "Tratamientos") 

m1 <- lm(prod~fert, data=duraznos) 

anova(m1)

par(mfrow=c(2,2))
plot(m1)

m1.gls <- gls(prod~fert, data=duraznos)
m2.gls <- gls(prod~fert, weights=varIdent(form=~1|fert), data=duraznos)
AICtab(m1.gls,m2.gls)
anova(m2.gls)


###################################################
### code chunk number 3: Examen.Rnw:244-247
###################################################
# Validacion del modelo gls
residm <- resid(m2.gls, type="normalized")
qqnorm(residm)


###################################################
### code chunk number 4: Examen.Rnw:248-250
###################################################
plot(fitted(m2.gls), residm, xlab="Valores predichos", 
ylab= "Residuos estandarizados")


###################################################
### code chunk number 5: Examen.Rnw:252-256
###################################################
#cuales grupos difieren entre si
summary(m2.gls) 
m3.gls <- gls(prod~fert-1, weights=varIdent(form=~1|fert)) 
confint(m3.gls)


###################################################
### code chunk number 6: Examen.Rnw:288-315
###################################################
# Vainas con un numero de semillas, variables
set.seed(3333) 
lambd <- 5
# numero de vainas colectadas, una por árbol 
nv <- 130
# numero de semillas por vaina 
nsem <- rpois(nv,5) 
table(nsem)
# Grosor de la vaina 
gv <- runif(nv,0.1, 0.8)
# valores de los parametros
a <- 4
b <- -8

logistica <- function(x) { exp(x)/(1 + exp(x)) }
y_det <- logistica(a+b*gv) 
# plot(gv,y_det)
# Semillas depredadas
nsemdep <- rbinom(nv, size=nsem,prob=y_det)
# Base de datos
acacia <- data.frame(gv, nsem, nsemdep)
# Ajuste del modelo
m1 <- glm(cbind(nsemdep, nsem-nsemdep)~gv,binomial)
summary(m1)
# Validacion del modelo
par(mfrow=c(2,2)) 
plot(m1)


###################################################
### code chunk number 7: Examen.Rnw:317-330
###################################################
coefi <- coef(m1) 
a.est <- coefi[1] 
b.est <- coefi[2]
par(las=1)
plot(gv, nsemdep/nsem, cex = nsem/10, xlab= "Grosor vaina (mm)", 
ylab= "Probabilidad de depredación") 
curve(plogis(a.est + b.est*x), add=TRUE, col="red")
curve(plogis(a + b*x), add=TRUE, col="blue")
legend("topright", c("coef. sim.", "coef. est."), 
col=c("blue", "red"), lty=c(1,1))
plogis(a.est + b.est*0.7)
vp07 <- list(gv=c(0.3,0.5,0.7))
predict(m1, vp07, type="response", se=T)


###################################################
### code chunk number 8: Examen.Rnw:375-394
###################################################
## linear mixed effects, mucho suenho 8hrs, poco suenho 4hrs, vigilia 
set.seed(1111)
nind <- 100 
est <- gl(3,1,300, labels=c("ms","ns", "ps")) 
ind <- rep(1:nind, each=3)
mms <- 28 
dems <- 2
ind.means <- rnorm(nind,mms,dems) 
eps <- rnorm(300, 0,3) 
x <- rep(1:nind, rep(3,100)) 
X <- as.matrix(model.matrix(~as.factor(x)-1))
yal <- as.numeric( X %*% as.matrix(ind.means) + eps) 
efs <- rep(c(0,-5,-1),100) 
om <- yal + efs
om <- round(om, digits=1)
vigilia <- data.frame(ind,est, om)
m1 <- lme(om ~ est, random=~1|ind, data=vigilia) 
summary(m1)
anova(m1)


###################################################
### code chunk number 9: Examen.Rnw:396-398
###################################################
#Validación del modelo
plot(m1)


###################################################
### code chunk number 10: Examen.Rnw:400-401
###################################################
qqnorm(m1)


###################################################
### code chunk number 11: Examen.Rnw:404-413
###################################################
library(multcomp)
library(effects)
efec <- allEffects(m1)
summary(efec)
cm <-glht(m1, linfct = mcp(est= "Tukey")) 
iccm<- confint(cm)
par(mar=c(5,8,3,2), las=1) 
plot(iccm)
save(acacia,duraznos,vigilia, file="bdatos.RData")