# R/select.R In ctmm: Continuous-Time Movement Modeling

#### Documented in ctmm.select

```# small sample size adjustment for ctmm.select to be more agressive
alpha.ctmm <- function(CTMM,alpha)
{
z <- stats::qnorm(alpha)
z <- sqrt(z^2 + (CTMM\$AICc-CTMM\$AIC))
alpha <- 1-stats::pnorm(z)
if(!length(alpha)) { alpha <- 0 }
return(alpha)
}

#########
get.MSPE <- function(CTMM,MSPE="position")
{
if(!is.na(MSPE)) { MSPE <- CTMM\$MSPE[MSPE] }
else { MSPE <- Inf }
return(MSPE)
}

########
get.IC <- function(CTMM,IC="AICc")
{
if(!is.na(IC)) { IC <- CTMM[[IC]] }
else { IC <- Inf }
return(IC)
}

##################
# function to simplify complexity of models
simplify.ctmm <- function(M,par)
{
if("minor" %in% par)
{
if("COV" %in% names(M))
{
M\$COV <- axes2var(M)
NAMES <- colnames(M\$COV)
VAR <- which(NAMES=="variance")
colnames(M\$COV)[VAR] <- rownames(M\$COV)[VAR] <- "major"
}

M\$isotropic <- TRUE
M\$sigma <- covm(M\$sigma,isotropic=TRUE,axes=M\$axes)
par <- c(par,'angle')
}

if("major" %in% par)
{
M\$isotropic <- TRUE
M\$sigma <- covm(0,isotropic=TRUE,axes=M\$axes)
M\$tau <- NULL
par <- c(par,c('circle','tau position','tau velocity','tau','omega'))
}

if("circle" %in% par)
{ M\$circle <- FALSE }

if("range" %in% par && length(M\$tau))
{
# convert to diffusion matrix
if(M\$tau[1]>0) { M\$sigma <- scale.covm(M\$sigma,1/M\$tau[1]) }
M\$tau[1] <- Inf
M\$range <- FALSE
}

# # reduce time-link parameters
# if("timelink" %in% par) { M <- timelink.simplify(M) }

# autocorrelation timescales can't be distinguished
if("diff.tau" %in% par)
{
M\$tau <- c(1,1)*mean(M\$tau)
M\$omega <- FALSE
}

if('tau velocity' %in% par) { M\$tau <- M\$tau[1] }
if(any(c('tau position','tau') %in% par)) { M\$tau <- NULL }
if('omega' %in% par) { M\$omega <- FALSE }
if('circle' %in% par) { M\$circle <- FALSE }
TEST <- grepl('error',par)
if(any(TEST))
{
PAR <- par[TEST]
PAR <- substr(PAR,nchar("error ?"),nchar(PAR))
M\$error[PAR] <- FALSE # why would you ever do this?
}

# re-identify features
UERE.FIT <- paste("error",names(M\$error)) %in% M\$features
names(UERE.FIT) <- names(M\$error)

M\$features <- id.parameters(M,profile=FALSE,linear=FALSE,UERE.FIT=UERE.FIT)\$NAMES

# don't think we need this
if("MLE" %in% names(M)) { M\$MLE <- simplify.ctmm(M\$MLE,par) }

return(M)
}

#############
# function to make autocorrelation models more complex
complexify.ctmm <- function(M,par,TARGET)
{
# consider non-zero eccentricity
if(("minor" %in% par) && M\$isotropic)
{
M\$isotropic <- FALSE
# copy over target angle, but leave eccentricity zero to start (featureless)
sigma <- attr(M\$sigma,"par")
sigma["angle"] <- attr(TARGET\$sigma,"par")['angle']
sigma <- covm(sigma,isotropic=FALSE,axes=TARGET\$axes)
M\$sigma <- sigma
}

# consider circulation
if(("circle" %in% par) && !M\$circle)
{  M\$circle <- 2 * .Machine\$double.eps * sign(TARGET\$circle) }

# consider finite range
if(("range" %in% par) && !M\$range)
{
# convert from diffusion matrix
M\$tau[1] <- TARGET\$tau[1]
M\$sigma <- scale.covm(M\$sigma,TARGET\$tau[1])
M\$range <- TRUE
}

# # increase time-link parameters
# if("timelink" %in% par) { M <- timelink.complexify(M) }

if("tau velocity" %in% par) { M\$tau[2] <- TARGET\$tau[2] }
if("omega" %in% par) { M\$omega <- TARGET\$omega }

# re-identify features
UERE.FIT <- M\$error>0 & (paste('error',names(M\$error)) %in% M\$features)
# if(M\$error && ('error' %in% M\$features)) { UERE <- TRUE } else { UERE <- 3 } # calibrated or not?
M\$features <- id.parameters(M,profile=FALSE,linear=FALSE,UERE.FIT=UERE.FIT)\$NAMES

# don't think we need this
if("MLE" %in% names(M)) { M\$MLE <- complexify.ctmm(M\$MLE,par,TARGET) }

return(M)
}

# initial guess in case of pREML/HREML/pHREML (better for optimization)
get.mle <- function(FIT)
{
if("MLE" %in% names(FIT))
{
axes <- FIT\$axes
MLE <- FIT\$MLE
FIT\$MLE <- NULL
# have parameters been altered after the fact?
if(MLE\$checksum==digest::digest(FIT,algo='md5')) { FIT <- MLE }
FIT\$axes <- axes # goes missing in MLE?
}
return(FIT)
}

###############
# keep removing uncertain parameters until AIC stops improving
ctmm.select <- function(data,CTMM,verbose=FALSE,level=1,IC="AICc",MSPE="position",trace=FALSE,cores=1,...)
{
M <- ctmm.iterate(data,CTMM,verbose=TRUE,level=level,IC=IC,MSPE=MSPE,trace=trace,cores=cores,recurse=FALSE,TRYS=NULL,...)

for(i in 1:length(M))
{
if(grepl('anisotropic',names(M)[i]))
{
ISO <- gsub(" anisotropic","",names(M)[i],fixed=TRUE)
if(ISO %in% names(M)) { M[[i]]\$ISO <- M[[ISO]] }
}
}

if(verbose) { return(M) }
else { return(M[[1]]) }
}

## internal function being called recursively by ctmm.select
# recurse: FALSE is outer-most call, TRUE are inner calls
# TRYS: models that we have tried
ctmm.iterate <- function(data,CTMM,verbose=FALSE,level=1,IC="AICc",MSPE="position",trace=FALSE,cores=1,recurse=FALSE,TRYS=NULL,...)
{
LIST <- list(...)
if(!is.null(LIST\$control\$message)) { message <- LIST\$control\$message }

CV <- c("LOOCV","HSCV")
IC <- match.arg(IC,c("AICc","AIC","BIC",CV,NA))
MSPE <- match.arg(MSPE,c("position","velocity",NA))

# CV: cross validation
# requires expensive post-fit calculation of IC
# can select between IID/OU/OUF and BM/IOU
CV <- IC %in% CV

alpha <- 1-level
trace2 <- if(trace) { trace-1 } else { 0 }

# accept list as completed candidates after first
# for internal recursion, not for end users
if(class(CTMM)[1]=="ctmm")
{
# TRYS <- NULL
MODELS <- list()
}
else if(class(CTMM)[1]=="list")
{
TRYS <- c(TRYS,attr(CTMM,"attempted")) # models that we have tried
MODELS <- CTMM[-1]
CTMM <- CTMM[[1]]
}
CAND <- length(MODELS) # number of extra candidate models included for features
if(CAND) { names(MODELS) <- sapply(MODELS,name.ctmm) }

TYPE <- DOP.match(CTMM\$axes)
UERE.DOF <- attr(data,"UERE")\$DOF[,TYPE]
names(UERE.DOF) <- rownames(attr(data,"UERE")\$DOF)
UERE.FIT <- CTMM\$error & !is.na(UERE.DOF) & UERE.DOF<Inf # will we be fitting any error parameters?

# best non-zero variance estimates -- avoid zero collapse in intermediate model propagating to best model
axes <- CTMM\$axes
AXES <- length(axes)
ERROR <- CTMM\$error # best non-zero error estimate
VAR <- CTMM\$sigma@par[1:AXES]
fix.vars <- function(M)
{
if(length(M))
{
for(i in 1:length(M))
{
TEST <- M[[i]]\$error <= .Machine\$double.eps
if(any(TEST)) { M[[i]]\$error[TEST] <- ERROR[TEST] }
M.VAR <- M[[i]]\$sigma@par
for(j in 1:AXES) { if(M.VAR[j]<=.Machine\$double.eps) { M.VAR[j] <- VAR[j] } }
M[[i]]\$sigma <- covm(M.VAR,axes=axes,isotropic=M[[i]]\$isotropic)
}
}
return(M)
}
update.fix <- function(CTMM)
{
TEST <- CTMM\$error > .Machine\$double.eps
if(any(TEST)) { ERROR[TEST] <<- CTMM\$error[TEST] }
for(i in 1:AXES) { if(CTMM\$sigma@par[i]>.Machine\$double.eps) { VAR[i] <<- CTMM\$sigma@par[i] } }
}

drift <- get(CTMM\$mean)
if(CTMM\$mean=="periodic")
{
Nyquist <- CTMM\$period/stats::median(diff(data\$t))/2
message("Nyquist frequency estimated at harmonic ",paste(Nyquist,collapse=" ")," of the period.")
}

# consider a bunch of new models and update best model without duplication
iterate <- function(DROP,REFINE=list(),phase=1)
{
# make sure everything tries errors
DROP <- fix.vars(DROP)
REFINE <- fix.vars(REFINE)

# name the proposed models
names(DROP) <- sapply(DROP,name.ctmm)
names(REFINE) <- sapply(REFINE,name.ctmm)

# remove models already attempted
DROP <- DROP[!(names(DROP) %in% TRYS)]
REFINE <- REFINE[!(names(REFINE) %in% TRYS)]

# copy over best initial parameter guess for refined drops
if(!drift@is.stationary(CTMM) && length(DROP))
{
for(i in 1:length(DROP))
{
for(j in 1:length(MODELS))
{
# take last/best model of same covariance structure
F1 <- DROP[[i]]\$features
F2 <- MODELS[[j]]\$features
if(length(F1)==length(F2) && all(F1==F2))
{
# copy over covariance parameters
DROP[[i]] <- copy.parameters(DROP[[i]],get.mle(MODELS[[j]]))
break # out of MODELS loop
}
} # end MODELS loop
} # end DROP loop
} # end refined drop adjustments

GUESS <- c(DROP,REFINE)
# add to TRYS before potential feature collapse
TRYS.OLD <- TRYS # don't count new attempts for next select
if(length(GUESS))
{
TRYS <<- c(TRYS,names(GUESS))
TRYS <<- unique(TRYS)
}

if(length(GUESS)>1) # keep selecting (recursive)
{
# fit every model
if(trace && FALSE) { message("* Fitting models ",paste(names(GUESS),collapse=", "),".") }
# GUESS <- plapply(GUESS,function(g){M<-c(list(g),MODELS); attr(M,"attempted")<-TRYS.OLD; ctmm.select(data,M,verbose=TRUE,level=0,IC=IC,MSPE=MSPE,trace=trace,...)},cores=cores)
GUESS <- plapply(GUESS,function(g){M<-c(list(g),MODELS); ctmm.iterate(data,M,verbose=TRUE,level=0,IC=IC,MSPE=MSPE,trace=trace,recurse=TRUE,TRYS=TRYS.OLD,...)},cores=cores)
GUESS[sapply(GUESS,is.null)] <- NULL # delete collapsed repeats
if(length(GUESS)) # concatenate list of lists
{
ATTEMPT <- do.call(c, lapply(GUESS,function(g){attr(g,"attempted")}) ) # do.call-lapply in case of NULL
ATTEMPT[sapply(ATTEMPT,is.null)] <- NULL # delete NULL results... ? does this still happen ?
TRYS <- c(TRYS,ATTEMPT)
TRYS <- unique(TRYS)
GUESS <- do.call(c,GUESS)
}
}
else if(length(GUESS)==1) # fit last model
{
if(trace) { message("* Fitting model ",names(GUESS)) }
GUESS[[1]] <- ctmm.fit(data,GUESS[[1]],trace=trace2,...)

# cross-validate
if(CV)
{
if(trace) { message("** Cross validating model ",names(GUESS)[1]) }
GUESS[[1]][[IC]] <- do.call(IC,list(data=data,CTMM=GUESS[[1]],cores=cores,...))
}
}
# corrected name after potential feature collapse
names(GUESS) <- sapply(GUESS,name.ctmm)
# add to TRYS after potential feature collapse
if(length(GUESS))
{
TRYS <<- c(TRYS,names(GUESS))
TRYS <<- unique(TRYS)
}
# TRYS can be longer than MODELS

MODELS <<- c(GUESS,MODELS)

# only sort newer models, not input candidate models, which we will not return
N <- length(MODELS) - CAND
if(N>1) { MODELS[1:N] <<- sort.ctmm(MODELS[1:N],IC=IC,MSPE=MSPE) }

# what is the new best model?
OLD <<- CTMM
CTMM <<- MODELS[[1]]
update.fix(CTMM) # update best non-zero error estimate
# CTMM <<- min.ctmm(c(GUESS,list(CTMM)),IC=IC,MSPE=MSPE)
} # end iterate()

########################
# PHASE 0: pair down to essential features for 'compatibility'
# all of the features we need to fit numerically
FEATURES <- id.parameters(CTMM,UERE.FIT=UERE.FIT)\$NAMES
# consider only features unnecessary for Stein "compatibility"
FEATURES <- FEATURES[!(FEATURES=="major")]
FEATURES <- FEATURES[!grepl("error",FEATURES)]
FEATURES <- FEATURES[!grepl("tau",FEATURES)]
FEATURES <- FEATURES[!(FEATURES=="omega")]
# if((CV || is.na(IC)) && CTMM\$range && length(CTMM\$tau)) { FEATURES <- c(FEATURES,"range") }

# # clean up non-parametrics
# {
#   FEATURES <- FEATURES[!grepl("timelink",FEATURES)]
#   FEATURES <- c(FEATURES,"timelink")
# }

TARGET <- CTMM
# start with the most basic "compatible" model, if not included in candidates
GUESS <- simplify.ctmm(CTMM,FEATURES)
NAME <- name.ctmm(GUESS)
if(NAME %in% TRYS) # we have tried the paired down model
{ CTMM <- MODELS[[1]] }
else # fit paired down model
{
if(trace) { message("* Fitting model ",NAME) }
TRYS <- c(NAME,TRYS)
CTMM <- ctmm.fit(data,GUESS,trace=trace2,...)
if(CV)
{
if(trace) { message("** Cross validating model ",name.ctmm(GUESS)) }
CTMM[[IC]] <- do.call(IC,list(data=data,CTMM=CTMM,cores=cores,...))
}
TRYS <- c(name.ctmm(CTMM),TRYS)
TRYS <- unique(TRYS)
MODELS <- c(list(CTMM),MODELS) # update initial guess to fit
update.fix(CTMM)
}
names(MODELS) <- sapply(MODELS,name.ctmm)

##############################
# PHASE 1: work our way up to complicated autocorrelation models
OLD <- ctmm()
while(!identical(CTMM,OLD))
{
MLE <- get.mle(CTMM)
GUESS <- lapply(FEATURES,function(feat){complexify.ctmm(MLE,feat,TARGET)})

# consider a bunch of new models and update best model without duplication
iterate(GUESS)
}

#############################
# PHASE 2: work our way down to simpler autocorrelation models & work our way up to more complex trend models
# we only do phase 2 if(level), so we can call ctmm.select recursively to build up features again
OLD <- ctmm()
# CTMM <- min.ctmm(MODELS)
while(level && !identical(CTMM,OLD))
{
GUESS <- list()
MLE <- get.mle(CTMM)
beta <- alpha.ctmm(CTMM,alpha)

VARS <- dimnames(CTMM\$COV)[[1]] # could be NULL

# consider if smallest timescale is zero
CI <- confint.ctmm(CTMM,alpha=beta)
if(length(CTMM\$tau)==2 && !is.na(IC)) # OUX -> OUx
{
if("tau velocity" %in% VARS) # OUF -> OU / IOU -> BM
{
Q <- CI["tau velocity",1]
if(is.na(Q) || (Q<=0))
{ GUESS <- c(GUESS,list(simplify.ctmm(MLE,'tau velocity'))) }
}
else if("omega" %nin% VARS) # OUf -> OU, IID
{
Q <- CI["tau",1]
if(is.na(Q) || (Q<=0))
{
GUESS <- c(GUESS,list(simplify.ctmm(MLE,'tau'))) # OUf -> IID
GUESS <- c(GUESS,list(simplify.ctmm(MLE,'tau velocity'))) # OUf -> OU
}
}
else # OUO -> OUf
{
Q <- 1/CI["tau period",3]
if(is.na(Q) || (Q<=0))
{ GUESS <- c(GUESS,list(simplify.ctmm(MLE,'omega'))) }
}
} # end # OUX -> OUx
else if(CTMM\$range && length(CTMM\$tau)==1 && !is.na(IC)) # OU -> IID
{
Q <- CI["tau position",1]
if(is.na(Q) || (Q<=0))
{ GUESS <- c(GUESS,list(simplify.ctmm(MLE,'tau position'))) }
} # end # OU -> IID

# can autocorrelation timescales be distinguished?
if(all(c('tau position','tau velocity') %in% VARS) || "omega" %in% VARS) # isn't omega redundant?
{
TEMP <- get.taus(CTMM,zeroes=TRUE)
nu <- TEMP\$f.nu[2] # frequency/difference
J <- TEMP\$J.nu.tau[2,] # Jacobian of nu WRT canonical parameters
Q <- TEMP\$tau.names
Q <- c(J %*% CTMM\$COV[Q,Q] %*% J) # variance of nu
Q <- ci.tau(nu,Q,alpha=beta)[1]

if(is.na(Q) || Q<=0 || level==1 || is.na(IC))
{ GUESS <- c(GUESS,list(simplify.ctmm(MLE,"diff.tau"))) }
}
else if("tau" %in% VARS) # try other side if boundary if chosen model is critically damped
{
# try overdamped
TEMP <- TARGET
if(length(TEMP\$tau)<2 || !TEMP\$tau[2] || TEMP\$tau[2] >= MLE\$tau[1]) { TEMP\$tau <- MLE\$tau * exp(c(1,-1)*sqrt(.Machine\$double.eps)) }
GUESS <- c(GUESS,list(complexify.ctmm(MLE,"tau velocity",TEMP)))

# try underdamped
TEMP <- TARGET
if(!TEMP\$omega) { TEMP\$omega <- sqrt(.Machine\$double.eps) }
GUESS <- c(GUESS,list(complexify.ctmm(MLE,"omega",TEMP)))
}
else if("tau position" %in% VARS && level==1) # OU -> OUf (bimodal likelihood)
{ GUESS <- c(GUESS,list(simplify.ctmm(MLE,"diff.tau"))) }

# consider if there is no circulation
if("circle" %in% VARS)
{
Q <- CI["circle",3]
if(is.na(Q) || (Q==Inf) || is.na(IC))
{ GUESS <- c(GUESS,list(simplify.ctmm(MLE,"circle"))) }
}

# consider if eccentricity is zero --- log(major)==log(minor)
if(!CTMM\$isotropic)
{
Q <- c("major","minor")
SD <- ifelse(all(Q %in% CTMM\$features),sqrt(c(GRAD %*% CTMM\$COV[Q,Q] %*% GRAD)),Inf) # variance could collapse early
Q <- log(CTMM\$sigma@par['major']/CTMM\$sigma@par['minor'])
Q <- stats::qnorm(beta/2,mean=Q,sd=SD)
if(is.na(Q) || Q<=0 || is.na(IC))
{ GUESS <- c(GUESS,list(simplify.ctmm(MLE,"minor"))) }
}

# is the animal even moving?
if(!CTMM\$sigma@par['major'] && any(grepl("error",VARS)))
{ GUESS <- c(GUESS,list(simplify.ctmm(MLE,"major"))) }

# consider if we can relax range residence (non-likelihood comparison only)
if(CTMM\$range && length(CTMM\$tau) && (is.na(IC) || CV))
{ GUESS <- c(GUESS,list(simplify.ctmm(MLE,"range"))) }
else if(!CTMM\$range && (is.na(IC) || CV) && TARGET\$tau[1]<Inf)
{ GUESS <- c(GUESS,list(complexify.ctmm(MLE,"range",TARGET))) }

# consider if time link could be more or less detailed

# What was special about these cases?
# consider if the mean could be more detailed or less detailed
REFINE <- c(drift.complexify(MLE),drift.simplify(MLE))

# consider a bunch of new models and update best model without duplication
iterate(GUESS,REFINE,phase=2)
}

# remove any input candidate models
N <- length(MODELS) - CAND
MODELS <- MODELS[1:N]

# return the best or return the full list of models
if(verbose)
{
if(N>1) { MODELS <- sort.ctmm(MODELS,IC=IC,MSPE=MSPE) }
NAMES <- sapply(MODELS,name.ctmm) -> names(MODELS)

if(recurse) # store all model attempts in attribute
{ attr(MODELS,'attempted') <- TRYS }
else # finishing stuff
{
# remove redundant models if outer most run
if(N>1)
{
IN <- rep(TRUE,N)
for(i in 2:N) { if(NAMES[i] %in% NAMES[1:(i-1)]) { IN[i] <- FALSE } }
MODELS <- MODELS[IN]
}

# store model selection criterion
attr(MODELS,"IC") <- IC
attr(MODELS,"MSPE") <- MSPE
}

return(MODELS)
} # end verbose
else
{ return(CTMM) }
}

################
name.ctmm <- function(CTMM,whole=TRUE)
{
if(is.null(CTMM)) { return(NULL) }

FEATURES <- CTMM\$features

# base model
tau <- CTMM\$tau
if(length(tau)==2)
{
if(tau[1]==Inf) { NAME <- "IOU" }
else if(tau[1]>tau[2]) { NAME <- "OUF" }
else if(CTMM\$omega) { NAME <- "OU\u03A9" } # underdamped
else { NAME <- "OUf" } # identical timescales
}
else if(length(tau)==1)
{ if(tau[1]<Inf) { NAME <- "OU" } else { NAME <- "BM" } }
else if(length(tau)==0)
{
POS <- "sigma" %in% names(CTMM) && CTMM\$sigma@par['major']<=0
if(POS || "major" %in% FEATURES || "sigma" %nin% names(CTMM))
{ NAME <- "IID" }
else
{ NAME <- "inactive" }
}

# isotropy
isotropic <- CTMM\$isotropic
if(length(CTMM\$axes)>1 && !is.null(isotropic))
{
if(all(!isotropic))
{ NAME <- c(NAME,"anisotropic") }
else if(!isotropic[1])
{ NAME <- c(NAME,"anisotropic-\u03C3") }
else if(length(isotropic)>1 && !isotropic[2])
{ NAME <- c(NAME,"anisotropic-\u03BC") }
# else
# { NAME <- c(NAME,"isotropic") }
}

# circulation
if(CTMM\$circle || "circle" %in% FEATURES)
{ NAME <- c(NAME,"circulation") }

# error
if(any(CTMM\$error>0) || any(grepl("error",FEATURES)))
{ NAME <- c(NAME,"error") }

{
NAME <- c(NAME,fn(CTMM))
}

# mean
drift <- get(CTMM\$mean)
DNAME <- drift@name(CTMM)

NAME <- paste(NAME,sep="",collapse=" ")

if(whole && !is.null(DNAME))
{ NAME <- paste(NAME,DNAME) }
else if(!whole)
{
if(is.null(DNAME)) { DNAME <- "stationary" }
NAME <- c(NAME,DNAME)
}

return(NAME)
}

########
sort.ctmm <- function(x,decreasing=FALSE,IC="AICc",MSPE="position",flatten=TRUE,INF=FALSE,...)
{
CV <- c("LOOCV","HSCV")
CV <- IC %in% CV

if(is.na(MSPE))
{ ICS <- sapply(x,function(m){get.IC(m,IC)}) }
else if(is.na(IC))
{ ICS <- sapply(x,function(m){get.MSPE(m,MSPE)}) }

if(is.na(MSPE) || is.na(IC))
{
ICS <- nant(ICS,Inf)

IND <- sort(ICS,index.return=TRUE,decreasing=decreasing)\$ix
x <- x[IND]

if(flatten) { return(x) }

# structure the same as below
if(is.na(IC)) { x <- list(x) }
if(is.na(MSPE)) { x <- lapply(x,list) }

return(x)
}

# model type names
NAMES <- sapply(x,function(fit) name.ctmm(fit,whole=FALSE) )
ACOV <- NAMES[1,]
MEAN <- NAMES[2,]

# group by ACOV
ACOVS <- unique(ACOV)
MEANS <- unique(MEAN)

# partition into ACF-identical blocks for MSPE sorting, and then all-identical blocks for likelihood sorting
y <- list()
ICS <- numeric(length(ACOVS)) # ICs of best MSPE models
for(i in 1:length(ACOVS))
{
# ACF-identical block
SUB <- (ACOV==ACOVS[i])
y[[i]] <- x[SUB]
MEAN.SUB <- MEAN[SUB]
MEANS.SUB <- unique(MEAN.SUB)

z <- list()
MSPES <- numeric(length(MEANS.SUB))
for(j in 1:length(MEANS.SUB)) # sort exactly same models by IC
{
# all-identical block
SUB <- (MEAN.SUB==MEANS.SUB[j])
z[[j]] <- sort.ctmm(y[[i]][SUB],IC=IC,MSPE=NA)
MSPES[j] <- get.MSPE(z[[j]][[1]],MSPE) # associate block with best's MSPE
}
IND <- sort(MSPES,index.return=TRUE,decreasing=decreasing)\$ix
y[[i]] <- do.call(c,z[IND]) # flatten to ACF blocks
ICS[i] <- get.IC(y[[i]][[1]],IC) # associate block with best's IC
}
# sort blocks by IC and flatten
IND <- sort(ICS,index.return=TRUE,decreasing=decreasing)\$ix
y <- y[IND]

# BM/IOU log-likelihood is infinitely lower than OU/OUF log-likelihood
RANGE <- sapply(y,function(Y){Y[[1]]\$range})
if(!is.na(IC) && !CV && any(RANGE) && any(!RANGE))
{
if(INF) { for(i in which(!RANGE)) { for(j in 1:length(y[[i]])) { y[[i]][[j]][[IC]] <- Inf } } }
y <- c( y[RANGE] , y[!RANGE] )
}

if(flatten) { y <- do.call(c,y) }

return(y)
}

############
min.ctmm <- function(x,IC="AICc",MSPE="position",...)
{ sort.ctmm(x,IC=IC,MSPE=MSPE,...)[[1]] }

########
summary.ctmm.list <- function(object, IC=NULL, MSPE=NULL, units=TRUE, ...)
{
# use what ctmm.select used
if(is.null(IC)) { IC <- attr(object,"IC") }
if(is.null(MSPE)) { MSPE <- attr(object,"MSPE") }

# defaults for manual lists
if(is.null(IC)) { IC <- "AICc" }
if(is.null(MSPE)) { MSPE <- "position" }

CV <- c("LOOCV","HSCV")
IC <- match.arg(IC,c("AICc","AIC","BIC",CV,NA))
MSPE <- match.arg(MSPE,c("position","velocity",NA))
CV <- IC %in% CV

N <- length(object)
object <- sort.ctmm(object,IC=IC,MSPE=MSPE,flatten=FALSE,INF=TRUE)
M <- length(object)

# if(N==M) { MSPE <- NA } # don't need to sort MSPE
# if(M==1) { IC <- NA } # don't need to sort IC
object <- do.call(c,object)

if(!is.na(IC))
{
ICS <- sapply(object,function(m){get.IC(m,IC)})

# show relative IC
ICS <- ICS - ICS[1]
ICS <- cbind(ICS)
colnames(ICS) <- paste0("\u0394",IC)
}
else { ICS <- NULL }

if(!is.na(MSPE))
{
MSPES <- sapply(object,function(m){get.MSPE(m,MSPE)})

# convert to meters/kilometers
CNAME <- paste0("\u0394","RMSPE")
MIN <- which.min(MSPES)
MSPES <- sqrt(MSPES)
if(MSPES[1]<Inf) { MSPES <- MSPES - MSPES[MIN] }
# MIN <- min(c(MSPES[MSPES>0],Inf))
AVE <- stats::median(MSPES[MSPES>0])
if(MSPE=="position") { UNIT <- "length" } else { UNIT <- "speed" }
UNIT <- unit(AVE,UNIT,concise=TRUE,SI=!units)
MSPES <- MSPES/UNIT\$scale
CNAME <- paste0(CNAME," (",UNIT\$name,")")
MSPES <- cbind(MSPES)
colnames(MSPES) <- CNAME
}
else { MSPES <- NULL }

ICS <- cbind(ICS,MSPES)
rownames(ICS) <- names(object)

if(is.na(MSPE))
{
DOF <- sapply(object,DOF.mean)
DOF <- cbind(DOF)
colnames(DOF) <- "DOF[mean]"
}
else if(MSPE=="position")
{
DOF <- sapply(object,DOF.area)
DOF <- cbind(DOF)
colnames(DOF) <- "DOF[area]"
}
else if(MSPE=="velocity")
{
DOF <- sapply(object,DOF.speed)
DOF <- cbind(DOF)
colnames(DOF) <- "DOF[speed]"
}

METH <- sapply(object,function(m){m\$method})
if(FALSE) # only prints correctly in unicode locale (Windows R bug)
{
DOF <- data.frame(DOF,METH)
colnames(DOF) <- c("DOF[mean]","method")
}

ICS <- cbind(ICS,DOF)

return(ICS)
}
```

## Try the ctmm package in your browser

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

ctmm documentation built on Sept. 24, 2023, 1:06 a.m.