.Rapp.history

1-pexp(a, rate=1/E.x)
plot(x, 1/10*exp(-(1/10)*x)); abline(v=1-pexp(a, rate=1/E.x))
plot(x, 1/10*exp(-(1/10)*x)); abline(v=pexp(a, rate=1/E.x))
runif(1, 1, 6)
p.die
sum(p.die)
x.die = 1:6
p.die = rep(1/6, 6)
x.die* p.die
mean(x.die* p.die)
sum(x.die* p.die)
mean.die <- sum(x.die* p.die)
mean.die/6
mean.die/10
out.die <- 1:6
f.x(out.die)
f.x <- function(x){#
	return(x*0.5)#
}#
out.die <- 1:6#
f.x(out.die)
mean(f.x(out.die))
f.x(mean(out.die))
fnonlin.x <- function(x){#
	return(x*2)#
}
plot(out.die, fnonlin.x(out.die))
## Linear function#
flin.x <- function(x){#
	return(x*0.5)#
}#
out.die <- 1:6#
mean(flin.x(out.die))#
flin.x(mean(out.die))#
#
## Nonlinear function#
fnonlin.x <- function(x){#
	return(x^2)#
}#
plot(out.die, fnonlin.x(out.die))
plot(out.die, fnonlin.x(out.die), type="o")
mean(fnonlin.x(out.die))
fnonlin.x(mean(out.die))
fnonlin.x(mean(out.die))==fnonlin.x(mean(out.die))
fnonlin.x(mean(out.die))==mean(fnonlin.x(out.die))
flin.x(mean(out.die))==mean(flin.x(out.die))
out.urv <- runif(50,1,6)
out.urv <- round(runif(50,1 , 6))
out.urv <- sample(1:6, 1, 100))
out.urv <- sample(1:6, 1, 100)
out.urv <- sample(1:6, 50, 100)
out.urv <- sample(1:6, 50, replace=T)
out.urv
## generate 50 unif rv#
out.urv <- sample(1:6, 50, replace=T)#
n_trials <- 100#
n_samples <- 100#
for(i in 1:n_samples){#
	outcomes = sample(1:6, n_samples, replace=T)#
	res      = fnonlin.x(outcomes)#
	v1       = mean(res)#
	v2       = fnonlin.x(mean(outcomes))#
	paste0(i, " = ", v1, " > ", v2)#
}
n_trials <- 100#
n_samples <- 100#
for(i in 1:n_samples){#
	outcomes = sample(1:6, n_samples, replace=T)#
	res      = fnonlin.x(outcomes)#
	v1       = mean(res)#
	v2       = fnonlin.x(mean(outcomes))#
	print(paste0(i, " = ", v1, " > ", v2))#
}
## generate 50 unif rv#
out.urv <- sample(1:6, 50, replace=T)#
n_trials <- 100#
n_samples <- 100#
for(i in 1:n_samples){#
	outcomes = sample(1:6, n_samples, replace=T)#
	res      = fnonlin.x(outcomes)#
	v1       = mean(res)#
	v2       = fnonlin.x(mean(outcomes))#
	print(paste0(i, " = ", v1, " > ", v2, " is ", v1> v2))#
}
## generate 50 unif rv#
out.urv <- sample(1:6, 50, replace=T)#
n_trials <- 10#
n_samples <- 100#
for(i in 1: n_trials){#
	outcomes = sample(1:6, n_samples, replace=T)#
	res      = fnonlin.x(outcomes)#
	v1       = mean(res)#
	v2       = fnonlin.x(mean(outcomes))#
	print(paste0(i, " = ", v1, " > ", v2, " is ", v1> v2))#
}
### arithmetic mean >= geometric mean#
flog.x <- function(x){#
	return(log(x))#
}
### arithmetic mean >= geometric mean#
flog.x <- function(x){#
	return(log(x))#
}#
out.die <- 1:6#
mean(flog.x(out.die))#
flog.x(mean(out.die))
flog.x(mean(out.die))==mean(flog.x(out.die))
plot(out.die, flog.x(out.die), type="o")
library(extRemes)
install.packages("extRemes")
library(fitdistrplus)#
library(moments)
? gofstat
# how many seeds each bird gets on a simulation run, which is equivalent to how many fruits in a foraging session in the same tree#
#
# nseeds <- round(runif(1, 3, 6))#
nseeds <- 5#
#
# Sample a random Gut Retention Time for each one of the seeds#
grt_tt <- round(rgamma(nseeds, shape = 4, scale = 5) + 8)#
#
# Maximum gut retention time determines the simulation time#
tt <- max(grt_tt)#
#
movrate <- 28#
#
# Sample movement distances from an exponential distribution. These are the lengths that the animal moves at each time step in minutes.#
movedist <- rexp(tt, 1/movrate)#
#
# Movement algorithm for this random walk. Sample a random angle of movement from a uniform distribution. Determine the movement in x and y, the animal's trajectory is the cumulative sum of those. We are not imposing any limits to the area of movement.#
angle <- runif(tt, min = 0, max = 360)#
distx <- movedist*cos(angle)#
xloc <- c(0, cumsum(distx))#
disty <- movedist*sin(angle)#
yloc <- c(0, cumsum(disty))#
animalTraj <- data.frame(time = 0:tt, xloc = xloc, yloc = yloc)#
#
# Determine where the seeds are dropped based on the animal's location.#
seed_loc <- data.frame(seed_id = 1:nseeds, time = grt_tt)#
seed_loc <- merge(animalTraj, seed_loc, by = "time")#
seed_loc$dispersal <- sqrt(seed_loc$xloc^2 + seed_loc$yloc^2)
# Determine where the seeds are dropped based on the animal's location.#
seed_loc <- data.frame(seed_id = 1:nseeds, time = grt_tt)#
seed_loc <- merge(animalTraj, seed_loc, by = "time")#
seed_loc$dispersal <- sqrt(seed_loc$xloc^2 + seed_loc$yloc^2)#
#
outsim <- list(lengths = movedist,#
               movement = animalTraj,#
              seed = seed_loc)#
#
# Calculate the mean location of seeds#
x_m <- mean(outsim$seed$xloc)#
y_m <- mean(outsim$seed$yloc)#
mean_dispersal <- mean(outsim$seed$dispersal)#
se_dispersal <- sd(outsim$seed$dispersal)/sqrt(length(outsim$seed$dispersal))#
#
# Calculate seed dispersion#
xi <- (x_m - outsim$seed$xloc)^2#
yi <- (y_m - outsim$seed$yloc)^2#
seed_dispersion <- sum(sqrt(xi + yi))/length(outsim$seed$time)#
par(mfrow = c(1,3))#
hist(movedist, main = "MD")#
plot(x = xloc, y = yloc, type = "l", main = "Movement")#
points(x = seed_loc$xloc, y = seed_loc$yloc, col = "blue", pch = 16)#
hist(seed_loc$dispersal, main = "Dispersal")#
abline(v = mean_dispersal, col = "red")
sim_run <- function(movrate){#
#
  nseeds <- 5#
  grt_tt <- round(rgamma(nseeds, shape = 4, scale = 5) + 8)#
  tt <- max(grt_tt)#
#
  movedist <- rexp(tt, 1/movrate)#
#
  angle <- runif(tt, min = 0, max = 360)#
  distx <- movedist*cos(angle)#
  xloc <- c(0, cumsum(distx))#
  disty <- movedist*sin(angle)#
  yloc <- c(0, cumsum(disty))#
  animalTraj <- data.frame(time = 0:tt, xloc = xloc, yloc = yloc)#
#
  seed_loc <- data.frame(seed_id = 1:nseeds, time = grt_tt)#
  seed_loc <- merge(animalTraj, seed_loc, by = "time")#
  seed_loc$dispersal <- sqrt(seed_loc$xloc^2 + seed_loc$yloc^2)#
#
  x_m <- mean(seed_loc$xloc)#
  y_m <- mean(seed_loc$yloc)#
  mean_xy_dispersal <- sqrt((x_m^2) + (y_m^2))#
  mean_dispersal <- mean(seed_loc$dispersal)#
  se_dispersal <- sd(seed_loc$dispersal)/sqrt(length(seed_loc$dispersal))#
#
  # Calculate seed dispersion#
  xi <- (x_m - seed_loc$xloc)^2#
  yi <- (y_m - seed_loc$yloc)^2#
  seed_dispersion <- sum(sqrt(xi + yi))/length(seed_loc$time)#
#
  outsim <- list(lengths = movedist,#
                 movement = animalTraj,#
                 seed = seed_loc,#
                 mean_xy_dispersal = mean_xy_dispersal,#
                 mean_dispersal = mean_dispersal,#
                 se_dispersal = se_dispersal,#
                 seed_dispersion = seed_dispersion)#
#
  return(outsim)#
}#
onerun <- sim_run(28)#
#
############# Simulations for the population#
#
library(dplyr)#
library(purrr)#
library(ggplot2)#
library(purrr)#
library(cowplot)#
#
set.seed(227)#
#
data(ptpl)#
#
null_moverate <- data.frame(Bird_ID = unique(ptpl$Bird_ID), movrate = mean(ptpl$mpm))#
#
indiv_moverate <- ptpl %>%#
  group_by(Bird_ID) %>%#
  summarise(movrate = mean(mpm))#
#
fam_moverate <- ptpl %>%#
  group_by(fam_g) %>%#
  summarise(movrate = mean(mpm))#
