#R Code for UWisc Example 1
#Updated 6.16.17
#With Questions Contact:
#	Robin White
#	rrwhite@vt.edu
#	509-701-9290

#___________________________________________Open Packages_______________________________________________
library(ggplot2)
library(reshape2)
library(plyr)
library(lme4)
library(lmerTest)

#______________________________________Open Google Drive____________________________________________
d <- read.csv("WorkshopData.csv")

#_______________________________________Visualize The Data______________________________________________
d1 <- d[c("TrtID","BW", "Milk", "DMI", "OMI", "Nkg", "NDFkg", "fMicNdu", "fMicNduSE", "rpH", "rNH3", "rVFA", "rAcet", "rProp", "rButr")]

mlt <- melt(d1, id="TrtID")
ggplot(mlt, aes(x=value))+geom_density()+facet_wrap(~variable, ncol=8, scales="free")

#_______________________________________Correct Outliers_________________________________________________
summary(d$rVFA)
d[d$rVFA>200,"PubID"]
d[d$PubID==42,]

#Assume this is a slipped decimal
d[d$TrtID==37,"rVFA"] <- 74.64
d[d$PubID==42,]

#_______________________________________Correct SEM Errors_______________________________________
SEM <- tapply(d$fMicNduSE, d$StatMethod, mean, na.rm=TRUE)[2]
SEF <- tapply(d$fMicNduSE, d$StatMethod, mean, na.rm=TRUE)[1]
d$StatSEM <- ifelse(d$StatMethod=="Mixed",  d$fMicNduSE/SEM, d$fMicNduSE/SEF)
d$Trunc <- ifelse(d$StatSEM<mean(d$StatSEM, na.rm=TRUE)*0.5,mean(d$StatSEM, na.rm=TRUE)*0.5,d$StatSEM)
d$MicN_Wt <- 1/d$Trunc
pct <- ifelse(d$StatSEM<mean(d$StatSEM, na.rm=TRUE)*0.5,1,0)
nonna <- ifelse(is.na(d$StatSEM)==TRUE, 0, 1)
sum(pct, na.rm=TRUE)/sum(nonna, na.rm=TRUE)

#Reduce cutoff to prevent excessive curtailing
d$Trunc <- ifelse(d$StatSEM<mean(d$StatSEM, na.rm=TRUE)*0.25,mean(d$StatSEM, na.rm=TRUE)*0.25,d$StatSEM)
d$MicN_Wt <- 1/d$Trunc
pct <- ifelse(d$StatSEM<mean(d$StatSEM, na.rm=TRUE)*0.25,1,0)
nonna <- ifelse(is.na(d$StatSEM)==TRUE, 0, 1)
sum(pct, na.rm=TRUE)/sum(nonna, na.rm=TRUE)


#__________________________________________Derive Model, Phase 1_________________________________________________
#Include all variables
m1 <- lmer(fMicNdu~DMI+BW+Nkg+NDFkg+LiqMarker+SolMarker+MicrMarker+SMPLoc+(1|PubID), data=d, weight=MicN_Wt)
#Remove DMI with p-value 0.955
m1 <- lmer(fMicNdu~BW+Nkg+NDFkg+LiqMarker+SolMarker+MicrMarker+SMPLoc+(1|PubID), data=d, weight=MicN_Wt)
#Remove NDF intake with p-value 0.641
m1 <- lmer(fMicNdu~BW+Nkg+LiqMarker+SolMarker+MicrMarker+SMPLoc+(1|PubID), data=d, weight=MicN_Wt)
#Remove SolMarker with p-value 0.815
m1 <- lmer(fMicNdu~BW+Nkg+LiqMarker+MicrMarker+SMPLoc+(1|PubID), data=d, weight=MicN_Wt)
#Remove BW with p-value 0.383
m1 <- lmer(fMicNdu~Nkg+LiqMarker+MicrMarker+SMPLoc+(1|PubID), data=d, weight=MicN_Wt)
#Remove Liq Marker with p-value 0.36 and 0.20
m1 <- lmer(fMicNdu~Nkg+MicrMarker+SMPLoc+(1|PubID), data=d, weight=MicN_Wt)

