Skip to content

Commit

Permalink
updates for 2024
Browse files Browse the repository at this point in the history
  • Loading branch information
kholsman committed Jan 8, 2024
1 parent 2132aa0 commit 6bf69a4
Show file tree
Hide file tree
Showing 10 changed files with 464 additions and 403 deletions.
254 changes: 126 additions & 128 deletions R/sub_fun/AIC_selection.R
View file
Original file line number Diff line number Diff line change
@@ -1,129 +1,127 @@
#-------------------------------------
# AICselection
#-------------------------------------
#' AICselection runs AIC information criterion analysis on resulting set of recruitment models
#' For more information contact Kirstin Holsman (kirstin.holsman@noaa.gov)
#' @weblink
#' @param data data.frame with models = rows and the following data for each model (columns):
#' * data$LL vector of logliklihoods of each model
#' * data$npar number of parameteris in a given model
#' * data$n number of data observations (nobs)
#' * data$mnames1 names for each model
#' * data$R2 pearson correlation coefficent value for the model (hat:obs); default = NULL
#' * data$type2 number 1= AIC;2 = AICc ; default = NULL
#' * data$covnm default = NULL
#' @param rsType default = NULL
#' @param LnDet default = NULL
#' @param CE Climate enhanced model? default = NULL
#' list of input data from makeDat() function: "parameters" "rs_dat" "maplist" "estparams" "phases"
#' @param returnAll return all phases? TRUE/FALSE
#' @param quiet print out phases
#' @return tmp1 a dataframe of summary statistics for AICc of submodels
#' @examples
#' @export
AICselection <- function (LL,
npar,
n,
mnames1 = legend.nm,
R2 = R2[, s],
type2 = 2,
covnm = NULL,
rsType = NULL,
hypoth = NULL,
LnDet = NULL,
CE = NULL)
{
# LL<--1*LL
#
aicfun<-function (npar, LL, n, type = 2)
{
if (type == 1)
return((2 * npar - 2 * (LL)))
if (type == 2)
return((2 * npar - 2 * (LL)) + (2 * npar * (npar + 1))/(n - npar - 1))
}


aicfun_marg<-function (npar, LL, n, type = 2,LnDet)
{
#LnDet = determinant(HESS, logarithm=TRUE)$modulus[[1]]
Ln_Integral <- log(2*pi) + (-1/2)*(LnDet) + LL #this is the MARGINAL likelihood

if (type == 1)
return((2 * npar - 2 * (Ln_Integral)))
if (type == 2)
return((2 * npar - 2 * (Ln_Integral)) + (2 * npar * (npar + 1))/(n - npar - 1))
}

nn <- length(LL)
tmp1 <- data.frame(LL = LL, npar = npar, lab = 1:nn, R2 = R2)
tmp1$names <- mnames1
tmp1$aicc <- tmp1$aicc_marg <- rep(0, nn)
for (i in 1:nn) {
tmp1$aicc[i] <- aicfun(npar[i], -LL[i], n[i], type = type2)
if(!is.null(LnDet))
tmp1$aicc_marg[i] <- aicfun_marg(npar=npar[i], LL=-LL[i],n= n[i], type = type2,LnDet=LnDet[i])
# tmp1$aicc_marg[i] <- GET_HESS_AIC(HESS[[i]], npar = npar[i],NLL = -1 * LL[i])[[1]]
}
tmp1$deltaAIC <- (tmp1$aicc - min(tmp1$aicc))
tmp1$AICweight <- exp(-0.5 * tmp1$deltaAIC)
tmp1$rank <- rank(tmp1$aicc)
if(!is.null(covnm))
tmp1$covnm = covnm
if(!is.null(rsType))
tmp1$rsType = rsType
if(!is.null(hypoth))
tmp1$hypoth = hypoth
if(!is.null(CE))
tmp1$CE = CE

tmp1 <- tmp1[order(tmp1$aicc), ]
tmp1$AICw_std <- tmp1$AICweight/sum(tmp1$AICweight)
tmp1$cumlAIC <- cumsum(tmp1$AICw_std)

cutoff <- which(tmp1$cumlAIC > 0.95)[1]
if (is.na(cutoff))
cutoff <- nn
cutoff2 <- which(tmp1$deltaAIC > 2)[1]
if (is.na(cutoff2))
cutoff2 <- nn
cutoff4 <- which(tmp1$deltaAIC > 4)[1]
if (is.na(cutoff4))
cutoff4 <- nn

t1 <- rep("", nn)
t1[1:cutoff] <- "o"
t2 <- rep("", nn)
t2[1:cutoff2] <- "*"
t4 <- rep("", nn)
t4[1:cutoff4] <- "*"


tmp1$topSet <- paste(t2, t4, t1, sep = "")

tmp1$aicc_marg[tmp1$aicc_marg==-Inf] <- NA
tmp1$deltaAIC_marg <- (tmp1$aicc_marg - min(tmp1$aicc_marg,na.rm=T))
tmp1$AICweight_marg <- exp(-0.5 * tmp1$deltaAIC_marg)
tmp1$rank_marg <- rank(tmp1$aicc_marg)
tmp1$AICw_std_marg <- tmp1$AICweight_marg/sum(tmp1$AICweight_marg,na.rm=T)
tmp1$AICw_std_marg[is.na(tmp1$AICw_std_marg)]<-0
tmp1$cumlAIC_marg <- cumsum(tmp1$AICw_std_marg)
cutoff_marg <- which(tmp1$cumlAIC_marg > 0.95)[1]

if (is.na(cutoff_marg))
cutoff_marg <- nn
cutoff2_marg <- which(tmp1$deltaAIC_marg > 2)[1]
if (is.na(cutoff2_marg))
cutoff2_marg <- nn
cutoff4_marg <- which(tmp1$deltaAIC_marg > 4)[1]
if (is.na(cutoff4_marg))
cutoff4_marg <- nn
t1_marg <- rep("", nn)
t1_marg[1:cutoff_marg] <- "o"
t2_marg <- rep("", nn)
t2_marg[1:cutoff2_marg] <- "*"
t4_marg <- rep("", nn)
t4_marg[1:cutoff4_marg] <- "*"
tmp1$topSet_marg <- paste(t2_marg, t4_marg, t1_marg, sep = "")
return(tmp1)
#-------------------------------------
# AICselection
#-------------------------------------
#' AICselection runs AIC information criterion analysis on resulting set of recruitment models
#' For more information contact Kirstin Holsman (kirstin.holsman@noaa.gov)
#' @weblink
#' @param data data.frame with models = rows and the following data for each model (columns):
#' * data$LL vector of logliklihoods of each model
#' * data$npar number of parameteris in a given model
#' * data$n number of data observations (nobs)
#' * data$mnames1 names for each model
#' * data$type2 number 1= AIC;2 = AICc ; default = NULL
#' * data$covnm default = NULL
#' @param rsType default = NULL
#' @param LnDet default = NULL
#' @param CE Climate enhanced model? default = NULL
#' list of input data from makeDat() function: "parameters" "rs_dat" "maplist" "estparams" "phases"
#' @param returnAll return all phases? TRUE/FALSE
#' @param quiet print out phases
#' @return tmp1 a dataframe of summary statistics for AICc of submodels
#' @examples
#' @export
AICselection <- function (LL,
npar,
n,
mnames1 = legend.nm,
type2 = 2,
covnm = NULL,
rsType = NULL,
hypoth = NULL,
LnDet = NULL,
CE = NULL)
{
# LL<--1*LL
#
aicfun<-function (npar, LL, n, type = 2)
{
if (type == 1)
return((2 * npar - 2 * (LL)))
if (type == 2)
return((2 * npar - 2 * (LL)) + (2 * npar * (npar + 1))/(n - npar - 1))
}


aicfun_marg<-function (npar, LL, n, type = 2,LnDet)
{
#LnDet = determinant(HESS, logarithm=TRUE)$modulus[[1]]
Ln_Integral <- log(2*pi) + (-1/2)*(LnDet) + LL #this is the MARGINAL likelihood

if (type == 1)
return((2 * npar - 2 * (Ln_Integral)))
if (type == 2)
return((2 * npar - 2 * (Ln_Integral)) + (2 * npar * (npar + 1))/(n - npar - 1))
}

nn <- length(LL)
tmp1 <- data.frame(LL = LL, npar = npar, lab = 1:nn)
tmp1$name <- mnames1
tmp1$aicc <- tmp1$aicc_marg <- rep(0, nn)
for (i in 1:nn) {
tmp1$aicc[i] <- aicfun(npar[i], -LL[i], n[i], type = type2)
if(!is.null(LnDet))
tmp1$aicc_marg[i] <- aicfun_marg(npar=npar[i], LL=-LL[i],n= n[i], type = type2,LnDet=LnDet[i])
# tmp1$aicc_marg[i] <- GET_HESS_AIC(HESS[[i]], npar = npar[i],NLL = -1 * LL[i])[[1]]
}
tmp1$deltaAIC <- (tmp1$aicc - min(tmp1$aicc))
tmp1$AICweight <- exp(-0.5 * tmp1$deltaAIC)
tmp1$rank <- rank(tmp1$aicc)
if(!is.null(covnm))
tmp1$covnm = covnm
if(!is.null(rsType))
tmp1$rsType = rsType
if(!is.null(hypoth))
tmp1$hypoth = hypoth
if(!is.null(CE))
tmp1$CE = CE

tmp1 <- tmp1[order(tmp1$aicc), ]
tmp1$AICw_std <- tmp1$AICweight/sum(tmp1$AICweight, na.rm = T)
tmp1$cumlAIC <- cumsum(tmp1$AICw_std)

cutoff <- which(tmp1$cumlAIC > 0.95)[1]
if (is.na(cutoff))
cutoff <- nn
cutoff2 <- which(tmp1$deltaAIC > 2)[1]
if (is.na(cutoff2))
cutoff2 <- nn
cutoff4 <- which(tmp1$deltaAIC > 4)[1]
if (is.na(cutoff4))
cutoff4 <- nn

t1 <- rep("", nn)
t1[1:cutoff] <- "o"
t2 <- rep("", nn)
t2[1:cutoff2] <- "*"
t4 <- rep("", nn)
t4[1:cutoff4] <- "*"


tmp1$topSet <- paste(t2, t4, t1, sep = "")

tmp1$aicc_marg[tmp1$aicc_marg==-Inf] <- NA
tmp1$deltaAIC_marg <- (tmp1$aicc_marg - min(tmp1$aicc_marg,na.rm=T))
tmp1$AICweight_marg <- exp(-0.5 * tmp1$deltaAIC_marg)
tmp1$rank_marg <- rank(tmp1$aicc_marg)
tmp1$AICw_std_marg <- tmp1$AICweight_marg/sum(tmp1$AICweight_marg,na.rm=T)
tmp1$AICw_std_marg[is.na(tmp1$AICw_std_marg)]<-0
tmp1$cumlAIC_marg <- cumsum(tmp1$AICw_std_marg)
cutoff_marg <- which(tmp1$cumlAIC_marg > 0.95)[1]

if (is.na(cutoff_marg))
cutoff_marg <- nn
cutoff2_marg <- which(tmp1$deltaAIC_marg > 2)[1]
if (is.na(cutoff2_marg))
cutoff2_marg <- nn
cutoff4_marg <- which(tmp1$deltaAIC_marg > 4)[1]
if (is.na(cutoff4_marg))
cutoff4_marg <- nn
t1_marg <- rep("", nn)
t1_marg[1:cutoff_marg] <- "o"
t2_marg <- rep("", nn)
t2_marg[1:cutoff2_marg] <- "*"
t4_marg <- rep("", nn)
t4_marg[1:cutoff4_marg] <- "*"
tmp1$topSet_marg <- paste(t2_marg, t4_marg, t1_marg, sep = "")
return(tmp1)
}
10 changes: 5 additions & 5 deletions R/sub_fun/getFit.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,15 +14,15 @@
#' plot_rs(r_fit)
#' @export
getFit <- function(mm, nm = "recType = 2"){
out <- mm$pred %>%
out <- data.frame(mm$pred %>%
dplyr::select(def,pred)%>%
dplyr::group_by(def) %>%
dplyr::filter(def%in%c("S_hat","S_obs","R_hat","R_obs"))%>%
dplyr::mutate(id = 1:n(), model = nm) %>%
dplyr::mutate(id = 1:n()) %>%
tidyr::spread(def, pred)%>%
dplyr::mutate(R_obsIN = mm$input$R_obs,
S_obsIN = mm$input$S_obs,
year = mm$input$years)%>%
dplyr::ungroup()
return(data.frame(out))
year = mm$input$years,model = nm)%>%
dplyr::ungroup())
return(out)
}
47 changes: 47 additions & 0 deletions R/sub_fun/get_interaction_mat.R
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,47 @@
#'
#'
#'
#'get_interaction_mat.R
#'
#'@param covIN covars (matrix with cols for each covar) e.g., pre_spawning_covs[,-1]
#'@param maxIN max number of interactions
#'@param ADDTEMP2 add a temperature ^2 term to those covars with "temp" in the name
#'@param cor_cutoff cut off value for correlation matrix (e.g., 0.5), those with cor > than this value will be excluded
#' get full set of 4 interactions (max_interactions)
#' returns a list with data.frames of all possible combinations

get_interaction_mat<-function(covIN, maxIN,ADDTEMP2=TRUE,cor_cutoff=0.5){
mod_mat <- list()
for(i in 1:maxIN){
tt <- combn(covIN,i,simplify=F)
for(k in 1:length(tt)){
if(i>1){
ex <- cor(tt[[k]])
diag(ex)<-0
}else{
ex = 0
}
names(tt)[k] <- paste0(names(tt[[k]]),collapse="_PLUS_")
# if the covars are not related:
if(!any( abs(ex) > cor_cutoff ) ){
mod_mat <- c(mod_mat,tt[k])
}
rm(ex)
if(ADDTEMP2){
aa <- grep("temp",colnames(tt[[k]]))
if(length(aa)>0){
# add temperature squared
tmp <- data.frame(tt[[k]][,aa]^2)
colnames(tmp) <- paste0(colnames(tt[[k]])[aa],"x2")
ll <- list(data.frame( tt[[k]],tmp))
names(ll) <- paste0(names(ll[[1]]),collapse="_PLUS_")
mod_mat <- c(mod_mat,ll)
rm(tmp)
}
rm(aa)
}
}
rm(tt)
}
return(mod_mat)
}
Loading

0 comments on commit 6bf69a4

Please sign in to comment.