#Exercise_1 Script 5 - Examples of Regression using Meta Data #Execute Load Observed Data.R before running this script require(lme4) require(lmerTest) #############Linear mixed effect model with PubID as Random ######################## #data weighted using the N of obs per treatment sqrt(o$N_study) #RUP flow predictions based on a single DC lmod <- lmer(Obs_RUPIn ~ Dt_CPIn + (1|PubID), data=o, weights = sqrt(N_study), REML=FALSE) summary(lmod) ############## Nonlinear LMER Model with PubID as Random ########################### parms1 <- c("Int","KpA","KpB","KpC") inputs1 <- c("Dt_CPAIn","Dt_CPBIn","Dt_CPCIn","KdRUP") args1 <- append(inputs1, parms1) init1 <- c(Int = 0.0, KpA = 0.2, KpB = 2.0, KpC=1.0) #lowerbound <- c(KpC=0, Int=-.2, KpFor=0.0, KpConc=0) #upperbound <- c(KpC=1.0, Int=0.7, KpFor=10, KpConc=12) #formnlmer is the model to predict the observed values. Documentation indicates the Int is required. formnlmer <- ~ Int + Dt_CPAIn * KpA + Dt_CPBIn * KpB/(KdRUP + KpB) + Dt_CPCIn * KpC form.d <- deriv(formnlmer, parms1, function.arg = args1) #the derivative of formnlmer nlmod1 <- nlmer(Obs_RUPIn ~ form.d(Dt_CPAIn, Dt_CPBIn, Dt_CPCIn, KdRUP,Int, KpA, KpB, KpC) ~ (Int|PubID), data = o, #weights = sqrt(N_study), start = init1, verbose=1, nlmerControl(optimizer = "Nelder_Mead")) #Allows data weighting as for lmer, but wont' converge with weights = sqrt(N_study) summary(nlmod1) #Calculate RMSE and CCC for each prediction to determing goodness of fit obs <- o$Obs_RUPIn predl <- predict(lmod, re.form=NA) #predict w/out study effects (re.form=NA) predlran <- predict(lmod) #predict w study effects prednl1 <- predict(nlmod1) #re.form option does not work for nlin models #load and execute the RMSE_CCC Funciton script to make the RMSE fuction available RMSE(obs,predl) RMSE(obs,predlran) RMSE(obs,prednl1) plot(obs,prednl1) #use the resprgrplt function from Plot Residuals Functions.R resprgrpplt(obs,predl,o$PubID,"RUP (linear DC)") resprgrpplt(obs,predlran,o$PubID,"RUP (linear DC)") resprgrpplt(obs,prednl1,o$PubID,"RUP (Kd/Kp)")