# R/misseNet.R In hisse: Hidden State Speciation and Extinction

#### Documented in MiSSENet

```######################################################################################################################################
######################################################################################################################################
### Assessses model optimization and triages those that fail
######################################################################################################################################
######################################################################################################################################

#Takes a vector of free parameters and find the distance between them and gives back an
#edge matrix for models 1 freepar away from each other
GetEdges <- function(turn.free.p, eps.free.p, nodes){
#Step 1 get model distance based on turnover
distances.turn <- dist(turn.free.p)
dist.mat.turn <- as.matrix(distances.turn)
dist.mat.turn[lower.tri(dist.mat.turn)] <- max(dist.mat.turn)

#Step 2: Get model distance based on eps
distances.eps <- dist(eps.free.p)
dist.mat.eps <- as.matrix(distances.eps)
dist.mat.eps[lower.tri(dist.mat.eps)] <- max(dist.mat.eps)

#Step 3: Get model diff of 1 in turn direction, zero for eps
one.par.diff.turn <- which(dist.mat.turn == 1 & dist.mat.eps == 0, arr.ind=TRUE)
#Step 4: Get model diff of 1 in eps direction, zero for turn
one.par.diff.eps <- which(dist.mat.eps == 1 & dist.mat.turn == 0, arr.ind=TRUE)
#Step 5: Combine Step 3 and 4
one.par.diff <- rbind(one.par.diff.turn, one.par.diff.eps)
#Step 5: Loop through table to assign proper node names for use in igraph
tmp.edges.mat <- c()
for(index in 1:dim(one.par.diff)[1]){
tmp.edges.mat <- rbind(tmp.edges.mat, c(nodes[one.par.diff[index,1]], nodes[one.par.diff[index,2]]))
}
final.edge.mat <- data.frame(x=tmp.edges.mat[,1], y=tmp.edges.mat[,2], weight=rep(1,dim(tmp.edges.mat)[1]))
#Fin.
return(final.edge.mat)
}

for(node.index in 1:length(nodes)){
desc.nodes <- igraph::neighbors(graph, igraph::V(graph)[node.index])
for(neighbor.index in 1:length(desc.nodes)){
if(length(desc.nodes) > 0){
if(round(loglik.vec[node.index],3) - round(loglik.vec[desc.nodes[neighbor.index]],3) > 0.1){
}
}
}
}
return(0)
}else{
}
}

RerunBadOptim <- function(bad.fits, misse.list, sann, sann.its, sann.temp, bounded.search, starting.vals, turnover.upper, eps.upper, trans.upper, restart.obj, retries){

if(sann == FALSE){
sann = TRUE
bounded.search = TRUE
}
print(paste("current loglik:", current.loglik))
sann.its <- sann.its * 2
sann.temp <- sann.temp
sann.seed <- c(-487281, -391945, -149149, -919193, -682017)[retries[[bad.fits]]+1]
print(paste("new loglik:", redo.fit\$loglik))
if(redo.fit\$loglik > current.loglik){
return(redo.fit)
}
}

MiSSENet <- function(misse.list, n.tries=2, remove.bad=TRUE, dont.rerun=FALSE, save.file=NULL, n.cores=1, sann=TRUE, sann.its=5000, sann.temp=5230, bounded.search=TRUE, starting.vals=NULL, turnover.upper=10000, eps.upper=3, trans.upper=100, restart.obj=NULL){

#Step 1: Make igraph of examined model space
tmp <- c()
for(model.index in 1:length(misse.list)){
tmp <- rbind(tmp, c(length(unique(misse.list[[model.index]]\$turnover)), length(unique(misse.list[[model.index]]\$eps)), misse.list[[model.index]]\$loglik, misse.list[[model.index]]\$AICc))
}
model.space <- data.frame(turnover=tmp[,1], eps=tmp[,2], loglik=tmp[,3], aic=tmp[,4])
nodes <- paste("T", model.space\$turnover, "_", "E", model.space\$eps, sep="")
edges <- GetEdges(model.space\$turnover, model.space\$eps, nodes)
graph.df <- igraph::graph.data.frame(edges, vertices=nodes)

#Step 2: Locate the bad fits by traversing graph in step 1

if(dont.rerun == TRUE){
}
misse.list <- Filter(Negate(anyNA), misse.list)
return(misse.list)
}else{
}
}else{
model.retries <- as.list(as.list(numeric(length(misse.list))))
done.enough <- FALSE
while(done.enough == FALSE){
cat("Current set of model fits identified as poorly optimized:", bad.fits, "\n")
#Step 3: Rerun the bad fits by increasing simulated annealing temperature? Still trying to figure out
rerun.fits <- parallel::mclapply(bad.fits, RerunBadOptim, misse.list=misse.list, sann=sann, sann.its=sann.its, sann.temp=sann.temp, bounded.search=bounded.search, starting.vals=starting.vals, turnover.upper=turnover.upper, eps.upper=eps.upper, trans.upper=trans.upper, restart.obj=restart.obj, retries=model.retries, mc.cores=ifelse(is.null(n.cores),1,n.cores))
### Keeping this here for debugging purposes ###
################################################
for(rerun.index in 1:length(rerun.fits)){
}

#Repeat step 1:
tmp <- c()
for(model.index in 1:length(misse.list)){
tmp <- rbind(tmp, c(length(unique(misse.list[[model.index]]\$turnover)), length(unique(misse.list[[model.index]]\$eps)), misse.list[[model.index]]\$loglik, misse.list[[model.index]]\$AICc))
}
model.space <- data.frame(turnover=tmp[,1], eps=tmp[,2], loglik=tmp[,3], aic=tmp[,4])
nodes <- paste("T", model.space\$turnover, "_", "E", model.space\$eps, sep="")
edges <- GetEdges(model.space\$turnover, model.space\$eps, nodes)
graph.df <- igraph::graph.data.frame(edges, vertices=nodes)
done.enough = TRUE
}else{
cat("One or more models still identified as poorly optimized.", "\n")
}
}
cat("Maximum number of reassessments reached. These models should be discarded:", bad.fits, "\n")
if(!is.null(save.file)) {
save(misse.list, file=save.file)
}
}
misse.list <- Filter(Negate(anyNA), misse.list)
}
}
}
return(misse.list)
}else{
cat("All models appeared to have optimized well.", "\n")
return(misse.list)
}
}
}
```

## Try the hisse package in your browser

Any scripts or data that you put into this service are public.

hisse documentation built on Nov. 16, 2021, 9:09 a.m.