#
ids <- ptpl %>% distinct(., Bird_ID, fam_g)#
#
#### Null model#
# Assume all individuals share the same movement rate, which is the average movement rate for all data points#
#
# Run 10,000 simulation runs per individual or family group#
#
nruns <- 10000#
#
null_dispersal <- NULL#
null_dispersion <- NULL#
for(j in 1:length(null_moverate$Bird_ID)){#
  moverate <- null_moverate$movrate[j]#
#
  manysims <- vector("list", nruns)#
  for(i in 1:length(manysims)){#
    manysims[[i]] <- sim_run(moverate)#
    manysims[[i]]$sim_run <- paste("sim_", i)#
  }#
#
  dispersal <- map_df(manysims, "seed")$dispersal#
  id <- null_moverate$Bird_ID[j]#
  mean_xy_dispersal <- map_dbl(manysims, "mean_xy_dispersal")#
  mean_dispersal <- map_dbl(manysims, "mean_dispersal")#
  se_dispersal <- map_dbl(manysims, "se_dispersal")#
  seed_dispersion <- map_dbl(manysims, "seed_dispersion")#
#
  out <- data.frame(dispersal = dispersal, id = id)#
  null_dispersal <- rbind.data.frame(null_dispersal, out)#
#
  out2 <- data.frame(id = id,#
                     mean_xy_dispersal = mean_xy_dispersal,#
                     mean_dispersal = mean_dispersal,#
                     se_dispersal = se_dispersal,#
                     seed_dispersion = seed_dispersion)#
  null_dispersion <- rbind.data.frame(null_dispersion, out2)#
#
}#
#
# Plots#
plot_grid(#
  null_dispersal %>%#
    ggplot(., aes(x = dispersal)) +#
    geom_histogram(),#
  plot_grid(null_dispersion %>%#
            ggplot(., aes(y = mean_dispersal)) +#
            geom_boxplot(),#
          null_dispersion %>%#
            ggplot(., aes(y = seed_dispersion)) +#
            geom_boxplot()), nrow = 2)#
### Individual simulation#
#
indiv_dispersal <- NULL#
indiv_dispersion <- NULL#
for(j in 1:length(indiv_moverate$Bird_ID)){#
  moverate <- indiv_moverate$movrate[j]#
#
  manysims <- vector("list", nruns)#
  for(i in 1:length(manysims)){#
    manysims[[i]] <- sim_run(moverate)#
    manysims[[i]]$sim_run <- paste("sim_", i)#
  }#
#
  dispersal <- map_df(manysims, "seed")$dispersal#
  id <- indiv_moverate$Bird_ID[j]#
  mean_xy_dispersal <- map_dbl(manysims, "mean_xy_dispersal")#
  mean_dispersal <- map_dbl(manysims, "mean_dispersal")#
  se_dispersal <- map_dbl(manysims, "se_dispersal")#
  seed_dispersion <- map_dbl(manysims, "seed_dispersion")#
#
  out <- data.frame(dispersal = dispersal, id = id)#
  indiv_dispersal <- rbind.data.frame(indiv_dispersal, out)#
#
  out2 <- data.frame(id = id,#
                     mean_xy_dispersal = mean_xy_dispersal,#
                     mean_dispersal = mean_dispersal,#
                     se_dispersal = se_dispersal,#
                     seed_dispersion = seed_dispersion)#
  indiv_dispersion <- rbind.data.frame(indiv_dispersion, out2)#
#
}#
#
# plots#
plot_grid(indiv_dispersal %>%#
            ggplot(., aes(x = dispersal)) +#
            geom_histogram(),#
          plot_grid(indiv_dispersion %>%#
                      ggplot(., aes(y = mean_dispersal)) +#
                      geom_boxplot(),#
                    indiv_dispersion %>%#
                      ggplot(., aes(y = seed_dispersion)) +#
                      geom_boxplot()), nrow = 2)#
### Social or Family group simulation#
#
fam_dispersal <- NULL#
fam_dispersion <- NULL#
for(j in 1:length(fam_moverate$fam_g)){#
  moverate <- fam_moverate$movrate[j]#
#
  manysims <- vector("list", nruns)#
  for(i in 1:length(manysims)){#
    manysims[[i]] <- sim_run(moverate)#
    manysims[[i]]$sim_run <- paste("sim_", i)#
  }#
#
  dispersal <- map_df(manysims, "seed")$dispersal#
  id <- fam_moverate$fam_g[j]#
  mean_xy_dispersal <- map_dbl(manysims, "mean_xy_dispersal")#
  mean_dispersal <- map_dbl(manysims, "mean_dispersal")#
  se_dispersal <- map_dbl(manysims, "se_dispersal")#
  seed_dispersion <- map_dbl(manysims, "seed_dispersion")#
#
  out <- data.frame(dispersal = dispersal, id = id)#
  fam_dispersal <- rbind.data.frame(fam_dispersal, out)#
#
  out2 <- data.frame(id = id,#
                     mean_xy_dispersal = mean_xy_dispersal,#
                     mean_dispersal = mean_dispersal,#
                     se_dispersal = se_dispersal,#
                     seed_dispersion = seed_dispersion)#
  fam_dispersion <- rbind.data.frame(fam_dispersion, out2)#
#
}#
#
# Plots#
plot_grid(#
  fam_dispersal %>%#
    ggplot(., aes(x = dispersal)) +#
    geom_histogram(),#
  plot_grid(fam_dispersion %>%#
            ggplot(., aes(y = mean_dispersal)) +#
            geom_boxplot(),#
          fam_dispersion %>%#
            ggplot(., aes(y = seed_dispersion)) +#
            geom_boxplot()), nrow = 2)
fm1 <- glmmTMB(count~mined+(1|spp),#
                  ziformula=~mined,#
                  data=Salamanders,#
                  family=nbinom1)
fm1R
library(boot)
summary(boot.ci(b1,type="perc"))
boot.ci(b1,type="perc")
fm1R <- refit(fm1, simulate(fm1)[[1]])
?refit
refit(fm1, simulate(fm1)[[1]])
refit(fm1, simulate(fm1)[[1]], 100)
refit(fm1, simulate(fm1)[[1]], 1000)
summary(refit(fm1, simulate(fm1)[[1]], 1000))
b1 <- lme4::bootMer(fm1, FUN=function(x) fixef(x)$zi, nsim=20, .progress="txt")#
   if (requireNamespace("boot")) {#
      boot.ci(b1,type="perc")#
    }
x2 <- sort(runif(100))#
yb2 <- rbinom(100, prob = plogis(2*(x2-1)), size = 1)#
yb2 <- factor(1 + yb2, labels = c("failure", "success"))#
modb2 <- glm(yb2 ~ x2, family = binomial)#
S4 <- simulate(modb2, nsim = 4)
S4
x <- 1:5#
mod1 <- lm(c(1:3, 7, 6) ~ x)#
S1 <- simulate(mod1, nsim = 4)#
## repeat the simulation:#
.Random.seed <- attr(S1, "seed")#
identical(S1, simulate(mod1, nsim = 4))#
#
S2 <- simulate(mod1, nsim = 200, seed = 101)#
rowMeans(S2) # should be about the same as#
fitted(mod1)
library(glmmTMB)
## CREATE A MODEL#
fm1 <- glmmTMB(count~mined+(1|spp),#
                  ziformula=~mined,#
                  data=Salamanders,#
                  family=nbinom1)
