MJacobs1985

## Plot Power curve
ggplot(powervalue, 
       aes(x=N,
           y=Power,
           col=ES,
           group=ES))+
  geom_line() + 
  theme_bw()+
  scale_x_discrete(limits=seq(1,35,1))+
  geom_hline(yintercept=80, linetype="dashed", color = "red", size=0.5)+

MJacobs1985

rm(list = ls())
#### Libraries ####
library(lme4)
library(ggplot2)
library(nlme)
library(MASS)
library(nls2)
library(nlstools)
library(car)
library(lattice)

MJacobs1985

#### Datafile ####
Grass <- read.csv("..../Grass.csv")
Grass$id<-paste(Grass$Code,Grass$Replicate)  # Create unique ID's for each curve
plot(Grass$Time, Grass$Dmres)

MJacobs1985

#### NL REGRESSION ####
## Formula 
formuDMres<-as.formula(y~U+D*exp(-(Kd*x)))
f1<-function(x,U,D,Kd){
  U+D*exp(-(Kd*x))
}

## Preview plots
plot(Grass_compl$x, Grass_compl$y)
curve(f1(x, U=19.7, D=52, Kd=0.035), add=T, col="red", lwd=3) # values provided by researcher

MJacobs1985

## Include identification variable - sample. This will influence if and how to split U, D or Kd

# Simplest way --> Split analysis --> SMALL model is the pooled data. LARGE model is the separate data combined
try<-nls(formula=formuDMres, data=(Grass[complete.cases(Grass$y),]),start=list(U=19.7,D=52,Kd=0.035))
try.SS1<-nls(formula=formuDMres, data=Grass[Grass$Standard.sample=="1",],start=list(U=19.7,D=52,Kd=0.035))
try.SS2<-nls(formula=formuDMres, data=Grass[Grass$Standard.sample=="2",],start=list(U=19.7,D=52,Kd=0.035))

summary(try.SS1)
summary(try.SS2)

MJacobs1985

Grass_new<-groupedData(y~x|id, 
                       data=Grass2,
                       FUN=mean,
                       labels=list(x="Time", y="DMres"))
summary(gsummary(Grass_new)) 
plot(Grass_new) 

MJacobs1985

## Tests for Nested Models - Largest model needs to go last
# ANOVA F-Tests 
anova(try1,try3,try2) 
anova(try3,try2)
anova(try1,try2)

# Manual F-Test
resid.small<-resid(try1)
resid.big<-resid(try2)
df.small<-summary(try1)$df

MJacobs1985

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))

try.nlm3<-nlme(y~U+D*exp(-(Kd*x)),
               data=Grass_compl,

MJacobs1985

Grass_new<-groupedData(y~x|Trial/Code,
                       data=Grass2,
                       FUN=mean,
                       labels=list(x="Time", y="DMres"))   
try.nlm1<-nlme(y~U+D*exp(-(Kd*x)),
               data=Grass_new,
               fixed=U+D+Kd~1, 
               random=U+D~1|Trial/Code,
               start=c(U=20,D=50.88,Kd=0.02698),
               method="REML", 

MJacobs1985

colnames(Grass)[4]<-"x"
colnames(Grass)[9]<-"y"
unique(Grass$x)
length(unique(Grass$x))
table(Grass$x)

Grass<-Grass[order(Grass$id),]
Grass$Standard.sample2<-factor(Grass$Standard.sample)
levels(Grass$Standard.sample2)