Mixed models and lme4

Mixed Models and MCMC

Model selection codes
Model 1
Initial model specification
library(lme4)
Model1=lmer(yield.per.ha~n.per.ha*rain*soilclass+manure+seed+n.squared+(1|year)+(1|hhid)+(1|aezsmall)+(1|aez),data=isfm)
Model selection using fitLMER.fnc
library(LMERConvenienceFunctions)
optimum1=fitLMER.fnc(Model1,backfit.on="F",item=F,alpha=0.05,if.warn.not.add=TRUE,llrt=T,prune.ranefs=TRUE,p.value="upper",t.threshold=2,set.REML.FALSE=TRUE,reset.REML.TRUE=TRUE)
optimum1
Model criticism plots
mcp.fnc(optimum1)
library(languageR)
Computation of p-values
Using MCMC
pvals.fnc(optimum1)
Using log likelihood ratio test
Model 1
p.values.lmer<- function(x) {
summary.model <- summary(x)
data.lmer <- data.frame(model.matrix(x))
names(data.lmer) <- names(fixef(x))
names(data.lmer) <- gsub(pattern=":", x=names(data.lmer), replacement=".", fixed=T)
names(data.lmer) <- ifelse(names(data.lmer)=="(Intercept)", "Intercept", names(data.lmer))
string.call <- strsplit(x=as.character(x@call), split=" + (", fixed=T)
var.dep <- unlist(strsplit(x=unlist(string.call)[2], " ~ ", fixed=T))[1]
vars.fixef <- names(data.lmer)
formula.ranef <- paste("+ (", string.call[[2]][-1], sep="")
formula.ranef <- paste(formula.ranef, collapse=" ")
formula.full <- as.formula(paste(var.dep, "~ -1 +", paste(vars.fixef, collapse=" + "),
formula.ranef))
data.ranef <- data.frame(x@frame[,
which(names(x@frame) %in% names(ranef(x)))])
names(data.ranef) <- names(ranef(x))
data.lmer <- data.frame(x@frame[, 1], data.lmer, data.ranef)
names(data.lmer)[1] <- var.dep
out.full <- lmer(formula.full, data=data.lmer, REML=F)
p.value.LRT <- vector(length=length(vars.fixef))
for(i in 1:length(vars.fixef)) {
formula.reduced <- as.formula(paste(var.dep, "~ -1 +", paste(vars.fixef[-i],
collapse=" + "), formula.ranef))
out.reduced <- lmer(formula.reduced, data=data.lmer, REML=F)
print(paste("Reduced by:", vars.fixef[i]))
print(out.LRT <- data.frame(anova(out.full, out.reduced)))
p.value.LRT[i] <- round(out.LRT[2, 7], 3)
}
summary.model@coefs <- cbind(summary.model@coefs, p.value.LRT)
summary.model@methTitle <- c("\n", summary.model@methTitle,
"\n(p-values from comparing nested models fit by maximum likelihood)")
print(summary.model)
}
library(lme4)
lmer.final.p<-lmer(yield.per.ha~n.per.ha+rain+soilclass+manure+seed+n.squared+n.per.ha:rain+rain:soilclass+(1|year)+(1|hhid)+(1|aezsmall)+(1|aez),data=isfm)
p.values.lmer(lmer.final.p)
Model 2
Initial model specification
library(lme4)
Model2=lmer(yield.per.ha~p.per.ha*rain*soilclass+manure+seed+p.squared+(1|year)+(1|hhid)+(1|aezsmall)+(1|aez),data=isfm)
Model2
library(LMERConvenienceFunctions)
Model selection using fitLMER.fnc
optimum2=fitLMER.fnc(Model2,backfit.on="F",item=F,alpha=0.05,if.warn.not.add=TRUE,llrt=T,prune.ranefs=TRUE,p.value="upper",t.threshold=2,set.REML.FALSE=TRUE,reset.REML.TRUE=TRUE)
optimum2
mcp.fnc(optimum2)
library(languageR)
Computation of p-values
Using MCMC
pvals.fnc(optimum2)
Using log likelihood ratios
p.values.lmer<- function(x) {
summary.model <- summary(x)
data.lmer <- data.frame(model.matrix(x))
names(data.lmer) <- names(fixef(x))
names(data.lmer) <- gsub(pattern=":", x=names(data.lmer), replacement=".", fixed=T)
names(data.lmer) <- ifelse(names(data.lmer)=="(Intercept)", "Intercept", names(data.lmer))
string.call <- strsplit(x=as.character(x@call), split=" + (", fixed=T)
var.dep <- unlist(strsplit(x=unlist(string.call)[2], " ~ ", fixed=T))[1]
vars.fixef <- names(data.lmer)
formula.ranef <- paste("+ (", string.call[[2]][-1], sep="")
formula.ranef <- paste(formula.ranef, collapse=" ")
formula.full <- as.formula(paste(var.dep, "~ -1 +", paste(vars.fixef, collapse=" + "),
formula.ranef))
data.ranef <- data.frame(x@frame[,
which(names(x@frame) %in% names(ranef(x)))])
names(data.ranef) <- names(ranef(x))
data.lmer <- data.frame(x@frame[, 1], data.lmer, data.ranef)
names(data.lmer)[1] <- var.dep
out.full <- lmer(formula.full, data=data.lmer, REML=F)
p.value.LRT <- vector(length=length(vars.fixef))
for(i in 1:length(vars.fixef)) {
formula.reduced <- as.formula(paste(var.dep, "~ -1 +", paste(vars.fixef[-i],
collapse=" + "), formula.ranef))
out.reduced <- lmer(formula.reduced, data=data.lmer, REML=F)
print(paste("Reduced by:", vars.fixef[i]))
print(out.LRT <- data.frame(anova(out.full, out.reduced)))
p.value.LRT[i] <- round(out.LRT[2, 7], 3)
}
summary.model@coefs <- cbind(summary.model@coefs, p.value.LRT)
summary.model@methTitle <- c("\n", summary.model@methTitle,
"\n(p-values from comparing nested models fit by maximum likelihood)")
print(summary.model)
}