nrow  <- 100
nrow  <- dim(Salamanders)[1]
ncols <- 10#
N  <- dim(Salamanders)[1]#
boot.parm <- matrix(data=NA, ncol=ncols, nrow=N)
ncols <- 10
N     <- dim(Salamanders)[1]
head(Salamanders)
with(Salamanders, plot(count, Wtemp))
with(Salamanders, plot(DOY, count))
fm1 <- glmmTMB(count ~ DOY + (1|spp), data=Salamanders, family=nbinom1(link="log"))
exp(-0.05749)
with(Salamanders, plot(Wtemp, count))
with(Salamanders, plot(cover, count))
with(Salamanders, plot(cover, count))\
fm1 <- glmmTMB(count ~ cover + (1|spp), data=Salamanders, family=nbinom1(link="log"))
ncols     <- 10
N         <- dim(Salamanders)[1]
predict(fm1)
predict(fm1, data=nice.xs)
simulate(fm1)
library(glmmTMB)#
#
## CREATE A MODEL#
fm1 <- glmmTMB(count ~ cover + (1|spp), data=Salamanders, family=nbinom1(link="log"))#
#
## PARAMETRIC BOOTSTRAP#
ncols     <- 10#
N         <- dim(Salamanders)[1]#
boot.parm <- matrix(data=NA, ncol=ncols, nrow=N)#
nice.xs   <- sort(Salamanders$cover) #
for(i in 1:ncol){#
	refitted      <- refit(fm1, simulate(fm1)[[1]])#
	boot.parm[,i] <- predict(refitted, data=nice.xs)#
}
predict(refitted, data=nice.xs)
library(glmmTMB)#
#
## CREATE A MODEL#
fm1 <- glmmTMB(count ~ cover + (1|spp), data=Salamanders, family=nbinom1(link="log"))#
#
## PARAMETRIC BOOTSTRAP#
ncols     <- 10#
N         <- dim(Salamanders)[1]#
boot.parm <- matrix(data=NA, ncol=ncols, nrow=N)#
nice.xs   <- sort(Salamanders$cover)
refitted      <- refit(fm1, simulate(fm1)[[1]])
boot.parm[,i] <- predict(refitted, data=nice.xs)
1:ncol
ncols     <- 10#
N         <- dim(Salamanders)[1]#
boot.parm <- matrix(data=NA, ncol=ncols, nrow=N)#
nice.xs   <- sort(Salamanders$cover) #
for(i in 1:ncols){#
	refitted      <- refit(fm1, simulate(fm1)[[1]])#
	boot.parm[,i] <- predict(refitted, data=nice.xs)#
}
MLE.boot <- apply(boot.parm, 2, FUN=function(x){ exp(x[1]+x[2]*nice.xs) })
refitted
length(fixef(fm1))
fixef(fm1)
length(fixef(fm1)$cond)
## PARAMETRIC BOOTSTRAP#
ncols     <- length(fixef(fm1)$cond)#
N         <- dim(Salamanders)[1]#
boot.parm <- matrix(data=NA, ncol=ncols, nrow=N)#
nice.xs   <- sort(Salamanders$cover) #
for(i in 1:ncols){#
	refitted      <- refit(fm1, simulate(fm1)[[1]])#
	boot.parm[,i] <- fixef(refitted)$cond#
}
fm1
n.reps    <- 10000
## PARAMETRIC BOOTSTRAP#
ncols     <- length(fixef(fm1)$cond)#
N         <- dim(Salamanders)[1]#
n.reps    <- 10000 #
boot.parm <- matrix(data=NA, ncol=ncols, nrow=N)#
nice.xs   <- sort(Salamanders$cover) #
for(i in 1:ncols){#
	refitted      <- refit(fm1, simulate(fm1)[[1]], n.reps)#
	boot.parm[,i] <- fixef(refitted)$cond#
}
library(glmmTMB)#
#
## CREATE A MODEL#
fm1 <- glmmTMB(count ~ cover + (1|spp), data=Salamanders, family=nbinom1(link="log"))#
#
## PARAMETRIC BOOTSTRAP#
ncols     <- length(fixef(fm1)$cond)#
N         <- dim(Salamanders)[1]#
n.reps    <- 100000 #
boot.parm <- matrix(data=NA, ncol=ncols, nrow=N)#
nice.xs   <- sort(Salamanders$cover) #
for(i in 1:ncols){#
	refitted      <- refit(fm1, simulate(fm1)[[1]], n.reps)#
	boot.parm[,i] <- fixef(refitted)$cond#
}
? refit
MLE.boot <- apply(boot.parm, 1, FUN=function(x){ exp(x[1]+x[2]*nice.xs) })
MLE.boot <- apply(boot.parm, 1, FUN=function(x){ exp(x[1] + x[2] * nice.xs ) })
MLE.boot
Boot.cis <- apply(MLE.boot, 1,FUN=function(x){quantile(x,probs=c(0.025, 0.5, 0.975))})
nice.xs   <- sort(Salamanders$cover)
nice.xs   <- seq(min(Salamanders$cover), max(Salamanders$cover), length.out=length(Salamanders$cover))
nice.xs
Boot.cis
t(Boot.cis)
Boot.cis <- t(Boot.cis) # To put in the same format that Wald CI's matrix
library(glmmTMB)#
#
## CREATE A MODEL#
fm1 <- glmmTMB(count ~ cover + (1|spp), data=Salamanders, family=nbinom1(link="log"))#
#
## PARAMETRIC BOOTSTRAP#
ncols     <- length(fixef(fm1)$cond)#
N         <- dim(Salamanders)[1]#
n.reps    <- 100000 #
boot.parm <- matrix(data=NA, ncol=ncols, nrow=N)#
nice.xs   <- seq(min(Salamanders$cover), max(Salamanders$cover), length.out=length(Salamanders$cover)) #
for(i in 1:ncols){#
	refitted      <- refit(fm1, simulate(fm1)[[1]], n.reps)#
	boot.parm[,i] <- fixef(refitted)$cond#
}#
#
MLE.boot <- apply(boot.parm, 1, FUN=function(x){ exp(x[1] + x[2] * nice.xs ) })#
Boot.cis <- apply(MLE.boot, 1, FUN=function(x){quantile(x,probs=c(0.025, 0.5, 0.975))})#
Boot.cis <- t(Boot.cis) # To put in the same format that Wald CI's matrix
MLE
exp(MLE[1]+ MLE[2]*nice.xs)
beta.hat <- exp(MLE[1]+ MLE[2]*nice.xs)
beta.hat <- exp(MLE[1]+ MLE[2]*nice.xs)#
plot(Salamanders$cover, Salamanders$count, #
     type="p", lwd=2, bty="l", xlab="", ylab="", pch=16, cex.lab=1.75, xaxt = "n", cex.axis= 1.5)
plot(Salamanders$cover, Salamanders$count, type="p", lwd=2, bty="l", xlab="", ylab="", pch=16, cex.lab=1.75, cex.axis= 1.5)
lines(nice.xs, p.hat, col="grey10")
lines(nice.xs, beta.hat, col="grey10")
points(nice.xs, Boot.cis[,1], type="l", col="grey40", lwd=1.5, lty="dashed")
points(nice.xs, Boot.cis[,3], type="l", col="grey40", lwd=1.5, lty="dashed")
Salamanders
beta.hat <- exp(MLE[1]+ MLE[2]*nice.xs)#
plot(Salamanders$cover, Salamanders$count, type="p", lwd=2, bty="l", xlab="", ylab="", pch=16, cex.lab=1.75, cex.axis= 1.5)#
lines(nice.xs, beta.hat, col="grey10")#
points(nice.xs, Boot.cis[,1], type="l", col="grey40", lwd=1.5, lty="dashed")#
points(nice.xs, Boot.cis[,3], type="l", col="grey40", lwd=1.5, lty="dashed")
Salamanders$counts
Salamanders$count <-
Salamanders$count
Salamanders[which(Salamanders$count>20),]
Salamanders <- Salamanders[-which(Salamanders$count>20),]
## CREATE A MODEL#
fm1 <- glmmTMB(count ~ cover + (1|spp), data=Salamanders, family=nbinom1(link="log"))#
#
## PARAMETRIC BOOTSTRAP#
ncols     <- length(fixef(fm1)$cond)#
N         <- dim(Salamanders)[1]#
n.reps    <- 100000 #
boot.parm <- matrix(data=NA, ncol=ncols, nrow=N)#
nice.xs   <- seq(min(Salamanders$cover), max(Salamanders$cover), length.out=length(Salamanders$cover)) #
for(i in 1:ncols){#
	refitted      <- refit(fm1, simulate(fm1)[[1]], n.reps)#
	boot.parm[,i] <- fixef(refitted)$cond#
}#
#
MLE      <- fixef(fm1)$cond#
MLE.boot <- apply(boot.parm, 1, FUN=function(x){ exp(x[1] + x[2] * nice.xs ) })#
Boot.cis <- apply(MLE.boot, 1, FUN=function(x){quantile(x,probs=c(0.025, 0.5, 0.975))})#
Boot.cis <- t(Boot.cis) # To put in the same format that Wald CI's matrix#
beta.hat <- exp(MLE[1]+ MLE[2]*nice.xs)#
plot(Salamanders$cover, Salamanders$count, type="p", lwd=2, bty="l", xlab="", ylab="", pch=16, cex.lab=1.75, cex.axis= 1.5)#
lines(nice.xs, beta.hat, col="grey10")#
points(nice.xs, Boot.cis[,1], type="l", col="grey40", lwd=1.5, lty="dashed")#
points(nice.xs, Boot.cis[,3], type="l", col="grey40", lwd=1.5, lty="dashed")
## CREATE A MODEL#
fm1 <- glmmTMB(count ~ cover + (1|spp), data=Salamanders, family=nbinom1(link="log"))#
#
## PARAMETRIC BOOTSTRAP#
ncols     <- length(fixef(fm1)$cond)#
N         <- dim(Salamanders)[1]#
n.reps    <- 100000 #
boot.parm <- matrix(data=NA, ncol=ncols, nrow=N)#
nice.xs   <- seq(min(Salamanders$cover), max(Salamanders$cover), length.out=length(Salamanders$cover))
for(i in 1:ncols){#
	refitted      <- refit(fm1, simulate(fm1)[[1]], n.reps)#
	boot.parm[,i] <- fixef(refitted)$cond#
}
MLE      <- fixef(fm1)$cond
dim(Salamanders)[1]
dim(Salamanders)
library(glmmTMB)#
#
## remove the outlier counts#
Salamanders <- Salamanders[-which(Salamanders$count>20),]#
## CREATE A MODEL#
fm1 <- glmmTMB(count ~ cover + (1|spp), data=Salamanders, family=nbinom1(link="log"))#
#
## PARAMETRIC BOOTSTRAP#
ncols     <- length(fixef(fm1)$cond)#
N         <- dim(Salamanders)[1]#
n.reps    <- 100000 #
boot.parm <- matrix(data=NA, ncol=ncols, nrow=N)#
nice.xs   <- seq(min(Salamanders$cover), max(Salamanders$cover), length.out=length(Salamanders$cover)) #
for(i in 1:ncols){#
	refitted      <- refit(fm1, simulate(fm1)[[1]], n.reps)#
	boot.parm[,i] <- fixef(refitted)$cond#
}#
#
MLE      <- fixef(fm1)$cond#
MLE.boot <- apply(boot.parm, 1, FUN=function(x){ exp(x[1] + x[2] * nice.xs ) })#
Boot.cis <- apply(MLE.boot, 1, FUN=function(x){quantile(x,probs=c(0.025, 0.5, 0.975))})#
Boot.cis <- t(Boot.cis) # To put in the same format that Wald CI's matrix#
beta.hat <- exp(MLE[1]+ MLE[2]*nice.xs)#
plot(Salamanders$cover, Salamanders$count, type="p", lwd=2, bty="l", xlab="", ylab="", pch=16, cex.lab=1.75, cex.axis= 1.5)#
lines(nice.xs, beta.hat, col="grey10")#
points(nice.xs, Boot.cis[,1], type="l", col="grey40", lwd=1.5, lty="dashed")#
points(nice.xs, Boot.cis[,3], type="l", col="grey40", lwd=1.5, lty="dashed")
library(glmmTMB)#
#
## remove the outlier counts#
dat <- Salamanders[-which(Salamanders$count>20),]#
## CREATE A MODEL#
fm1 <- glmmTMB(count ~ cover + (1|spp), data=dat, family=nbinom1(link="log"))
library(glmmTMB)#
#
## remove the outlier counts#
dat <- Salamanders[-which(Salamanders$count>20),]#
## CREATE A MODEL#
fm1 <- glmmTMB(count ~ cover + (1|spp), data=dat, family=nbinom1(link="log"))#
#
## PARAMETRIC BOOTSTRAP#
ncols     <- length(fixef(fm1)$cond)#
N         <- dim(dat)[1]#
n.reps    <- 100000 #
boot.parm <- matrix(data=NA, ncol=ncols, nrow=N)#
nice.xs   <- seq(min(dat$cover), max(dat$cover), length.out=length(dat$cover)) #
for(i in 1:ncols){#
	refitted      <- refit(fm1, simulate(fm1)[[1]], n.reps)#
	boot.parm[,i] <- fixef(refitted)$cond#
}#
#
MLE      <- fixef(fm1)$cond#
MLE.boot <- apply(boot.parm, 1, FUN=function(x){ exp(x[1] + x[2] * nice.xs ) })#
Boot.cis <- apply(MLE.boot, 1, FUN=function(x){quantile(x,probs=c(0.025, 0.5, 0.975))})#
Boot.cis <- t(Boot.cis) # To put in the same format that Wald CI's matrix
n.reps    <- 100000
boot.parm <- matrix(data=NA, ncol=ncols, nrow=N)
nice.xs   <- seq(min(dat$cover), max(dat$cover), length.out=length(dat$cover))
dat
Salamanders[-which(Salamanders$count>20),]
Salamanders[!which(Salamanders$count>20),]
Salamanders[which(Salamanders$count<20),]
dat <- Salamanders[which(Salamanders$count<20),]
fm1 <- glmmTMB(count ~ cover + (1|spp), data=dat, family=nbinom1(link="log"))
ncols     <- length(fixef(fm1)$cond)
N         <- dim(dat)[1]
N
## remove the outlier counts#
dat <- Salamanders[which(Salamanders$count<20),]#
## CREATE A MODEL#
fm1 <- glmmTMB(count ~ cover + (1|spp), data=dat, family=nbinom1(link="log"))#
#
## PARAMETRIC BOOTSTRAP#
ncols     <- length(fixef(fm1)$cond)#
N         <- dim(dat)[1]#
n.reps    <- 100000 #
boot.parm <- matrix(data=NA, ncol=ncols, nrow=N)#
nice.xs   <- seq(min(dat$cover), max(dat$cover), length.out=length(dat$cover))
for(i in 1:ncols){#
	refitted      <- refit(fm1, simulate(fm1)[[1]], n.reps)#
	boot.parm[,i] <- fixef(refitted)$cond#
}#
#
MLE      <- fixef(fm1)$cond#
MLE.boot <- apply(boot.parm, 1, FUN=function(x){ exp(x[1] + x[2] * nice.xs ) })#
Boot.cis <- apply(MLE.boot, 1, FUN=function(x){quantile(x,probs=c(0.025, 0.5, 0.975))})#
Boot.cis <- t(Boot.cis) # To put in the same format that Wald CI's matrix
refit(fm1, simulate(fm1)[[1]], n.reps)
fixef(refitted)$cond
boot.parm
## PARAMETRIC BOOTSTRAP#
ncols     <- length(fixef(fm1)$cond)#
N         <- dim(dat)[1]#
n.reps    <- 100000 #
boot.parm <- matrix(data=NA, ncol=ncols, nrow=N)#
nice.xs   <- seq(min(dat$cover), max(dat$cover), length.out=length(dat$cover)) #
for(i in 1:N){#
	refitted      <- refit(fm1, simulate(fm1)[[1]], n.reps)#
	boot.parm[i,] <- fixef(refitted)$cond#
}
library(glmmTMB)#
#
## remove the outlier counts#
dat <- Salamanders[which(Salamanders$count<20),]#
## CREATE A MODEL#
fm1 <- glmmTMB(count ~ cover + (1|spp), data=dat, family=nbinom1(link="log"))#
lme4::bootMer(fm1,nsim=100,FUN=function(x) unlist(fixef(x)))
merBoot <- lme4::bootMer(fm1, nsim=100, FUN=function(x) unlist(fixef(x)))
library(arm)
display(fm1)
merBoot
I.lower = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.025, na.rm=TRUE)))
CI.upper = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=.975, na.rm=TRUE)))
I.lower
CI.upper
cbind(lower  = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=0.025, na.rm=TRUE))),#
	  med    = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=0.5, na.rm=TRUE)))#
	  upper = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=0.975, na.rm=TRUE))))
