###################################################################################
# Helper functions to create and work with a DLM model sensu Zuur
# x(t)=x(t-1) + w(t), W~MVN(0,1)
# y(t)=Z x(t) + A(t) + D(t) d(t) + v(t), V~MVN(0,R)
# x(t0) = x0 + l, L ~ MVN(0,5)
# MARSS.dfa: The conversion functions have 2 parts
# Part 1 Set up the DLM model in MARSS.marxss form
# Part 2 Call MARSS.marxss to finish the set-up and checking
# Functions that are called by the generic functions print.marssMLE, residuals.marssMLE,
# predict.marssMLE, coef.marssMLE, MARSSinits.marssMLE
# print_dfa, residuals_dfa, predict_dfa, coef_dfa, MARSSinits_dfa
###################################################################################
MARSS.dfa=function(MARSS.call){
# MARSS(data, model=list(), covariates=NULL, z.score=TRUE, demean=TRUE, control=list())
# model.defaults =list(A="zero", R="diagonal and equal", D="zero", x0="zero", V0=diag(5,1), tinitx=0, diffuse=FALSE, m=1)
#Part 1 Set up defaults and check that what the user passed in is allowed
# 1 Check for form dependent user inputs for method and reset defaults for inits, MCbounds, and control if desired
# 2 Specify the text shortcuts and whether factors or matrices can be passed in
# The names in the allowed list do not need to be A, B, Q .... as used in the marss form
# Other names can be used if you want the user to use those names; then in the MARSS.form function,
# you convert the user passed in names into the marss form with the A, B, Q, R, ... names
# checkModelList() will check what the user passes in against these allowed values, so
# so you need to make sure each name in model.defaults has a model.allowed value here
#Check that no args were passed into MARSS that are not allowed
dfa.allowed.in.MARSS.call=c("model","z.score","demean","covariates")
allowed.in.call=c(dfa.allowed.in.MARSS.call,common.allowed.in.MARSS.call)
if(any(!(names(MARSS.call) %in% allowed.in.call))){
bad.names=names(MARSS.call)[!(names(MARSS.call) %in% allowed.in.call)]
msg=paste("Argument ", paste(bad.names, collapse=", ")," not allowed MARSS call for form ", MARSS.call$form, ". See ?MARSS.dfa\n", sep="")
cat("\n","Errors were caught in MARSS.dfa \n", msg, sep="")
stop("Stopped in MARSS.dfa() due to problem(s) with model specification.\n", call.=FALSE)
}
#Set up some defaults
if(is.null(MARSS.call[["z.score"]])) MARSS.call$z.score=TRUE
if(is.null(MARSS.call[["demean"]])) MARSS.call$demean=TRUE
#Start error checking
problem=FALSE
msg=c()
#check that data and covariates elements are matrix or vector, no dataframes, and is numeric
for(el in c("data","covariates"[!is.null(MARSS.call[["covariates"]])])){
if( !(is.matrix(MARSS.call[[el]]) || is.vector(MARSS.call[[el]])) ) {
problem=TRUE
msg = c(msg, paste(el," must be a matrix or vector (not a data frame or 3D array).\n"))
}else{
if(is.vector(MARSS.call[[el]])) MARSS.call[[el]]=matrix(MARSS.call[[el]],1)
}
}
if(!is.null(MARSS.call[["covariates"]])){
if(dim(MARSS.call[["data"]])[2]!= dim(MARSS.call[["covariates"]])[2]){
problem=TRUE
msg = c(msg, "data and covariates must have the same number of time steps.\n")
}
}
if(!is.null(MARSS.call[["covariates"]])){
if(any(is.na(MARSS.call[["covariates"]]))){
problem=TRUE
msg = c(msg, "covariates cannot have any missing values in a standard DFA.\n See User Guide section on DFA for alternate approaches when covariates have missing values.\n")
}
}
if(!is.null(MARSS.call[["model"]][["m"]])){
if(length(MARSS.call[["model"]][["m"]])!=1){
problem=TRUE
msg = c(msg, "model$m must be an integer between 1 and n.\n")
}else{
if( !is.numeric(MARSS.call[["model"]][["m"]]) | !is.wholenumber(MARSS.call$model$m) ){
problem=TRUE
msg = c(msg, "model$m must be an integer between 1 and n.\n")
}else{
if(MARSS.call[["model"]][["m"]] > dim(MARSS.call[["data"]])[1]){
problem=TRUE
msg = c(msg, "model$m must be an integer between 1 and n.\n")
}
}
} }
if(problem){
cat("\n","Errors were caught in MARSS.dfa \n", msg, sep="")
stop("Stopped in MARSS.dfa() due to specification problem(s).\n", call.=FALSE)
}
n=dim(MARSS.call[["data"]])[1]
model.allowed = list(
#if it is a length 1 vector then the value must be one of these. All elements in your model list must be here
A=c("unequal","zero"),
R=c("identity", "zero", "unconstrained", "diagonal and unequal", "diagonal and equal", "equalvarcov"),
D=c("identity", "zero", "unconstrained", "diagonal and unequal", "diagonal and equal", "equalvarcov"),
x0=c("unconstrained", "unequal", "zero"),
B=c("identity","diagonal and equal","diagonal and unequal"),
Q=c("identity","diagonal and equal","diagonal and unequal"),
V0=c("identity", "zero"),
tinitx=c(0,1),
diffuse=c(TRUE,FALSE),
m=1:n,
#This line says what is allowed to be a matrix
matrices = c("Z","A","R","D","x0","V0")
)
if(is.null(MARSS.call[["model"]][["m"]])) m=1 else m=MARSS.call[["model"]][["m"]]
if(!is.null(MARSS.call[["covariates"]])) D="unconstrained" else D="zero"
if(is.null(MARSS.call[["model"]][["Z"]])){ #Set up default Z
Z = matrix(list(), nrow=n, ncol=m)
# insert row (i) & col (j) indices
for(i in seq(n)) {Z[i,] = paste(i, seq(m), sep="")}
# set correct i,j values in Z to numeric 0
if(m>1) for(i in 1:(m-1)){Z[i,(i+1):m] = 0}
}
#defaults for any missing model list elements
model.defaults =list(
A="zero",
R="diagonal and equal",
D=D,
B="identity",
Q="identity",
x0="zero",
V0=diag(5,m),
tinitx=0,
diffuse=FALSE,
m=1,
Z=Z)
#This checks that what user passed in model list can be interpreted and converted to form marss
#if no errors, it updates the model list by filling in missing elements with the defaults
MARSS.call$model=checkModelList( MARSS.call[["model"]], model.defaults, model.allowed )
model=MARSS.call[["model"]]
if(!(MARSS.call$z.score %in% c(TRUE,FALSE)))
stop("Stopped in MARSS.dlm: demean must be TRUE/FALSE.\n", call.=FALSE)
if(!(MARSS.call$demean %in% c(TRUE,FALSE)))
stop("Stopped in MARSS.dlm: demean must be TRUE/FALSE.\n", call.=FALSE)
#Set up U; always 0 for dfa
U=matrix(0,m,1)
#Set up D and d
if(is.null(MARSS.call[["covariates"]])) d=matrix(0,1,1) else d=MARSS.call$covariates
#Set up Q & B; always fixed for dfa
Q=diag(1,m); B=diag(1,m)
# set up list of model components for a marxss model
dfa.model = list(Z=Z, A=model$A, D=model$D, d=d, R=model$R, B=B, U=U, Q=Q, x0=model$x0, V0=model$V0, tinitx=model$tinitx)
dat=MARSS.call[["data"]]
if(MARSS.call$demean){
y.bar = apply(dat, 1, mean, na.rm=TRUE)
dat = (dat - y.bar)
}
if(MARSS.call[["z.score"]]){
Sigma = sqrt(apply(dat, 1, var, na.rm=TRUE))
dat = dat * (1/Sigma)
}
MARSS.call = list(data=dat, inits=MARSS.call$inits, MCbounds=MARSS.call$MCbounds, control=MARSS.call$control, method=MARSS.call$method, form="dlm", silent=MARSS.call$silent, fit=MARSS.call$fit, fun.kf=MARSS.call$fun.kf)
#dfa is a type of marxss model, so use MARSS.marxss to test it and set up the marss object
tmp = MARSS.marxss(list(data=dat,model=dfa.model,method=MARSS.call$method,silent=MARSS.call$silent))
#marss is the name for the form=marss model object that MARSS.form functions return
#need to add "dfa" to attribute form
marxss_object=tmp$model
attr(marxss_object, "form")=c("dfa", "marxss")
MARSS.call$model = marxss_object
MARSS.call$marss = tmp$marss
## Return MARSS inputs as list
MARSS.call
}
#This works since dfa just creates a marxss object so x$model is marssMODEL form=c("dfa", "marxss")
#probably want to customize for dfa later
print_dfa = function(x){ return(print_marxss(x)) }
coef_dfa = function(x){ return(coef_marxss(x)) }
MARSSinits_dfa = function(MLEobj, inits){ return(MARSSinits_marxss(MLEobj, inits)) }
predict_dfa = function(x, newdata, n.ahead, t.start){ predict_marxss(x, newdata, n.ahead, t.start) }
describe_dfa = function(modelObj){ describe_marss(modelObj) }
is.marssMODEL_dfa = function(modelObj, method="kem"){ is.marssMODEL_marxss(modelObj, method=method) }
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.