#__________________________________________Re-test Dropped Parameters_________________________________________________
#Test DMI again
m1 <- lmer(fMicNdu~Nkg+MicrMarker+SMPLoc+DMI+(1|PubID), data=d, weight=MicN_Wt)
#remove Nkg with p-value 0.392
m2 <- lmer(fMicNdu~MicrMarker+SMPLoc+DMI+(1|PubID), data=d, weight=MicN_Wt)
#remove MicrMarker with p-value 0.1915
m2 <- lmer(fMicNdu~SMPLoc+DMI+(1|PubID), data=d, weight=MicN_Wt)

#Test NDF again
m1 <- lmer(fMicNdu~Nkg+MicrMarker+SMPLoc+NDFkg+(1|PubID), data=d, weight=MicN_Wt)
#Test SolMarker again
m1 <- lmer(fMicNdu~Nkg+MicrMarker+SMPLoc+SolMarker+(1|PubID), data=d, weight=MicN_Wt)
#Test BW again
m1 <- lmer(fMicNdu~Nkg+MicrMarker+SMPLoc+BW+(1|PubID), data=d, weight=MicN_Wt)

#__________________________________________Check VIF_________________________________________________

#Make a function to return variance inflation factors
VIF <- function (fit) {
v <- vcov(fit)
nam <- names(fixef(fit))
ns <- sum(1 * (nam == "Intercept" | nam == "(Intercept)"))
if (ns > 0) {
v <- v[-(1:ns), -(1:ns), drop = FALSE]
nam <- nam[-(1:ns)]
}
d <- diag(v)^0.5
v <- diag(solve(v/(d %o% d)))
names(v) <- nam
return(v)
}

#Check VIF of model 1
VIF(m1)

#Check VIF of model 2
VIF(m2)

#__________________________________________Check Model Fit_________________________________________________

#Generate an equation to estimate model fit statistics
RMSE <- function(m) {
  m <- m
  o <- attributes(m)$frame[,1]
  p<-fitted(m)
  rand <- ranef(m)$PubID
  PubID <- attributes(ranef(m)$PubID)$row.names
  temp <- data.frame(PubID, rand)
  new <- merge(attributes(m)$frame, temp, by="PubID")
  p2 <- fitted(m)-new$X.Intercept.
  meano <- mean(o, na.rm=TRUE)
  meanp <- mean(p, na.rm=TRUE)
  res=o-p
  resad = o-p2
  res2=res^2
  resad2 = resad^2
  rm=sqrt(mean(res2, na.rm=TRUE));
  rmad=sqrt(mean(resad2, na.rm=TRUE))
  uss=sum(res2, na.rm=TRUE);
  lo <- ifelse(is.na(o)==FALSE, 1, 0)
  n=sum(lo);
  meanO=mean(o, na.rm=TRUE);
  mb=sum(res, na.rm=TRUE)/n;
  sse <- anova(lm(res~p))[2,2];
  msb <- mb^2;
  mspe <- rm^2;
  msre <- sse/n;
  msslope <- mspe-msre-msb;
  mean <- msb/mspe*100;
  slope <- msslope/mspe*100;
  residual <- msre/mspe*100;
  check <- mean+slope+residual
  rsr <- rmad/sd(o, na.rm=TRUE)
  ccc <- epi.ccc(o,p)$rho.c[1]
  cccad <- epi.ccc(o,p2)$rho.c[1]
  rmp = rm/meano*100
  rmpad = rmad/mean(p2, na.rm=TRUE)*100
  mb <- mean(res, na.rm=TRUE)
  sb <- coef(lm(res~p))[2]
  aicc <- AICc(m)
  Shats <- attributes(VarCorr(m)$PubID)$stddev
  Shate <- attributes(VarCorr(m))$sc
  Shats_e <- Shats/Shate
  output <- format(c(n,meano, meanp, rm, rmad, rmp, rmpad, mean, slope, residual, check, mb, sb, rsr,ccc[,1],cccad[,1], aicc, Shats, Shate, Shats_e ), scientific=FALSE)
  labels <- c("N", "Observed Mean", "Predicted Mean", "RMSE, units", "Unadjusted RMSE, units", "RMSE, % mean", "Unadjusted RMSE, % mean", "Mean Bias, % MSE", "Slope Bias, % MSE", "Residual Error", "Error Check", "Mean Bias", "Slope Bias", "RSR", "CCC", "Unadjusted CCC", "AICc", "Sigma Hat Study", "Sigma Hat Error", "Ratio of Sigma Hat Study/Error")
  out <- data.frame(labels,output)
  temp <- NULL
  new <- NULL
  return(out)
}