cbind(lower  = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=0.025, na.rm=TRUE))),#
	  med    = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=0.5, na.rm=TRUE))),#
	  upper = apply(merBoot$t, 2, function(x) as.numeric(quantile(x, probs=0.975, na.rm=TRUE))))
print(fm1)
merBoot$t
## PARAMETRIC BOOTSTRAP#
ncols     <- length(fixef(fm1)$cond)#
N         <- dim(dat)[1]#
n.reps    <- 100#
boot.parm <- matrix(data=NA, ncol=ncols, nrow=N)#
nice.xs   <- seq(min(dat$cover), max(dat$cover), length.out=length(dat$cover)) #
#
for(i in 1:N){#
	refitted      <- refit(fm1, simulate(fm1)[[1]], n.reps)#
	boot.parm[i,] <- fixef(refitted)$cond#
}#
#
MLE      <- fixef(fm1)$cond#
MLE.boot <- apply(boot.parm, 1, FUN=function(x){ exp(x[1] + x[2] * nice.xs ) })#
Boot.cis <- apply(MLE.boot, 1, FUN=function(x){quantile(x,probs=c(0.025, 0.5, 0.975))})#
Boot.cis <- t(Boot.cis)
library(glmmTMB)#
library(merTools)#
#
## remove the outlier count#
dat <- Salamanders[which(Salamanders$count<20),]#
#
fm1 <- glmmTMB(count ~ cover + (1|spp), #
               data=dat, family=nbinom2(link="log"))#
fm2 <- glmer.nb(count ~ cover + (1|spp), data=dat)#
summary(fm1); summary(fm2)
source("/Users/georgeglen/Documents/GitHub/Pfaller_et_al2021/functions.r")#
#
## Load each dataset in individually#
setwd("~/Documents/GitHub/Pfaller_et_al2021")#
FULL <- fun$importExcel(#
  fileName = "nesting_dat.xlsx",#
  sheetNames = c("Recruitment","OCF", "ECF","RMI","BF","Production"),#
  nskip = 0#
);
## Seperate the datasets#
## 1: Clutch Frequency#
OCF <- FULL[[2]] %>% gather("Beach", "OCF")#
OCF$Beach <- str_replace_all(OCF$Beach, "OCF.", "")#
OCF <- OCF %>% transform(Beach = as.factor(Beach),#
					   OCF 	 = as.numeric(OCF)) %>%#
			 filter(!is.na(OCF))#
#
ECF <- FULL[[3]] %>% gather("Beach", "ECF")#
ECF$Beach <- str_replace_all(ECF$Beach, "ECF.", "")#
ECF <- ECF %>% transform(Beach = as.factor(Beach),#
					   ECF 	 = as.numeric(ECF)) %>%#
			 filter(!is.na(ECF))#
#
## 2: RMI#
RMI <- FULL[[4]] %>% gather("Beach", "RMI")#
RMI$Beach <- str_replace_all(RMI$Beach, "RI.", "")#
RMI <- RMI %>% transform(Beach = as.factor(Beach),#
						 RMI 	 = as.numeric(RMI)) %>%#
			   filter(!is.na(RMI))#
## 3: BP#
BF <- FULL[[5]] %>% gather("Beach", "BF")#
BF$Beach <- str_replace_all(BF$Beach, "BF.", "")#
BF <- BF %>% transform(Beach = as.factor(Beach),#
					   BF 	 = as.numeric(BF)) %>%#
			   filter(!is.na(BF))
#### ----------------------------------------- #####
#### Models#
#### ----------------------------------------- #####
#
## Write out the model formulas#
## 1: OCF#
null.OCF.form <- formula(OCF ~ 1)#
beach.OCF.form <- formula(OCF ~ Beach - 1)#
#
null.ECF.form <- formula(ECF ~ 1)#
beach.ECF.form <- formula(ECF ~ Beach - 1)#
#
## 2: RMI #
null.RMI.form <- formula(RMI ~ 1)#
beach.RMI.form <- formula(RMI ~ Beach - 1)#
#
## 3: BP#
null.BF.form <- formula(BF ~ 1)#
beach.BF.form <- formula(BF ~ Beach - 1)#
#
#### ----------------------------------------- #####
#### Model for OCF
## Create initial models#
null.OCF.Model <- glm(null.OCF.form, data = OCF, family=poisson(link="log"))#
OCF.Model      <- glm(beach.OCF.form, data = OCF, family=poisson(link="log"))#
#
## check dispersion#
testDispersion(null.OCF.Model); testDispersion(OCF.Model)
testDispersion(null.OCF.Model); testDispersion(OCF.Model)
## model to account for overdispersion but underdispersion is present so lets #
## use a gen.pois model in the glmmTMB package and do distribution selection#
#
ICtab(null.OCF.Model,#
	  glmmTMB(formula = null.OCF.form, data= OCF, family=genpois(link="log"), na.action = na.exclude), ## UPDATE THE MODEL DISTRIBUTIONS#
	  type = "AICc", delta = TRUE, base = TRUE, logLik= TRUE,#
	  weights=T, #
	  mnames=c("Poisson", "Generalized Poisson"))#
