.Rhistory

MDSplot(rf_classifier)
?MDSplot
MDSplot(rf_classifier, (as.factor(ran.forest.data$maturity))
MDSplot(rf_classifier, as.factor(ran.forest.data$maturity))
MDSplot(rf_classifier, as.factor(ran.forest.data$maturity))
rf_classifier <- randomForest(as.factor(maturity) ~ ccl + ovilen + ovavol,
data = ran.forest.data, ntree=10000, importance=TRUE, proximity=TRUE)
varImp(rf_classifier, scale = FALSE)
varImpPlot(rf_classifier,
sort = T,
main = "Top 3 - Variable Importance")
plot(rf_classifier)
MDSplot(rf_classifier, as.factor(ran.forest.data$maturity))
?MDSplot
MDSplot(rf_classifier, as.factor(ran.forest.data$maturity), k=3)
MDSplot(rf_classifier, as.factor(ran.forest.data$maturity), k=2)
MDSplot(rf_classifier, as.factor(ran.forest.data$maturity), k=1)
MDSplot(rf_classifier, as.factor(ran.forest.data$maturity), k=21)
MDSplot(rf_classifier, as.factor(ran.forest.data$maturity), k=2)
MDSplot(rf_classifier, as.factor(ran.forest.data$maturity), k=3)
MDSplot(rf_classifier, as.factor(ran.forest.data$maturity), k=2)
?VSURF
mat.vsurf <- VSURF(as.factor(maturity) ~., data = ran.forest.data,
ntree = 500, nfor.thres = 20,
nfor.interp = 10, nfor.pred = 10)
summary(mat.vsurf)
library(VSURF)
varImpPlot(rf_classifier,
sort = T,
main = "Top 3 - Variable Importance")
plot(rf_classifier)
varImp(rf_classifier, scale = FALSE)
varImpPlot(rf_classifier,
sort = T,
main = "Top 3 - Variable Importance")
varImp(rf_classifier, scale = FALSE)
varImp(rf_classifier, scale = T)
varImp(rf_classifier, scale = FALSE)
varImpPlot(rf_classifier,
sort = T, scale = FALSE,
main = "Top 3 - Variable Importance")
varImpPlot(rf_classifier,
sort = T, scale = T,
main = "Top 3 - Variable Importance")
importance(rf_classifier)
library(lava)
lvm
?lvm
regression
?regression
colnames(BRED_Data)
regression(semMod) <- list(
maturity ~ ccl + ovilen + ovavol
ovilen ~ ovavol,
ovilen ~ ovavol)
?regression
regression(semMod) <- list(
maturity ~ ccl + ovilen + ovavol,
ovilen ~ ovavol,
ovilen ~ ovavol,
ccl    ~ ovilen + ovavol)
list(
maturity ~ ccl + ovilen + ovavol,
ovilen ~ ovavol,
ovilen ~ ovavol,
ccl    ~ ovilen + ovavol)
# This is a package lava SEM
semMod <- lvm(){
regression(semMod) <- list(
maturity ~ ccl + ovilen + ovavol,
ovilen ~ ovavol,
ccl    ~ ovilen
ccl    ~ ovavol)
latent(semMod) <- ~ ccl + ovilen + ovavol
regression(semMod, ovilen ~ ovavol) <- 1
regression(semMod, ccl ~ ovilen) <- 1
regression(semMod, ccl ~ ovavol) <- 1
covariance(semMod) <- ovilen ~ ovavol # error correlations between latents
covariance(semMod) <- ccl ~ ovavol # error correlations between latents
covariance(semMod) <- ccl ~ ovilen # error correlations between latents
}
latent(semMod) <- ~ ccl + ovilen + ovavol
# This is a package lava SEM
semMod <- lvm(){
regression(semMod) <- list(
maturity ~ ccl + ovilen + ovavol,
ovilen ~ ovavol,
ccl    ~ ovilen
ccl    ~ ovavol)
latent(semMod) <- ~ ccl + ovilen + ovavol
regression(semMod, ovilen ~ ovavol) <- 1
regression(semMod, ccl ~ ovilen) <- 1
regression(semMod, ccl ~ ovavol) <- 1
covariance(semMod) <- ovilen ~ ovavol # error correlations between latents
covariance(semMod) <- ccl ~ ovavol # error correlations between latents
covariance(semMod) <- ccl ~ ovilen # error correlations between latents
}
regression(semMod) <- list(
maturity ~ ccl + ovilen + ovavol,
ovilen ~ ovavol,
ccl    ~ ovilen
ccl    ~ ovavol)
# This is a package lava SEM
semMod <- lvm(){
regression(semMod) <- list(
maturity ~ ccl + ovilen + ovavol,
ovilen ~ ovavol,
ccl    ~ ovilen,
ccl    ~ ovavol)
latent(semMod) <- ~ ccl + ovilen + ovavol
regression(semMod, ovilen ~ ovavol) <- 1
regression(semMod, ccl ~ ovilen) <- 1
regression(semMod, ccl ~ ovavol) <- 1
covariance(semMod) <- ovilen ~ ovavol # error correlations between latents
covariance(semMod) <- ccl ~ ovavol # error correlations between latents
covariance(semMod) <- ccl ~ ovilen # error correlations between latents
}
lvm()
# This is a package lava SEM
semMod <- lvm()
semMod <- lvm(){
regression(semMod) <- list(
maturity ~ ccl + ovilen + ovavol,
ovilen ~ ovavol,
ccl    ~ ovilen,
ccl    ~ ovavol)
latent(semMod) <- ~ ccl + ovilen + ovavol
regression(semMod, ovilen ~ ovavol) <- 1
regression(semMod, ccl ~ ovilen) <- 1
regression(semMod, ccl ~ ovavol) <- 1
covariance(semMod) <- ovilen ~ ovavol # error correlations between latents
covariance(semMod) <- ccl ~ ovavol # error correlations between latents
covariance(semMod) <- ccl ~ ovilen # error correlations between latents
}
regression(semMod) <- list(
maturity ~ ccl + ovilen + ovavol,
ovilen ~ ovavol,
ccl    ~ ovilen,
ccl    ~ ovavol)
latent(semMod) <- ~ ccl + ovilen + ovavol
regression(semMod, ovilen ~ ovavol) <- 1
regression(semMod, ccl ~ ovilen) <- 1
regression(semMod, ccl ~ ovavol) <- 1
covariance(semMod) <- ovilen ~ ovavol # error correlations between latents
covariance(semMod) <- ccl ~ ovavol # error correlations between latents
covariance(semMod) <- ccl ~ ovilen # error correlations between latents
semMod
# This is a package lava SEM
semMod <- lvm()
regression(semMod) <- list(
maturity ~ ccl + ovilen + ovavol,
ovilen ~ ovavol,
ccl    ~ ovilen,
ccl    ~ ovavol)
latent(semMod) <- ~ ccl + ovilen + ovavol
regression(semMod, ovilen ~ ovavol) <- 1
regression(semMod, ccl ~ ovilen) <- 1
regression(semMod, ccl ~ ovavol) <- 1
covariance(semMod) <- ovilen ~ ovavol # error correlations between latents
covariance(semMod) <- ccl ~ ovavol # error correlations between latents
covariance(semMod) <- ccl ~ ovilen # error correlations between latents
semMod
plot(semMod)
estimate(semMod, BRED_Data)
lava.options(param = "hybrid")
estimate(semMod, BRED_Data)
# This is a package lava SEM
semMod <- binomial.lvm()
regression(semMod) <- list(
maturity ~ ccl + ovilen + ovavol,
ovilen ~ ovavol,
ccl    ~ ovilen,
ccl    ~ ovavol)
regression(semMod) <- list(
maturity ~ ccl + ovilen + ovavol,
ovilen ~ ovavol,
ccl    ~ ovilen,
ccl    ~ ovavol)?binomial.lvm
?binomial.lvm
estimate(semMod, BRED_Data, estimator="glm")
# This is a package lava SEM
semMod <- lvm()
regression(semMod) <- list(
maturity ~ ccl + ovilen + ovavol,
ovilen ~ ovavol,
ccl    ~ ovilen,
ccl    ~ ovavol)
latent(semMod) <- ~ ccl + ovilen + ovavol
regression(semMod, ovilen ~ ovavol) <- 1
regression(semMod, ccl ~ ovilen) <- 1
regression(semMod, ccl ~ ovavol) <- 1
distribution(semMod, ~ ccl + ovilen + ovavol) <- binomial.lvm("logit")
lava.options(param = "hybrid")
estimate(semMod, BRED_Data, estimator="glm")
distribution(semMod, maturity ~ ccl + ovilen + ovavol) <- binomial.lvm("logit")
# This is a package lava SEM
semMod <- lvm()
regression(semMod) <- list(
maturity ~ ccl + ovilen + ovavol,
ovilen ~ ovavol,
ccl    ~ ovilen,
ccl    ~ ovavol)
latent(semMod) <- ~ ccl + ovilen + ovavol
regression(semMod, ovilen ~ ovavol) <- 1
regression(semMod, ccl ~ ovilen) <- 1
regression(semMod, ccl ~ ovavol) <- 1
covariance(semMod) <- ovilen ~ ovavol # error correlations between latents
covariance(semMod) <- ccl ~ ovavol # error correlations between latents
covariance(semMod) <- ccl ~ ovilen # error correlations between latents
distribution(semMod, maturity ~ ccl + ovilen + ovavol) <- binomial.lvm("logit")
semMod
distribution(semMod, ccl ~ ovilen + ovavol) <- gaussian.lvm("log")
semMod
lava.options(param = "hybrid")
estimate(semMod, BRED_Data, estimator="glm")
# This is a package lava SEM
semMod <- lvm()
regression(semMod) <- list(
maturity ~ ccl + ovilen + ovavol,
ovilen ~ ovavol,
ccl    ~ ovilen,
ccl    ~ ovavol)
latent(semMod) <- ~ ccl + ovilen + ovavol
regression(semMod, ovilen ~ ovavol) <- 1
regression(semMod, ccl ~ ovilen) <- 1
regression(semMod, ccl ~ ovavol) <- 1
#covariance(semMod) <- ovilen ~ ovavol # error correlations between latents
#covariance(semMod) <- ccl ~ ovavol # error correlations between latents
#covariance(semMod) <- ccl ~ ovilen # error correlations between latents
distribution(semMod, maturity ~ ccl + ovilen + ovavol) <- binomial.lvm("logit")
distribution(semMod, ccl ~ ovilen + ovavol) <- gaussian.lvm("log")
lava.options(param = "hybrid")
estimate(semMod, BRED_Data, estimator="glm")
semMod
# This is a package lava SEM
semMod <- lvm()
regression(semMod) <- list(
maturity ~ ccl + ovilen + ovavol,
ovilen ~ ovavol,
ccl    ~ ovilen,
ccl    ~ ovavol)
semMod
latent(semMod) <- ~ ccl + ovilen + ovavol
semMod
#latent(semMod) <- ~ ccl + ovilen + ovavol
regression(semMod, ovilen ~ ovavol) <- 1
regression(semMod, ccl ~ ovilen) <- 1
regression(semMod, ccl ~ ovavol) <- 1
regression(semMod) <- list(
maturity ~ ccl + ovilen + ovavol,
ovilen ~ ovavol,
ccl    ~ ovilen,
ccl    ~ ovavol)
#latent(semMod) <- ~ ccl + ovilen + ovavol
regression(semMod, ovilen ~ ovavol) <- 1
regression(semMod, ccl ~ ovilen) <- 1
regression(semMod, ccl ~ ovavol) <- 1
semMod
regression(semMod) <- list(
maturity ~ ccl + ovilen + ovavol,
ovilen ~ ovavol,
ccl    ~ ovilen,
ccl    ~ ovavol)
semMod
# This is a package lava SEM
semMod <- lvm()
regression(semMod) <- list(
maturity ~ ccl + ovilen + ovavol,
ovilen ~ ovavol,
ccl    ~ ovilen,
ccl    ~ ovavol)
semMod
#latent(semMod) <- ~ ccl + ovilen + ovavol
#regression(semMod, ovilen ~ ovavol) <- 1
#regression(semMod, ccl ~ ovilen) <- 1
#regression(semMod, ccl ~ ovavol) <- 1
#covariance(semMod) <- ovilen ~ ovavol # error correlations between latents
#covariance(semMod) <- ccl ~ ovavol # error correlations between latents
#covariance(semMod) <- ccl ~ ovilen # error correlations between latents
distribution(semMod, maturity ~ ccl + ovilen + ovavol) <- binomial.lvm("logit")
distribution(semMod, ccl ~ ovilen + ovavol) <- gaussian.lvm("log")
semMod
lava.options(param = "hybrid")
estimate(semMod, BRED_Data, estimator="glm")
?
estimate
#latent(semMod) <- ~ ccl + ovilen + ovavol
#regression(semMod, ovilen ~ ovavol) <- 1
#regression(semMod, ccl ~ ovilen) <- 1
#regression(semMod, ccl ~ ovavol) <- 1
covariance(semMod) <- ovilen ~ ovavol # error correlations between latents
covariance(semMod) <- ccl ~ ovavol # error correlations between latents
covariance(semMod) <- ccl ~ ovilen # error correlations between latents
distribution(semMod, maturity ~ ccl + ovilen + ovavol) <- binomial.lvm("logit")
distribution(semMod, ccl ~ ovilen + ovavol) <- gaussian.lvm("log")
lava.options(param = "hybrid")
estimate(semMod, BRED_Data, estimator="glm")
semMod
# This is a package lava SEM
semMod <- lvm()
regression(semMod) <- list(
maturity ~ ccl + ovilen + ovavol,
ovilen ~ ovavol,
ccl    ~ ovilen,
ccl    ~ ovavol)
#latent(semMod) <- ~ ccl + ovilen + ovavol
regression(semMod, ovilen ~ ovavol) <- 1
regression(semMod, ccl ~ ovilen) <- 1
regression(semMod, ccl ~ ovavol) <- 1
#covariance(semMod) <- ovilen ~ ovavol # error correlations between latents
#covariance(semMod) <- ccl ~ ovavol # error correlations between latents
#covariance(semMod) <- ccl ~ ovilen # error correlations between latents
distribution(semMod, maturity ~ ccl + ovilen + ovavol) <- binomial.lvm("logit")
distribution(semMod, ccl ~ ovilen + ovavol) <- gaussian.lvm("log")
lava.options(param = "hybrid")
estimate(semMod, BRED_Data, estimator="glm")
semMod
# This is a package lava SEM
semMod <- lvm()
regression(semMod) <- list(
maturity ~ ccl + ovilen + ovavol,
ovilen ~ ovavol,
ccl    ~ ovilen,
ccl    ~ ovavol)
#latent(semMod) <- ~ ccl + ovilen + ovavol
regression(semMod, ovilen ~ ovavol) <- 1
regression(semMod, ccl ~ ovilen) <- 1
regression(semMod, ccl ~ ovavol) <- 1
covariance(semMod) <- ovilen ~ ovavol # error correlations between latents
covariance(semMod) <- ccl ~ ovavol # error correlations between latents
covariance(semMod) <- ccl ~ ovilen # error correlations between latents
distribution(semMod, maturity ~ ccl + ovilen + ovavol) <- binomial.lvm("logit")
distribution(semMod, ccl ~ ovilen + ovavol) <- gaussian.lvm("log")
lava.options(param = "hybrid")
estimate(semMod, BRED_Data, estimator="glm")
semMod
estimate(semMod, BRED_Data, estimator="glm")
semMod
gaussian.lvm
## Label new parameters
estimate(g,function(p) list("a1"=p[1]+p[2],"b1"=p[1]*p[2]))
##'
## Multiple group
m <- lvm(y~x)
m <- baptize(m)
d2 <- d1 <- sim(m,50)
e <- estimate(list(m,m),list(d1,d2))
estimate(e) ## Wrong
estimate(e,id=rep(seq(nrow(d1)),2))
estimate(lm(y~x,d1))
?lmv
?lvm
q()
## Load libraries
lib.fn <- function(){
library(readxl)
library(MASS)
library(stats)
library(tidyverse)
library(stringr)
library(bbmle)
library(AER)
library(sjPlot)
library(glmmTMB)
library(jtools)
library(multcomp)
library(ggplot2)
library(reshape2)
library(DHARMa)
library(broom.mixed)
library(AICcmodavg)
library(MuMIn)
library(effects)
library(emmeans)
library(ggthemes)
library(patchwork)
}; lib.fn()
## Import data from Excel files
fun <- new.env() # Create environments.
fun$importExcel <- function(fileName, sheetNames, nskip=0){ # Create a function to call in various tabs in Excel
df <- lapply(sheetNames, function(x){
read_excel(fileName, sheet=x, skip=nskip, col_types = "text")
})
names(df) = sheetNames
df
}
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")
## underdispersion still present in the BEACH model
## 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")
x <- emmeans(OCF.Model, ~ Beach, transform = "response")
str(x)
x[1]
x[2]
x[[2]]
x[[1]
]
x[[1]]
x
1.96*0.0503
3.54 + 1.96*0.0503
3.54 - 1.96*0.0503
q()