#Open additional required packages
library(epiR)
library(MuMIn)


#Output fit statistics on m1
RMSE(m1)

#Output fit statistics on m2
RMSE(m2)

#~~~~~~~~~~~~~~~~Cross Validate Models~~~~~~~~~~~~~~~~~~~~~
RMSE2 <- function(o,p) {
o <- o
p <- p
meano <- mean(o, na.rm=TRUE)
meanp <- mean(p, na.rm=TRUE)
res=o-p;
res2=res^2; 
rm=sqrt(mean(res2, na.rm=TRUE));
uss=sum(res2, na.rm=TRUE);
lo <- ifelse(is.na(o)==FALSE, 1, 0)
n=sum(lo);
meanO=mean(o, na.rm=TRUE);
mb=sum(res, na.rm=TRUE)/n;
sse <- anova(lm(res~p))[2,2];
msb <- mb^2;
mspe <- rm^2;
msre <- sse/n;
msslope <- mspe-msre-msb;
mean <- msb/mspe*100;
slope <- msslope/mspe*100;
residual <- msre/mspe*100;
check <- mean+slope+residual
rsr <- rm/sd(o, na.rm=TRUE)
ccc <- epi.ccc(o,p)$rho.c[1]
rmp = rm/meanp*100
mb <- mean(res, na.rm=TRUE)
sb <- coef(lm(res~p))[2]
output <- format(c(meano, meanp, rmp, mean, slope, residual, check, mb, sb, rsr,ccc[,1]), scientific=FALSE)
labels <- c("Observed Mean", "Predicted Mean", "RMSE, % mean", "Mean Bias, % MSE", "Slope Bias, % MSE", "Residual Error", "Error Check", "Mean Bias", "Slope Bias", "RSR", "CCC")
out <- data.frame(labels,output)
return(out)
}


CV_CP <- function(end,dpct) {
i <- 0
while(i<end) {

#Create a sequential ID vector
d$nID <- seq(1:length(d$TrtID))

#Count how many treatments to include in the derivaiton dataset
DerN <- round((dpct/100)*length(d$nID),0)

#Select data from the derivation dataset
DerID <- sample(1:length(d$TrtID), DerN)
d$dCheck <- ifelse(d$nID %in% DerID, 1,0)

#Subset data into derivation and evaluation data
dev <- subset(d, d$nID %in% DerID)
eval <- subset(d, d$dCheck==0)
m1 <- lmer(fMicNdu~Nkg+MicrMarker+SMPLoc+(1|PubID), data=d, weight=MicN_Wt)
parms <- fixef(m1)
eval$MicMarker1 <- ifelse(eval$MicrMarker== "NucAcids", 1, 0)
eval$SMPLoc1 <- ifelse(eval$SMPLoc=="Omasal", 1, 0)
pred <- fixef(m1)[1]+fixef(m1)[2]*eval$Nkg+fixef(m1)[3]*eval$MicMarker1+fixef(m1)[4]*eval$SMPLoc1

rmspe <- RMSE2(eval$fMicNdu, pred)[3,2]
ccc <- RMSE2(eval$fMicNdu, pred)[11,2]


#Save the parameter estimates and significance values
row <- c(fixef(m1), as.numeric(as.character(rmspe)), as.numeric(as.character(ccc)))
if(i==0) data=row else data=rbind(data, row)
i <- i+1
}

return(data)
}
o <- CV_CP(10,70)
o <- data.frame(o)


bettersum <- function(d) {
grid <- data.frame(id=NA, n=NA, mean=NA, median=NA, sd=NA, min=NA, Q1=NA, Q3=NA, max=NA)
for(i in 1:length(names(d))) {
id <- names(d)[i]
n <- ifelse(is.na(d[,i])==FALSE, 1,0)
n <- sum(n)
Mean <- mean(d[,i], na.rm=TRUE)
Median <- median(d[,i], na.rm=TRUE)
Q1 <- quantile(d[,i], na.rm=TRUE)[2]
Q3 <- quantile(d[,i], na.rm=TRUE)[4]
SD <- sd(d[,i], na.rm=TRUE)
min <- min(d[,i], na.rm=TRUE)
max <- max(d[,i], na.rm=TRUE)
row <- data.frame(id,n, Mean, Median, SD, min, Q1, Q3, max)
grid[i,] <- row
}
grid$id <- names(d)
return(grid) 
}


bettersum(o)