ICtab(OCF.Model,#
	  glmmTMB(formula = beach.OCF.form, data= OCF, family=genpois(link="log"),na.action = na.exclude,#
	          control=glmmTMBControl(optimizer=optim, optArgs=list(method="BFGS"))), ## UPDATE THE MODEL DISTRIBUTIONS#
	  type = "AICc", delta = TRUE, base = TRUE, logLik= TRUE,#
	  weights=T, #
	  mnames=c("Poisson", "Generalized Poisson"))
## The final null and full model for the OCF data is a genpois model#
null.OCF.Model <- glmmTMB(formula = null.OCF.form, data= OCF, family=genpois(link="log"), na.action = na.exclude)#
OCF.Model      <- glmmTMB(formula = beach.OCF.form, data= OCF, family=genpois(link="log"),#
                          control=glmmTMBControl(optimizer=optim, optArgs=list(method="BFGS")), na.action = na.exclude)#
summary(null.OCF.Model); sigma(null.OCF.Model) ## returns the dispersion parameter#
summary(OCF.Model); sigma(OCF.Model)#
#
## Checking the models#
simulationNULL <- simulateResiduals(fittedModel = null.OCF.Model, plot = F, n = 1000)#
simulationALL  <- simulateResiduals(fittedModel = OCF.Model, plot = F, n = 1000)#
plot(simulationNULL); plot(simulationALL)#
testDispersion(simulationOutput = simulationNULL, alternative ="less") #
testDispersion(simulationOutput = simulationALL, alternative ="less")
simulationNULL <- simulateResiduals(fittedModel = null.OCF.Model, plot = F, n = 1000)
simulationALL  <- simulateResiduals(fittedModel = OCF.Model, plot = F, n = 1000)
plot(simulationNULL); plot(simulationALL)
plot(simulationNULL);
## Final comparison#
ICtab(null.OCF.Model, OCF.Model,#
	  type = "AICc", delta = TRUE, base = TRUE, logLik= TRUE,#
	  weights=T, #
	  mnames=c("NULL", "~Beach"))
## Final comparison#
ICtab(null.OCF.Model, OCF.Model,#
	  type = "AICc", delta = TRUE, base = TRUE, logLik= TRUE,#
	  weights=T, #
	  mnames=c("NULL", "~Beach"))#
## Model with beach is best between the two#
#
## compute the evidence ratio#
evidence(aictab(cand.set = list(null.OCF.Model, OCF.Model), #
                modnames = c("NULL", "FULL"), #
                second.ord = TRUE))
## compare simulated data to the observed data #
simBEACH <- lapply(simulate(OCF.Model, seed = 1, nsim = 1000), function(x){ #
	cbind(x, OCF, dist="gen.pois")#
	})#
simBEACH           <- do.call(rbind, simBEACH)[,-3]#
OCF$dist           <- "observed"#
colnames(simBEACH) <- colnames(OCF[,c(2,1,3)]) #
ssd                <-  rbind(simBEACH, OCF)
## plot#
pal= c("#999999","#000000")#
ggplot(ssd, aes(x=OCF, fill= dist)) + facet_wrap(~Beach) + #
	  geom_histogram(aes(y = ..density..), position="dodge", binwidth=1, colour=1) +#
      xlab("Observed Clutch Frequency") + ylab(NULL) + #
      scale_fill_manual(values =pal) + #
      scale_x_continuous(limits = c(NA, 10)) +#
      theme_classic(base_size=20) + #
      theme(legend.title=element_blank(), #
            legend.position = "top")
## compute profile CI and mean differences#
confint(OCF.Model)#
emmeans(OCF.Model, ~ Beach, transform = "response")#
#
contrast(emmeans(OCF.Model, specs="Beach", transform = "response"),list(c(-1,1)))#
(exp(0.811) / exp(1.265) ) ## % underestimated is 0.635#
exp(0.811) - exp(1.265)    ## difference is on average 1.29
ull.ECF.Model <- glm(null.ECF.form, data = ECF, family=poisson(link="log"))#
ECF.Model      <- glm(beach.ECF.form, data = ECF, family=poisson(link="log"))#
#
## check dispersion#
testDispersion(null.ECF.Model); testDispersion(ECF.Model)#
#
## Overdispersoin not present so we do not need to fit a negative binomial #
## model to account for overdispersion but underdispersion is present so lets #
## use a gen.pois model in the glmmTMB package and do distribution selection#
ICtab(null.ECF.Model,#
	  glmmTMB(formula = null.ECF.form, data= ECF, family=genpois(link="log"), na.action = na.exclude), ## UPDATE THE MODEL DISTRIBUTIONS#
	  type = "AICc", delta = TRUE, base = TRUE, logLik= TRUE,#
	  weights=T, #
	  mnames=c("Poisson", "Generalized Poisson"))#
ICtab(ECF.Model,#
	  glmmTMB(formula = beach.ECF.form, data= ECF, family=genpois(link="log"),na.action = na.exclude,#
	          control=glmmTMBControl(optimizer=optim, optArgs=list(method="BFGS"))), ## UPDATE THE MODEL DISTRIBUTIONS#
	  type = "AICc", delta = TRUE, base = TRUE, logLik= TRUE,#
	  weights=T, #
	  mnames=c("Poisson", "Generalized Poisson"))
null.ECF.Model <- glm(null.ECF.form, data = ECF, family=poisson(link="log"))#
ECF.Model      <- glm(beach.ECF.form, data = ECF, family=poisson(link="log"))#
#
## check dispersion#
testDispersion(null.ECF.Model); testDispersion(ECF.Model)#
#
## Overdispersoin not present so we do not need to fit a negative binomial #
## model to account for overdispersion but underdispersion is present so lets #
## use a gen.pois model in the glmmTMB package and do distribution selection#
ICtab(null.ECF.Model,#
	  glmmTMB(formula = null.ECF.form, data= ECF, family=genpois(link="log"), na.action = na.exclude), ## UPDATE THE MODEL DISTRIBUTIONS#
	  type = "AICc", delta = TRUE, base = TRUE, logLik= TRUE,#
	  weights=T, #
	  mnames=c("Poisson", "Generalized Poisson"))#
ICtab(ECF.Model,#
	  glmmTMB(formula = beach.ECF.form, data= ECF, family=genpois(link="log"),na.action = na.exclude,#
	          control=glmmTMBControl(optimizer=optim, optArgs=list(method="BFGS"))), ## UPDATE THE MODEL DISTRIBUTIONS#
	  type = "AICc", delta = TRUE, base = TRUE, logLik= TRUE,#
	  weights=T, #
	  mnames=c("Poisson", "Generalized Poisson"))
## The final null and full model for the ECF data is a genpois model#
null.ECF.Model <- glmmTMB(formula = null.ECF.form, data= ECF, family=genpois(link="log"), na.action = na.exclude)#
ECF.Model      <- glmmTMB(formula = beach.ECF.form, data= ECF, family=genpois(link="log"),#
                          control=glmmTMBControl(optimizer=optim, optArgs=list(method="BFGS")), na.action = na.exclude)#
summary(null.ECF.Model); sigma(null.ECF.Model) ## returns the dispersion parameter#
summary(ECF.Model); sigma(ECF.Model)#
#
## Checking the models#
testDispersion(null.ECF.Model); #
testDispersion(ECF.Model); #
## underdispersion still present in the BEACH model
testDispersion(ECF.Model);
simulationNULL <- simulateResiduals(fittedModel = null.ECF.Model, plot = F, n = 1000)#
simulationALL  <- simulateResiduals(fittedModel = ECF.Model, plot = F, n = 1000)#
plot(simulationNULL)#
plot(simulationALL)
## Final comparison#
ICtab(null.ECF.Model, ECF.Model,#
	  type = "AICc", delta = TRUE, base = TRUE, logLik= TRUE,#
	  weights=T, #
	  mnames=c("NULL", "~Beach"))#
## Model with beach is best#
#
## compute the evidence ratio#
evidence(aictab(cand.set = list(null.ECF.Model, ECF.Model), #
                modnames = c("NULL", "FULL"), second.ord = TRUE))
## compare simulated data to the observed data #
simBEACH <- lapply(simulate(ECF.Model, seed = 1, nsim = 1000), function(x){ #
	cbind(x, ECF, dist="gen.pois")#
	})#
simBEACH           <- do.call(rbind, simBEACH)[,-3]#
ECF$dist           <- "observed"#
colnames(simBEACH) <- colnames(ECF[,c(2,1,3)]) #
ssd                <-  rbind(simBEACH, ECF)#
#
## plot#
ggplot(ssd, aes(x=ECF, fill= dist)) + facet_wrap(~Beach) + #
	  geom_histogram(aes(y = ..density..), position="dodge", binwidth=1, colour=1) +#
      xlab("Estimated Clutch Frequency") + ylab(NULL) + #
      scale_fill_manual(values =pal) + #
      scale_x_continuous(limits = c(NA, 10)) +#
      theme_classic(base_size=20) + #
      theme(legend.title=element_blank(), #
            legend.position = "top")
