MJacobs1985

try<-nlsList(y~f1(x,U,D,Kd), 
              data=Grass_new,
              start=list(U=rep(19.7,2),D=rep(52,2),Kd=0.035)) # the (x,2) splits the variable in two
summary(try)
plot(try)
plot(intervals(try))
plot(try, Trial ~ resid(.), abline = 0 ) # Because a single concentration curve is used for all subjects, the individual differences noticed in Figure 6.3 are incorporated in the residuals, thus inflating the residual standard error
plot(try, as.factor(Cow) ~ resid(.), abline = 0 )
plot(try, as.factor(Replicate) ~ resid(.), abline = 0 )
plot(fixef(try))

MJacobs1985

#### NLIN - Use identification variable to already look into NLMixed
# Simple model, ignoring groups
try1<-nls(y~f1(x,U,D,Kd), data=(Grass[complete.cases(Grass$y),]),start=list(U=19.7,D=52,Kd=0.035))
summary(try1)

# Full model with different parameters for the two groups (which makes it the larger model, because of more parameters)
try2<-nls(y~f1(x,
               U[Standard.sample2],
               D[Standard.sample2],
               Kd[Standard.sample2]), 

MJacobs1985

Grass_new<-groupedData(y~x|Trial/Code/id,
                       data=Grass2,
                       outer=~Standard.sample2,
                       FUN=mean,
                       labels=list(x="Time", y="DMres"))
                       
fit<-nlme(y~U+D*exp(-(Kd*x)),
                data=Grass_new, 
                fixed=list(U~Standard.sample2+Trial,D+Kd~1), 
                random=U+D~1,

MJacobs1985

fit2<-nlme(y~U+D*exp(-(Kd*x)),
                data=Grass_new, 
                fixed=U+D+Kd~1, 
                random=U+D+Kd~1|Trial,
                start=c(U=20,D=50.88,Kd=0.02698),
                method="ML", 
                control=nlmeControl(maxIter=100000, msMaxIter = 100))
plot(augPred(fit2, level=0:1)) 

fit3<-nlme(y~U+D*exp(-(Kd*x)),

MJacobs1985

Grass_new<-groupedData(y~x|Trial,
                       data=Grass2,
                       FUN=mean,
                       labels=list(x="Time", y="DMres"))      
                       
nlin<-nlsList(y~f1(x,U,D,Kd), 
              data=Grass_new,
              start=list(U=19.7,D=52,Kd=0.035))
              
nlmixed<-nlme(y~U+D*exp(-(Kd*x)),

MJacobs1985

Grass_new<-groupedData(y~x|Trial/Code,
                       data=Grass2,
                       outer=~Standard.sample2,
                       FUN=mean,
                       labels=list(x="Time", y="DMres"))                                    
                       
fit6<-nlme(y~U+D*exp(-(Kd*x)),
                data=Grass_new, 
                fixed=list(U+D+Kd~1), 
                random=pdDiag(U+D~1), 

MJacobs1985

# Individual variancies for the standard samples
# weights = varIdent(form = ~ 1 | Standard.sample2)

fit8<-nlme(y~U+D*exp(-(Kd*x)),
                data=Grass_new, 
                fixed=list(U+D+Kd~Standard.sample2), 
                random=pdDiag(U+D~1),
                weights=varIdent(form=~1|Standard.sample2),
                method="REML",
                start=c(20,0,50.88,0,0.02698,0))

MJacobs1985

Grass_new<-groupedData(y~x|id,
                       data=Grass2,
                       FUN=mean,
                       labels=list(x="Time", y="DMres"))    
nlm_auto<-nlme(y~U+D*exp(-(Kd*x)),
          data=Grass_new,
          fixed=U+D+Kd~1, 
          random=U+D+Kd~1,
          start=c(U=20,D=50.88,Kd=0.02698),
          correlation = corARMA(c(0.5, -0.4,-0.2,0.3), p=1, q=3, form=~1|id),

MJacobs1985

#### PREDICTIONS ####
Grass_new<-groupedData(y~x|Trial/Cow,
                       data=Grass2,
                       outer=~Standard.sample2,
                       FUN=mean,
                       labels=list(x="Time", y="DMres")) 
                       
fm1<-nlme(y~U+D*exp(-(Kd*x)),
          data=Grass_new,
          fixed=U+D+Kd~1, 

MJacobs1985

#### Plotting NLMIXED Model #####
try.nlm2<-nlme(y~U+D*exp(-(Kd*x)),
               data=Grass_compl,
               fixed=U+D+Kd~1, 
               random=U+D+Kd~1|Trial/Cow,
               start=c(U=20,D=50.88,Kd=0.02698),
               method="ML", 
               control=nlmeControl(maxIter=1000, msMaxIter = 100))
               
pdata<-expand.grid(x=as.numeric(0:350),