---
title: "Crossed-effects model for orange tree circumference"
output: pdf_document
---
```{r,eval=F,include=FALSE}
#Doesn't work. Error appears to occur in MakeADFun
#gdbsource("gdbScript.R",TRUE)
```  
  
## Scaled up version of the Orange Tree example
```{r}
require(TMB)
compile("orangeCrossed.cpp")
dyn.load(dynlib("orangeCrossed"))
```

# Load data and create obj using muliple copies
```{r Input data}
# Read data
source("orange_data.R")
#pdf("OrangeData.pdf",width=5,height=3)
par(mar=c(4,4,1,4),las=1)
plot(y~t,data=Oranges,xlab="Days",ylab="Circumference (cm)")
Tree=rep(1:Oranges$M,Oranges$ngroup)
for(i in 1:Oranges$M) 
  points(Oranges$t[Tree==i],Oranges$y[Tree==i],type="l",lty=2)
#graphics.off()
```

```{r Create extra trees}
#Create jittered copies of data
nextra=0
if(nextra>0) {
ymatrix=matrix(Oranges$y,nrow=Oranges$M,byrow=T)
tmatrix=matrix(Oranges$t,nrow=Oranges$M,byrow=T)
yextra=matrix(NA,nrow=nextra,ncol=7,byrow=T)
textra=matrix(NA,nrow=nextra,ncol=7,byrow=T)
for(i in 1:nextra) {
  tree=sample(1:Oranges$M,1)
  yextra[i,]=round(ymatrix[tree,]*runif(1,0.9,1.1)+rnorm(7,0,7))
  textra[i,]=tmatrix[tree,]                 
}
N=Oranges$M+nextra
yN=rbind(ymatrix,yextra)
tN=rbind(tmatrix,textra)
Oranges<-list(n=N*7,y=as.vector(t(yN)),t=as.vector(t(tN)),M=N,ngroup=rep(7,N))
plot(y~t,data=Oranges)
Tree=rep(1:Oranges$M,Oranges$ngroup)
for(i in 1:Oranges$M) 
  points(Oranges$t[Tree==i],Oranges$y[Tree==i],type="l",lty=2)
}
```

```{r Model fitting}
obj <- MakeADFun(data=Oranges,DLL="orangeCrossed",
              parameters=list(beta=c(0,0,0),log_sigma=1,log_sigma_u=2,log_sigma_v=2,
                              u = rep(0,Oranges$M),v=rep(0,7)),random=c("u","v"),silent=T)
obj$control <- list(trace=0,parscale=c(1,1,1,1,1,1),reltol=1e-12,maxit=100)
obj$hessian <- T
newtonOption(obj,smartsearch=TRUE)
opt <- do.call("optim",obj)
opt


obj$report()
rep <- sdreport(obj)
summary(rep)

```

```{r,eval=F,include=F}
# Bayesian fit (assuming flat priors)
obj <- MakeADFun(data=Oranges,DLL="orangeCrossed",
              parameters=list(beta=c(0,0,0),log_sigma=1,log_sigma_u=2,log_sigma_v=2,
                              u = rep(0,Oranges$M),v=rep(0,7)),silent=T)
obj$control <- list(trace=0,reltol=1e-12,maxit=100)
obj$hessian <- T
newtonOption(obj,smartsearch=TRUE)
opt <- do.call("optim",obj)
opt

par(mfrow=c(2,3)) 
samps1 <- mcmc(obj,1000,"RWM",alpha=0.25)
samps1=samps2[seq(1,nrow(samps1),10),] #Thinning
for(i in 1:6) acf(samps1[,i])
for(i in 1:6) plot(density(samps1[,i]))

par(mfrow=c(2,3)) 
samps2 <- mcmc(obj,10000,"HMC",L=1)
samps2=samps2[seq(1,nrow(samps2),10),] #Thinning
for(i in 1:6) acf(samps2[,i])
for(i in 1:6) plot(density(samps2[,i]))

#samps3 <- mcmc(obj,1000,"NUTS") #Very slow (20 mins for 1000)
#for(i in 1:5) acf(samps3[,i])

```