## compute profile CI and mean differences#
confint(ECF.Model)#
emmeans(ECF.Model, ~ Beach, transform = "response")#
contrast(emmeans(ECF.Model, specs="Beach", transform = "response"),list(c(-1,1)))#
(exp(0.864) / exp(1.364) ) ## % underestimated is 0.6065307#
exp(0.864) - exp(1.364)    ## difference is on average 1.539177#
#
#### -----------------------------------------
confint(ECF.Model)
emmeans(ECF.Model, ~ Beach, transform = "response")
contrast(emmeans(ECF.Model, specs="Beach", transform = "response"),list(c(-1,1)))
source("/Users/georgeglen/Documents/GitHub/Pfaller_et_al2021/functions.r")#
#
## Load each dataset in individually#
setwd("~/Documents/GitHub/Pfaller_et_al2021")#
FULL <- fun$importExcel(#
  fileName = "nesting_dat.xlsx",#
  sheetNames = c("Recruitment","OCF", "ECF","RMI","BF","Production"),#
  nskip = 0#
);#
#
## Seperate the datasets#
## 1: Clutch Frequency#
OCF <- FULL[[2]] %>% gather("Beach", "OCF")#
OCF$Beach <- str_replace_all(OCF$Beach, "OCF.", "")#
OCF <- OCF %>% transform(Beach = as.factor(Beach),#
					   OCF 	 = as.numeric(OCF)) %>%#
			 filter(!is.na(OCF))#
#
ECF <- FULL[[3]] %>% gather("Beach", "ECF")#
ECF$Beach <- str_replace_all(ECF$Beach, "ECF.", "")#
ECF <- ECF %>% transform(Beach = as.factor(Beach),#
					   ECF 	 = as.numeric(ECF)) %>%#
			 filter(!is.na(ECF))#
#
## 2: RMI#
RMI <- FULL[[4]] %>% gather("Beach", "RMI")#
RMI$Beach <- str_replace_all(RMI$Beach, "RI.", "")#
RMI <- RMI %>% transform(Beach = as.factor(Beach),#
						 RMI 	 = as.numeric(RMI)) %>%#
			   filter(!is.na(RMI))#
## 3: BP#
BF <- FULL[[5]] %>% gather("Beach", "BF")#
BF$Beach <- str_replace_all(BF$Beach, "BF.", "")#
BF <- BF %>% transform(Beach = as.factor(Beach),#
					   BF 	 = as.numeric(BF)) %>%#
			   filter(!is.na(BF))#
#
#### ----------------------------------------- #####
#### Models#
#### ----------------------------------------- #####
#
## Write out the model formulas#
## 1: OCF#
null.OCF.form <- formula(OCF ~ 1)#
beach.OCF.form <- formula(OCF ~ Beach - 1)#
#
null.ECF.form <- formula(ECF ~ 1)#
beach.ECF.form <- formula(ECF ~ Beach - 1)#
#
## 2: RMI #
null.RMI.form <- formula(RMI ~ 1)#
beach.RMI.form <- formula(RMI ~ Beach - 1)#
#
## 3: BP#
null.BF.form <- formula(BF ~ 1)#
beach.BF.form <- formula(BF ~ Beach - 1)#
#
#### ----------------------------------------- #####
#### Model for OCF#
#### ----------------------------------------- #####
#
## Create initial models#
null.OCF.Model <- glm(null.OCF.form, data = OCF, family=poisson(link="log"))#
OCF.Model      <- glm(beach.OCF.form, data = OCF, family=poisson(link="log"))#
#
## check dispersion#
testDispersion(null.OCF.Model); testDispersion(OCF.Model)#
#
## Overdispersoin not present so we do not need to fit a negative binomial #
## model to account for overdispersion but underdispersion is present so lets #
## use a gen.pois model in the glmmTMB package and do distribution selection#
#
ICtab(null.OCF.Model,#
	  glmmTMB(formula = null.OCF.form, data= OCF, family=genpois(link="log"), na.action = na.exclude), ## UPDATE THE MODEL DISTRIBUTIONS#
	  type = "AICc", delta = TRUE, base = TRUE, logLik= TRUE,#
	  weights=T, #
	  mnames=c("Poisson", "Generalized Poisson"))#
ICtab(OCF.Model,#
	  glmmTMB(formula = beach.OCF.form, data= OCF, family=genpois(link="log"),na.action = na.exclude,#
	          control=glmmTMBControl(optimizer=optim, optArgs=list(method="BFGS"))), ## UPDATE THE MODEL DISTRIBUTIONS#
	  type = "AICc", delta = TRUE, base = TRUE, logLik= TRUE,#
	  weights=T, #
	  mnames=c("Poisson", "Generalized Poisson"))	  #
## The final null and full model for the OCF data is a genpois model#
null.OCF.Model <- glmmTMB(formula = null.OCF.form, data= OCF, family=genpois(link="log"), na.action = na.exclude)#
OCF.Model      <- glmmTMB(formula = beach.OCF.form, data= OCF, family=genpois(link="log"),#
                          control=glmmTMBControl(optimizer=optim, optArgs=list(method="BFGS")), na.action = na.exclude)#
summary(null.OCF.Model); sigma(null.OCF.Model) ## returns the dispersion parameter#
summary(OCF.Model); sigma(OCF.Model)#
#
## Checking the models#
simulationNULL <- simulateResiduals(fittedModel = null.OCF.Model, plot = F, n = 1000)#
simulationALL  <- simulateResiduals(fittedModel = OCF.Model, plot = F, n = 1000)#
plot(simulationNULL); plot(simulationALL)#
testDispersion(simulationOutput = simulationNULL, alternative ="less") #
testDispersion(simulationOutput = simulationALL, alternative ="less") #
#
## Final comparison#
ICtab(null.OCF.Model, OCF.Model,#
	  type = "AICc", delta = TRUE, base = TRUE, logLik= TRUE,#
	  weights=T, #
	  mnames=c("NULL", "~Beach"))#
## Model with beach is best between the two#
#
## compute the evidence ratio#
evidence(aictab(cand.set = list(null.OCF.Model, OCF.Model), #
                modnames = c("NULL", "FULL"), #
                second.ord = TRUE))#
#
## compare simulated data to the observed data #
simBEACH <- lapply(simulate(OCF.Model, seed = 1, nsim = 1000), function(x){ #
	cbind(x, OCF, dist="gen.pois")#
	})#
simBEACH           <- do.call(rbind, simBEACH)[,-3]#
OCF$dist           <- "observed"#
colnames(simBEACH) <- colnames(OCF[,c(2,1,3)]) #
ssd                <-  rbind(simBEACH, OCF)#
#
## plot#
pal= c("#999999","#000000")#
ggplot(ssd, aes(x=OCF, fill= dist)) + facet_wrap(~Beach) + #
	  geom_histogram(aes(y = ..density..), position="dodge", binwidth=1, colour=1) +#
      xlab("Observed Clutch Frequency") + ylab(NULL) + #
      scale_fill_manual(values =pal) + #
      scale_x_continuous(limits = c(NA, 10)) +#
      theme_classic(base_size=20) + #
      theme(legend.title=element_blank(), #
            legend.position = "top")#
#
## compute profile CI and mean differences#
confint(OCF.Model)#
emmeans(OCF.Model, ~ Beach, transform = "response")#
#
contrast(emmeans(OCF.Model, specs="Beach", transform = "response"),list(c(-1,1)))#
(exp(0.811) / exp(1.265) ) ## % underestimated is 0.635#
exp(0.811) - exp(1.265)    ## difference is on average 1.29
contrast(emmeans(OCF.Model, specs="Beach", transform = "response"),list(c(-1,1)))
exp(0.811) - exp(1.265)    ## difference is on average 1.29
emmeans(OCF.Model, ~ Beach, transform = "response")
OCF.Model
coef(OCF.Model)
fixef(OCF.Model)
fixef(OCF.Model)$cond
fixef(OCF.Model)$cond[1]
(exp(fixef(OCF.Model)$cond[2]) / exp(fixef(OCF.Model)$cond[1]) ) ## % underestimated is 0.635
exp(fixef(OCF.Model)$cond[2]) - exp(fixef(OCF.Model)$cond[1])    ## difference is on average 1.29
null.RMI.Model <- glm(null.RMI.form, data = RMI, family=poisson(link="log"))#
RMI.Model      <- glm(beach.RMI.form, data = RMI, family=poisson(link="log"))#
#
## check dispersion#
testDispersion(null.RMI.Model); testDispersion(RMI.Model)#
#
## Overdispersoin not present so we do not need to fit a negative binomial #
## model to account for overdispersion but underdispersion is present so lets #
## use a gen.pois model in the glmmTMB package and do distribution selection#
ICtab(null.RMI.Model,#
	  glmmTMB(formula = null.RMI.form, data= RMI, family=genpois(link="log"), na.action = na.exclude, #
	          control=glmmTMBControl(optimizer=optim, optArgs=list(method="BFGS"))), ## UPDATE THE MODEL DISTRIBUTIONS#
	  type = "AICc", delta = TRUE, base = TRUE, logLik= TRUE,#
	  weights=T, #
	  mnames=c("Poisson", "Generalized Poisson"))#
ICtab(RMI.Model,#
	  glmmTMB(formula = beach.RMI.form, data= RMI, family=genpois(link="log"), na.action = na.exclude,#
	          control=glmmTMBControl(optimizer=optim, optArgs=list(method="BFGS"))), ## UPDATE THE MODEL DISTRIBUTIONS#
	  type = "AICc", delta = TRUE, base = TRUE, logLik= TRUE,#
	  weights=T, #
	  mnames=c("Poisson", "Generalized Poisson"))	  #
