library(MASS)
library(glmmTMB)
library(ggplot2)
library(readr)

timber2 <- read_delim("https://sciences.ucf.edu/biology/d4lab/wp-content/uploads/sites/23/2017/01/timber2.txt", 
                delim = "\t", escape_double = FALSE, col_types = cols(volume = col_number(), 
                girth = col_number(), height = col_number()), trim_ws = TRUE)
attach(timber2)

model<-glmmTMB(volume~species+height*girth,data=timber2,family=gaussian(link="log"))
summary(model)
    

##### SIMPLE PREDICTION FOR ONE SPECIES (MEAN HEIGHT)
nd<-data.frame(species=rep("A",100),
               height=rep(mean(height[species=="A"]),100),
               girth=seq(0,max(girth),length.out=100))

nd$vol<-predict(model,newdata=nd)

plot(girth[species=="A"],volume[species=="A"],xlab="girth",ylab="volume")
lines(nd$girth,exp(nd$vol),lwd=2)

#### PREDICT BOTH SPECIES (MEAN HEIGHT FOR EACH SPECIES)

# Only means
nd<-data.frame(species=c(rep("A",100),rep("B",100)),
               height=c(rep(mean(height[species=="A"]),100),rep(mean(height[species=="B"]),100)),
               girth=rep(seq(0,max(girth),length.out=100),2))

nd$vol<-predict(model,newdata=nd)

par(mfrow=c(1,2))

plot(girth[species=="A"],volume[species=="A"],xlab="girth",ylab="volume",main="Species A")
lines(nd$girth[nd$species=="A"],exp(nd$vol[nd$species=="A"]),lwd=2)

plot(girth[species=="B"],volume[species=="B"],xlab="girth",ylab="volume",main="Species B")
lines(nd$girth[nd$species=="B"],exp(nd$vol[nd$species=="B"]),lwd=2)

# Means and CI

preds<-predict(model,newdata=nd,se.fit=T)
nd$vol<-exp(preds$fit)
nd$upCI<-exp(preds$fit+qnorm(0.975)*preds$se.fit)
nd$loCI<-exp(preds$fit-qnorm(0.975)*preds$se.fit)

par(mfrow=c(1,2))

plot(girth[species=="A"],volume[species=="A"],xlab="girth",ylab="volume",main="Species A")
lines(nd$girth[nd$species=="A"],nd$vol[nd$species=="A"],lwd=2)
lines(nd$girth[nd$species=="A"],nd$upCI[nd$species=="A"],lty=2)
lines(nd$girth[nd$species=="A"],nd$loCI[nd$species=="A"],lty=2)

plot(girth[species=="B"],volume[species=="B"],xlab="girth",ylab="volume",main="Species B")
lines(nd$girth[nd$species=="B"],nd$vol[nd$species=="B"],lwd=2)
lines(nd$girth[nd$species=="B"],nd$upCI[nd$species=="B"],lty=2)
lines(nd$girth[nd$species=="B"],nd$loCI[nd$species=="B"],lty=2)



#### PREDICT BOTH SPECIES WITH CI FOR DIFFERENT HEIGHTS (GGPLOT example)

hght<-rep(as.numeric(rep(quantile(timber2$height),each=50)),2)
grth<-rep(seq(0,max(girth),length.out=50),10)
spc<-rep(c("A","B"),each=250)


nd<-data.frame(species=spc,height=hght,girth=grth)

preds<-predict(model,newdata=nd,se.fit = TRUE)

nd$volume<-exp(preds$fit)
nd$lowCI<-exp(preds$fit-qnorm(0.975)*preds$se.fit)
nd$upCI<-exp(preds$fit+qnorm(0.975)*preds$se.fit)
nd$Height<-rep(rep(as.character(quantile(timber2$height)),each=50),2) # This is to let ggplot use height as factors


p<-ggplot(timber2, aes(girth,volume))+
  geom_smooth(data=nd,mapping=aes(girth,volume,colour=Height))+
  geom_point()
p+facet_grid(cols=vars(species))


p<-ggplot(timber2, aes(girth,volume))+
  geom_smooth(data=nd,mapping=aes(girth,volume,colour=species))+
  geom_point()+
  geom_ribbon(data=nd,aes(ymin = lowCI, ymax = upCI, fill = species), 
              alpha=0.4)
p+facet_grid(cols=vars(Height))


##### WITH sjPlot

library(sjPlot)
library(sjmisc)

plot_model(model,type="pred",terms = "girth")
plot_model(model,type="pred",terms = c("girth","species"))
plot_model(model,type="pred",terms = c("girth","height"))
plot_model(model,type="pred",terms = c("girth","height","species"))