#Derive Milk Protien Prediction
require(lme4)
#Data centering function; use, c.(Year) to center Year
c. <- function (x) scale(x, scale = FALSE)			#Center the x var yielding c.x

o <- subset(o, !is.na(o$Obs_MilkTP_g))  #reduce the data set to only those with milk protein
ls(o)

#Create a sum of squared individual AA terms
o$Abs_EAA2_g <- o$Abs_Arg_g^2 + o$Abs_His_g^2 + o$Abs_Ile_g^2 + o$Abs_Leu_g^2 + 
      o$Abs_Lys_g^2 + o$Abs_Met_g^2 + o$Abs_Phe_g^2 + o$Abs_Thr_g^2 + o$Abs_Trp_g^2 + 
      o$Abs_Val_g^2
 
#An_DEInp is the DEI from all nutrients except protein
mTPmod <- lmer(Obs_MilkTP_g ~ An_DEInp + Abs_Arg_g + Abs_His_g + Abs_Ile_g + Abs_Leu_g + 
               Abs_Lys_g + Abs_Met_g + Abs_Phe_g + Abs_Thr_g + Abs_Trp_g + Abs_Val_g + 
               Abs_EAA2_g + Parity_rl + c.(Year) + (1|PubID), 
               data=o, weights = sqrt(N_study), REML=FALSE)
summary(mTPmod)

#Reduce the model by backward elimination
mTPmod <- lmer(Obs_MilkTP_g ~ An_DEInp + Abs_Arg_g + Abs_His_g + Abs_Ile_g + Abs_Leu_g + 
                    Abs_Lys_g + Abs_Met_g + Abs_Trp_g + Abs_Val_g + 
                    Abs_EAA2_g + (1|PubID), 
               data=o, weights = sqrt(N_study), REML=FALSE)
summary(mTPmod)

pr_mlk_TP_g <- predict(mTPmod, re.form=NA)     #predict w/out study effects (re.form=NA)
length(pr_mlk_TP_g)
new <- cbind(o, pr_mlk_TP_g)  #this is a dangerous merge as there is no check on data correspondence

ls(new)

#load and execute the RMSE_CCC Funciton script to make the RMSE fuction available
RMSE(new$Obs_MilkTP_g, new$pr_mlk_TP_g)

plot(new$pr_mlk_TP_g,new$Obs_MilkTP_g)
resprgrpplt(new$Obs_MilkTP_g, new$pr_mlk_TP_g, new$PubID,"Milk Protein, g/d")