## The final null and full model for the OCF data is a genpois model#
null.RMI.Model <- glmmTMB(formula = null.RMI.form, data= RMI, family=genpois(link="log"),#
                          control=glmmTMBControl(optimizer=optim, optArgs=list(method="BFGS")),#
                          na.action = na.exclude)#
RMI.Model      <-  glmmTMB(formula = beach.RMI.form, data= RMI, family=genpois(link="log"),#
						   control=glmmTMBControl(optimizer=optim, optArgs=list(method="BFGS")), na.action = na.exclude)#
summary(null.RMI.Model); sigma(null.RMI.Model) ## returns the dispersion parameter#
summary(RMI.Model); sigma(RMI.Model) ## returns the dispersion parameter
## Checking the models#
testDispersion(null.RMI.Model); #
testDispersion(RMI.Model )#
simulationNULL <- simulateResiduals(fittedModel = null.RMI.Model, plot = F)#
simulationALL  <- simulateResiduals(fittedModel = RMI.Model , plot = F)#
plot(simulationNULL); #
plot(simulationALL)
## Final comparison#
ICtab(null.RMI.Model, RMI.Model ,#
	  type = "AICc", delta = TRUE, base = TRUE, logLik= TRUE,#
	  weights=T, #
	  mnames=c("NULL", "~Beach"))#
## NULL model is best#
#
## compute the evidence ratio#
evidence(aictab(cand.set = list(null.RMI.Model, RMI.Model), #
                modnames = c("NULL", "FULL"), second.ord = TRUE))#
#
## compare simulated data to the observed data #
simBEACH <- lapply(simulate(RMI.Model, seed = 1, nsim = 1000), function(x){ #
	cbind(x, RMI, dist="gen.pois")#
	})#
simBEACH           <- do.call(rbind, simBEACH)[,-3]#
RMI$dist           <- "observed"#
colnames(simBEACH) <- colnames(RMI[,c(2,1,3)]) #
ssd                <-  rbind(simBEACH, RMI)#
#
## plot#
ggplot(ssd, aes(x=RMI, fill= dist)) + facet_wrap(~Beach) + #
	  geom_histogram(aes(y = ..density..), position="dodge", binwidth=1, colour=1) +#
      xlab("Remigration Interval") + ylab(NULL) + #
      scale_fill_manual(values =pal) + #
      scale_x_continuous(limits = c(NA, 10)) +#
      theme_classic(base_size=20) + #
      theme(legend.title=element_blank(), #
            legend.position = "top")
confint(RMI.Model)
emmeans(RMI.Model, ~ 1, transform = "response")
emmeans(RMI.Model, ~ Beach, transform = "response")
contrast(emmeans(RMI.Model, specs="Beach", transform = "response"),list(c(-1,1)))
(exp(fixef(ECF.Model)$cond[2]) / exp(fixef(ECF.Model)$cond[1]) ) ## % underestimated is 0.6065307
exp(fixef(ECF.Model)$cond[2]) - exp(fixef(ECF.Model)$cond[1])    ## difference is on average 1.539177
(exp(fixef(RMI.Model)$cond[2]) / exp(fixef(RMI.Model)$cond[1]) ) ## % underestimated is 0.6065307
exp(fixef(RMI.Model)$cond[2]) - exp(fixef(RMI.Model)$cond[1])    ## difference is on average 1.539177
#### ----------------------------------------- #####
#### Model for BP#
#### ----------------------------------------- #####
#### The number of breeding seasons for turtles across the 8-year study period#
#
null.BF.Model <- glm(null.BF.form, data = BF, family=poisson(link="log"))#
BF.Model      <- glm(beach.BF.form, data = BF, family=poisson(link="log"))#
#
## check dispersion#
testDispersion(null.BF.Model); testDispersion(BF.Model )#
#
## Overdispersoin not present so we do not need to fit a negative binomial #
## model to account for overdispersion but underdispersion is present so lets #
## use a gen.pois model in the glmmTMB package and do distribution selection#
ICtab(null.BF.Model,#
	  glmmTMB(formula = null.BF.form, data= BF, family=genpois(link="log"), na.action = na.exclude,#
	          control=glmmTMBControl(optimizer=optim, optArgs=list(method="BFGS"))), ## UPDATE THE MODEL DISTRIBUTIONS#
	  type = "AICc", delta = TRUE, base = TRUE, logLik= TRUE,#
	  weights=T, #
	  mnames=c("Poisson", "Generalized Poisson"))#
ICtab(BF.Model,#
	  glmmTMB(formula = beach.BF.form, data= BF, family=genpois(link="log"), na.action = na.exclude,#
	          control=glmmTMBControl(optimizer=optim, optArgs=list(method="BFGS"))), ## UPDATE THE MODEL DISTRIBUTIONS#
	  type = "AICc", delta = TRUE, base = TRUE, logLik= TRUE,#
	  weights=T, #
	  mnames=c("Poisson", "Generalized Poisson"))
## The final null and full model for the OCF data is a genpois model#
null.BF.Model <- glmmTMB(formula = null.BF.form, data= BF, family=genpois(link="log"), na.action = na.exclude,#
                          control=glmmTMBControl(optimizer=optim, optArgs=list(method="BFGS")))#
#
BF.Model      <-  glmmTMB(formula = beach.BF.form, data= BF, family=genpois(link="log"),  na.action = na.exclude,#
                          control=glmmTMBControl(optimizer=optim, optArgs=list(method="BFGS")))#
summary(null.BF.Model); sigma(null.BF.Model) ## dispersion parameter#
summary(BF.Model); sigma(BF.Model) ## dispersion parameter#
#
## by hand, the dispersion parameter is computed as #
exp(BF.Model$sdr$par.fixed[3])
summary(BF.Model); sigma(BF.Model) ## dispersion parameter
## Checking the models#
testDispersion(null.BF.Model); #
testDispersion(BF.Model )#
simulationNULL <- simulateResiduals(fittedModel = null.BF.Model, plot = F)#
simulationALL  <- simulateResiduals(fittedModel = BF.Model , plot = F)#
plot(simulationNULL); #
plot(simulationALL)
## Final comparison#
ICtab(null.BF.Model, BF.Model ,#
	  type = "AICc", delta = TRUE, base = TRUE, logLik= TRUE,#
	  weights=T, #
	  mnames=c("NULL", "~Beach"))#
## FULL model is best#
#
## compute the evidence ratio#
evidence(aictab(cand.set = list(null.BF.Model, BF.Model), #
                modnames = c("NULL", "FULL"), second.ord = TRUE))
## compare simulated data to the observed data #
simBEACH <- lapply(simulate(BF.Model, seed = 1, nsim = 1000), function(x){ #
	cbind(x, BF, dist="gen.pois")#
	})#
simBEACH           <- do.call(rbind, simBEACH)[,-3]#
BF$dist           <- "observed"#
colnames(simBEACH) <- colnames(BF[,c(2,1,3)]) #
ssd                <-  rbind(simBEACH, BF)#
#
## plot#
ggplot(ssd, aes(x=BF, fill= dist)) + facet_wrap(~Beach) + #
	  geom_histogram(aes(y = ..density..), position="dodge", binwidth=1, colour=1) +#
      xlab("Breeding Frequency") + ylab(NULL) + #
      scale_fill_manual(values =pal) + #
      scale_x_continuous(limits = c(NA, 10)) +#
      theme_classic(base_size=20) + #
      theme(legend.title=element_blank(), #
            legend.position = "top")