library(lme4)
lmer.final.p<-lmer(yield.per.ha~p.per.ha+rain+soilclass+manure+seed+p.per.ha:rain+p.per.ha:soilclass+rain:soilclass+p.per.ha:rain:soilclass+(1|year)+(1|hhid)+(1|aezsmall)+(1|aez),data=isfm)
p.values.lmer(lmer.final.p)
Model 3
Initial model specification
library(lme4)
Model3=lmer(yield.per.ha~k.per.ha*rain*soilclass+manure+seed+k.squared+(1|year)+(1|hhid)+(1|aezsmall)+(1|aez),data=isfm)
Model3
Model selection fitLMER.fnc
library(LMERConvenienceFunctions)
optimum3=fitLMER.fnc(Model3,backfit.on="F",item=F,alpha=0.05,if.warn.not.add=TRUE,llrt=FALSE,prune.ranefs=TRUE,p.value="upper",t.threshold=2,set.REML.FALSE=TRUE,reset.REML.TRUE=TRUE)
optimum3
library(languageR)
Computation of p-values
Using MCMC
pvals.fnc(optimum3)
Using log likelihood ratios
p.values.lmer<- function(x) {
summary.model <- summary(x)
data.lmer <- data.frame(model.matrix(x))
names(data.lmer) <- names(fixef(x))
names(data.lmer) <- gsub(pattern=":", x=names(data.lmer), replacement=".", fixed=T)
names(data.lmer) <- ifelse(names(data.lmer)=="(Intercept)", "Intercept", names(data.lmer))
string.call <- strsplit(x=as.character(x@call), split=" + (", fixed=T)
var.dep <- unlist(strsplit(x=unlist(string.call)[2], " ~ ", fixed=T))[1]
vars.fixef <- names(data.lmer)
formula.ranef <- paste("+ (", string.call[[2]][-1], sep="")
formula.ranef <- paste(formula.ranef, collapse=" ")
formula.full <- as.formula(paste(var.dep, "~ -1 +", paste(vars.fixef, collapse=" + "),
formula.ranef))
data.ranef <- data.frame(x@frame[,
which(names(x@frame) %in% names(ranef(x)))])
names(data.ranef) <- names(ranef(x))
data.lmer <- data.frame(x@frame[, 1], data.lmer, data.ranef)
names(data.lmer)[1] <- var.dep
out.full <- lmer(formula.full, data=data.lmer, REML=F)
p.value.LRT <- vector(length=length(vars.fixef))
for(i in 1:length(vars.fixef)) {
formula.reduced <- as.formula(paste(var.dep, "~ -1 +", paste(vars.fixef[-i],
collapse=" + "), formula.ranef))
out.reduced <- lmer(formula.reduced, data=data.lmer, REML=F)
print(paste("Reduced by:", vars.fixef[i]))
print(out.LRT <- data.frame(anova(out.full, out.reduced)))
p.value.LRT[i] <- round(out.LRT[2, 7], 3)
}
summary.model@coefs <- cbind(summary.model@coefs, p.value.LRT)
summary.model@methTitle <- c("\n", summary.model@methTitle,
"\n(p-values from comparing nested models fit by maximum likelihood)")
print(summary.model)
}

library(lme4)
lmer.final.p<-lmer(yield.per.ha~k.per.ha+rain+manure+seed+k.squared+k.per.ha:rain+(1|year)+(1|hhid)+(1|aezsmall)+(1|aez),data=isfm)
p.values.lmer(lmer.final.p)