## Playing around with the Soybean data

library(nlme)
data(Soybean)

fit1 <- nlsList(weight ~ SSlogis(Time, Asym, xmid, scal) ,data = Soybean)

fit.nlme1 <- nlme(fit1)

## Adding the effect of Variety

fix <- fixef(fit.nlme1)
fit.nlme2 <- update(fit.nlme1, fixed = Asym + xmid + scal ~ Variety,
                    start = c(fix[1], 0, fix[2], 0, fix[3], 0), data = Soybean)

## Checking assumptions of mixed model
## I'll pick fit.nlme2

## Normality
pdf('../figs/normal-check.pdf')
qqnorm(fit.nlme2)
dev.off()

## Variances
pdf('../figs/variance-check.pdf')
plot(fit.nlme2)
dev.off()

## Independence of random effects
pairs(ranef(fit.nlme2))

intervals(fit.nlme2)

## Fixing the variance assumption

fit.nlme3 <- update(fit.nlme2, weights = varPower())

plot(ranef(fit.nlme3))

pairs(ranef(fit.nlme3))

## Simpler model
fit.nlme4 <- update(fit.nlme3 , random = pdDiag(Asym + xmid + scal ~ 1))

intervals(fit.nlme4 , which = "var-cov")

fit.nlme5 <- update(fit.nlme4 , random = pdDiag(Asym + xmid ~ 1))

## We also need to check for independence

pdf('../figs/auto-correlation.pdf')
plot(ACF(fit.nlme5), alpha = 0.05)
dev.off()

plot(augPred(fit.nlme5, level = 0:1))