(exp(fixef(BF.Model)$cond[2]) / exp(fixef(BF.Model)$cond[1]) ) ## % underestimated is 1.002551
(exp(fixef(BF.Model)$cond[2])
exp(fixef(BF.Model)
)
exp(fixef(BF.Model)$cond[1])
exp(fixef(BF.Model)$cond[2]) - exp(fixef(BF.Model)$cond[1])    ## difference is on average 0.007159968
RECRUIT <- FULL[[1]] %>% gather("Beach", "NestingStatus", -Year)#
RECRUIT$Beach <- str_replace_all(RECRUIT$Beach, "Status.", "")#
RECRUIT <- RECRUIT %>% transform(Year 			= as.numeric(Year), #
                                 Beach 			= as.factor(Beach),#
                                 NestingStatus 	= as.factor(NestingStatus))#
RECRUIT <- as.data.frame(table(RECRUIT[,1:3])) %>% filter(Freq>0)#
## FOCAL = FALSE + TRUE#
## ALL = NEO#
## FOCAL overestimate proportion of neophytes#
#
## count the number of nesters annually#
nester.counter <- function(dat){#
  by(dat, dat$Year, FUN=function(x){#
  ##FILTER BY BEACH #
  beach.filt <- x[which(x$Beach == "All"),]#
  ## SUM BY Year#
  dat <- sum(beach.filt$Freq)#
  return(dat)#
  })#
}
empAcounts <- data.frame(Year = as.factor(seq(2013, 2017, 1)),#
                          NN   = apply(nester.counter(RECRUIT), 1, unlist))#
RECRUIT <- left_join(RECRUIT, tempAcounts)#
#
## Count the number of recruits on the focal nestinng beach#
recruit.counter <- function(dat){#
  by(dat, dat$Year, FUN=function(x){#
  ##FILTER BY BEACH #
  beach.filt <- x[which(x$Beach == "Focal" & x$NestingStatus != "REM"), -5]#
  ## SUM BY Year#
  dat <- sum(beach.filt$Freq)#
  return(dat)#
  })#
}#
tempRcounts <- data.frame(Year  = as.factor(seq(2013, 2017, 1)),#
                          Beach = "Focal",#
                          REC   = apply(recruit.counter(RECRUIT), 1, unlist))
RECRUIT     <- left_join(RECRUIT, tempRcounts)#
RECRUIT$REC <- ifelse(is.na(RECRUIT$REC), RECRUIT$Freq, RECRUIT$RE) #
Mod.dat <- RECRUIT %>% filter(NestingStatus != "REM") %>% dplyr::select(Year, Beach, REC, NN) %>% unique()#
Mod.dat$Year <- as.numeric(as.character(Mod.dat$Year))#
Mod.dat$Beach <- as.factor(Mod.dat$Beach)
tempAcounts <- data.frame(Year = as.factor(seq(2013, 2017, 1)),#
                          NN   = apply(nester.counter(RECRUIT), 1, unlist))#
RECRUIT <- left_join(RECRUIT, tempAcounts)#
#
## Count the number of recruits on the focal nestinng beach#
recruit.counter <- function(dat){#
  by(dat, dat$Year, FUN=function(x){#
  ##FILTER BY BEACH #
  beach.filt <- x[which(x$Beach == "Focal" & x$NestingStatus != "REM"), -5]#
  ## SUM BY Year#
  dat <- sum(beach.filt$Freq)#
  return(dat)#
  })#
}#
tempRcounts <- data.frame(Year  = as.factor(seq(2013, 2017, 1)),#
                          Beach = "Focal",#
                          REC   = apply(recruit.counter(RECRUIT), 1, unlist))#
RECRUIT     <- left_join(RECRUIT, tempRcounts)#
RECRUIT$REC <- ifelse(is.na(RECRUIT$REC), RECRUIT$Freq, RECRUIT$RE) #
Mod.dat <- RECRUIT %>% filter(NestingStatus != "REM") %>% dplyr::select(Year, Beach, REC, NN) %>% unique()#
Mod.dat$Year <- as.numeric(as.character(Mod.dat$Year))#
Mod.dat$Beach <- as.factor(Mod.dat$Beach)
null.REC.form <- formula(cbind(REC, NN-REC) ~ 1)
FULL.REC.form <- formula(cbind(REC, NN-REC) ~ Beach -1)
## write out model formulas#
null.REC.form <- formula(cbind(REC, NN-REC) ~ 1)#
FULL.REC.form <- formula(cbind(REC, NN-REC) ~ Beach -1)#
#
## fit binomial models#
options(na.action = "na.fail") #
null.REC.Model <- glm(null.REC.form, data = Mod.dat, family=binomial(link="logit"))#
REC.Model      <- glm(FULL.REC.form, data = Mod.dat, family=binomial(link="logit"))#
#
summary(null.REC.Model);#
summary(REC.Model);  #
#
## Final comparison#
ICtab(null.REC.Model, REC.Model,#
	  type = "AICc", delta = TRUE, base = TRUE, logLik= TRUE,#
	  weights=T, #
	  mnames=c("NULL", "~Beach"))#
## beach model is best out of the set#
#
## compute the evidence ratio#
evidence(aictab(cand.set = list(null.REC.Model, REC.Model), #
                modnames = c("NULL", "RECRUIT"), second.ord = TRUE))
## compute profile CI and mean differences#
confint(REC.Model)#
emmeans(REC.Model, ~ 1, transform = "response")#
emmeans(REC.Model, specs= c("Beach"), transform = "response")#
contrast(emmeans(REC.Model, specs="Beach", transform = "response"), list(c(-1,1)))
(exp(fixef(REC.Model)$cond[2]) / exp(fixef(REC.Model)$cond[1]) ) ## % underestimated is 0.6640397
(exp(coef(REC.Model)$cond[2]) / exp(coef(REC.Model)$cond[1]) ) ## % underestimated is 0.6640397
contrast(emmeans(REC.Model, specs="Beach", transform = "response"), list(c(-1,1)))
tab_model(null.BF.Model,  BF.Model,#
          show.stat = T, show.se = T, show.dev = T, show.loglik = T, #
          show.icc =T, show.ngroups = T, show.re.var = T, p.style = "numeric",#
          CSS = css_theme("regression"), collapse.ci = T, #transform = NULL,#
          #bootstrap = T, seed = 123, iterations = 1000,#
          vcov.type = "HC0", digits = 2, p.val = "wald")#
tab_model(null.RMI.Model, RMI.Model,#
          show.stat = T, show.se = T, show.dev = T, show.loglik = T, #
          show.icc =T, show.ngroups = T, show.re.var = T, p.style = "numeric",#
          CSS = css_theme("regression"), collapse.ci = T, #transform = NULL,#
          #bootstrap = T, seed = 123, iterations = 1000,#
          vcov.type = "HC0", digits = 2, p.val = "wald")#
tab_model(null.OCF.Model,  OCF.Model, #
          show.stat = T, show.se = T, show.dev = T, show.loglik = T, #
          show.icc =T, show.ngroups = T, show.re.var = T, p.style = "numeric",#
          CSS = css_theme("regression"), collapse.ci = T, #transform = NULL,#
          #bootstrap = T, seed = 123, iterations = 1000,#
          vcov.type = "HC0", digits = 2, p.val = "wald")#
tab_model(null.ECF.Model,  ECF.Model, #
          show.stat = T, show.se = T, show.dev = T, show.loglik = T, #
          show.icc =T, show.ngroups = T, show.re.var = T, p.style = "numeric",#
          CSS = css_theme("regression"), collapse.ci = T, #transform = NULL,#
          #bootstrap = T, seed = 123, iterations = 1000,#
          vcov.type = "HC0", digits = 2, p.val = "wald")#
tab_model(null.REC.Model,  REC.Model, #
          show.stat = T, show.se = T, show.dev = T, show.loglik = T, #
          show.icc =T, show.ngroups = T, show.re.var = T, p.style = "numeric",#
          CSS = css_theme("regression"), collapse.ci = T, #transform = NULL,#
          #bootstrap = T, seed = 123, iterations = 1000,#
          vcov.type = "HC0", digits = 2, p.val = "wald")
## Plot model coefficients#
## count models#
plot_summs(OCF.Model, ECF.Model, RMI.Model, BF.Model, #
            ci_level = 0.95, scale = TRUE, robust = TRUE, #
             confint = TRUE, digits = 3, inner_ci_level = .9, exp=T, grid = TRUE, #
            model.names = c("OCF Model w/ BEACH","ECF Model w/ BEACH",#
                             "RMI Model w/ BEACH","BF Model w/ BEACH")) +#
             scale_x_continuous(breaks=seq(1, 5, 0.3)) +#
             theme_sjplot2() + labs(x="Parameter Estimates", y="") + theme(text = element_text(size=25))
## proportion model#
plot_summs(REC.Model, ci_level = 0.95, robust = TRUE, plot.distributions=T,#
            confint = TRUE, digits = 3, inner_ci_level = .9, grid = TRUE) +#
           theme_sjplot2() + labs(x="Parameter Estimates", y="") + theme(text = element_text(size=25))
## proportion model#
plot_summs(REC.Model, ci_level = 0.95, robust = TRUE, plot.distributions=T,#
            confint = TRUE, digits = 3, inner_ci_level = .9, grid = TRUE) +#
           theme_sjplot2() + labs(x="Parameter Estimates", y="") + theme(text = element_text(size=25))#
#
# Compute disgnostic stats (square root of the variance of the residuals) for each model #
summaryALL <- data.frame(models=c("NULL OCF","OCF w/ BEACH",#
                                  "NULL ECF","ECF w/ BEACH",#
                                  "NULL RMI","RMI w/ BEACH",#
                                  "NULL BF","BF w/ BEACH", #
                                  "NULL REC","REC w/ BEACH"), #
                         rmse_value = NA, AICc_value = NA, Deviance_Exp = NA, RRS_value = NA)#
mod_names <- list(null.OCF.Model, OCF.Model, null.ECF.Model, ECF.Model, #
                  null.RMI.Model, RMI.Model, null.BF.Model, BF.Model, #
                  null.REC.Model, REC.Model)#
library(ie2misc); library(qpcR)#
for (i in seq_along(mod_names)) {#
  summaryALL$rmse_value[i]   <- sjstats::rmse(mod_names[[i]])#
  summaryALL$AICc_value[i]   <- AICc(mod_names[[i]])#
  summaryALL$Deviance_Exp[i] <- AICc(mod_names[[i]])#
  summaryALL$RRS_value[i]    <- sum(residuals(mod_names[[i]])^2) #RSS(mod_names[[i]]) # Weighted sum of squares#
}#
summaryALL$delta_AICc <- c(#
		ICtab(null.OCF.Model, OCF.Model, k=2, sort=F, base=TRUE, weights=TRUE, logLik=TRUE,  type = "AICc")$dAICc,#
		ICtab(null.ECF.Model, ECF.Model, k=2, sort=F, base=TRUE, weights=TRUE, logLik=TRUE, type = "AICc")$dAICc, #
		ICtab(null.RMI.Model, RMI.Model, k=2, sort=F, base=TRUE, weights=TRUE, logLik=TRUE,  type = "AICc")$dAICc,#
		ICtab(null.BF.Model, BF.Model, k=2, sort=F, base=TRUE, weights=TRUE, logLik=TRUE,  type = "AICc")$dAICc,#
		ICtab(null.REC.Model, REC.Model, k=2, sort=F, base=TRUE, weights=TRUE, logLik=TRUE,  type = "AICc")$dAICc)#
summaryALL
q()