monogenea

lmm1 <- lme(total.fruits ~ rack + nutrient + amd + status,
            random = ~1|reg, method = "ML")
lmm2 <- lme(total.fruits ~ rack + nutrient + amd + status,
            random = ~1|popu, method = "ML")
lmm3 <- lme(total.fruits ~ rack + nutrient + amd + status,
            random = ~1|gen, method = "ML")
lmm4 <- lme(total.fruits ~ rack + nutrient + amd + status,
            random = ~1|reg/popu, method = "ML")
lmm5 <- lme(total.fruits ~ rack + nutrient + amd + status,
            random = ~1|reg/gen, method = "ML")

monogenea

# Set optimization pars
ctrl <- lmeControl(opt="optim")
lmm6.2 <- update(lmm6, .~., random = ~nutrient|popu/gen, control = ctrl)
lmm7.2 <- update(lmm7, .~., random = ~nutrient|reg/popu/gen, control = ctrl)
anova(lmm6, lmm6.2, lmm7, lmm7.2) # both models improved
anova(lmm6.2, lmm7.2) # similar fit; lmm6.2 more parsimonious
summary(lmm6.2)

monogenea

plot(ranef(lmm6.2, level = 2))

monogenea

# QQ plots (drawn to the same scale!)
par(mfrow = c(1,2))
lims <- c(-3.5,3.5)
qqnorm(resid(GLM, type = "normalized"),
       xlim = lims, ylim = lims,main = "GLM")
abline(0,1, col = "red", lty = 2)
qqnorm(resid(lmm6.2, type = "normalized"),
       xlim = lims, ylim = lims, main = "lmm6.2")
abline(0,1, col = "red", lty = 2)

monogenea

summary(lmm6.2)

monogenea

lmm8 <- update(lmm6.2, .~. + nutrient:amd)
summary(lmm8)
anova(lmm8, lmm6.2)

monogenea

finalModel <- update(lmm6.2, .~., method = "REML")
summary(finalModel)

monogenea

dev.off() # Reset previous graphical pars
# New GLM, updated from the first by estimating with REML
GLM2 <- update(GLM, .~., method = "REML")
# Plot side by side, beta with respective SEs
plot(coef(GLM2), xlab = "Fixed Effects", ylab = expression(beta), axes = F,
     pch = 16, col = "black", ylim = c(-.9,2.2))
stdErrors <- coef(summary(GLM2))[,2]
segments(x0 = 1:6, x1 = 1:6, y0 = coef(GLM2) - stdErrors, y1 = coef(GLM2) + stdErrors,
         col = "black")
axis(2)

monogenea

plot(finalModel)

monogenea

library(snpStats)
load("conversionTable.RData")

pathM <- paste("public/Genomics/108Malay_2527458snps", c(".bed", ".bim", ".fam"), sep = "")
SNP_M <- read.plink(pathM[1], pathM[2], pathM[3])

pathI <- paste("public/Genomics/105Indian_2527458snps", c(".bed", ".bim", ".fam"), sep = "")
SNP_I <- read.plink(pathI[1], pathI[2], pathI[3])

pathC <- paste("public/Genomics/110Chinese_2527458snps", c(".bed", ".bim", ".fam"), sep = "")