###################################################################################
# marss form definition file
###################################################################################
###################################################################################
# MARSS.marss create marss from MARSS inputs
# marss_to_marss not needed
# print_marss
# coef_marss
# MARSSintis_marss
# predict_marss
# describe_marss
# is.marssMODEL_marss
###################################################################################
MARSS.marss=function(MARSS.call){
#Part 1 Set up defaults and check that what the user passed in is allowed
#Check that no args were passed into MARSS that are not allowed
marss.allowed.in.MARSS.call=c("model")
allowed.in.call=c(marss.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.marxss\n", sep="")
cat("\n","Errors were caught in MARSS.marss \n", msg, sep="")
stop("Stopped in MARSS.marss() due to problem(s) with model specification.\n", call.=FALSE)
}
# 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 form=marss object
# 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 form marss names 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
model.allowed = list(
A=c("scaling", "unconstrained", "unequal", "equal", "zero"),
B=c("identity", "zero", "unconstrained", "unequal", "diagonal and unequal", "diagonal and equal", "equalvarcov"),
Q=c("identity", "zero", "unconstrained", "diagonal and unequal", "diagonal and equal", "equalvarcov"),
R=c("identity", "zero", "unconstrained", "diagonal and unequal", "diagonal and equal", "equalvarcov"),
U=c("unconstrained", "equal", "unequal", "zero"),
x0=c("unconstrained", "equal", "unequal", "zero","diagonal and unequal","diagonal and equal"),
V0=c("identity", "zero", "unconstrained", "diagonal and unequal", "diagonal and equal", "equalvarcov"),
Z=c("identity","unconstrained", "diagonal and unequal", "diagonal and equal", "equalvarcov", "onestate"),
tinitx=c(0,1),
diffuse=c(TRUE,FALSE),
factors = c("Z"),
matrices = c("A","B","Q","R","U","x0","Z","V0")
)
#model.defaults is form dependent so you must specify it
model.defaults =list(Z="identity", A="scaling", R="diagonal and equal", B="identity", U="unconstrained", Q="diagonal and unequal", x0="unconstrained", V0="zero", tinitx=0, diffuse=FALSE)
#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 )
# Part 2 Convert the model list to a marssMODEL object, form=marxss
## set up fixed and free elements
fixed = free = list()
model = MARSS.call$model
model.elem = c("Z","A","R","B","U","Q","x0","V0")
dat = MARSS.call[["data"]]
if(is.vector(dat)) dat=matrix(dat,1,length(dat))
n = dim(dat)[1]; TT = dim(dat)[2]
if(is.null(rownames(dat))){
Y.names = paste("Y",seq(1, n),sep="") #paste(seq(1, n), sep="") #
rownames(dat)=Y.names
}else{
Y.names = rownames(dat)
if(any(duplicated(Y.names))){
for(i in Y.names[duplicated(Y.names)]){
loc = (Y.names==i)
nn = sum(loc)
Y.names[loc]=paste(i,"-",1:nn,sep="")
rownames(dat)[loc] = Y.names[loc]
MARSS.call[["data"]] = dat
}
}
}
## Set m based on Z specification IF Z was specified; errors will be reported later if m conflicts with other parameters
m = NA
if (identical(model$Z, "unconstrained")) m = n
if (identical(model$Z, "equalvarcov")) m = n
if (identical(model$Z, "diagonal and equal")) m = n
if (identical(model$Z, "diagonal and unequal")) m = n
if (identical(model$Z, "onestate")) m = 1
if (identical(model$Z, "identity")) m = n
if (is.factor(model$Z)) m = length(levels(model$Z))
if (is.array(model$Z)) m = dim(model$Z)[2]
X.names=NULL
if(!(is.null(model[["X.names"]]))) X.names = model[["X.names"]]
if(is.null(X.names) & identical(model$Z,"identity"))
X.names=paste("X.",Y.names,sep="")
if( is.null(X.names) & is.array(model$Z)){
if(length(dim(model$Z))==3)
if(dim(model$Z)[3]==1)
if(is.design(model$Z) & !is.null(colnames(model$Z)))
X.names=colnames(model$Z)
}
if( is.null(X.names) & is.factor(model$Z)){
X.names=unique(as.character(model$Z))
}
if(is.null(X.names) & is.matrix(model$Z))
if(is.design(model$Z) & !is.null(colnames(model$Z)))X.names=colnames(model$Z)
if(is.null(X.names) & is.matrix(model$Z))
if(is.identity(model$Z))
X.names = paste("X.",Y.names,sep="")
if(is.null(X.names)) X.names = paste("X",seq(1, m),sep="") #paste(seq(1, m),sep="") #
#3rd dim of params are set to 1 and will be reset to correct value at end
model.dims = list(data=c(n,TT),x=c(m,TT),y=c(n,TT),w=c(m,TT),v=c(n,TT),Z=c(n,m,1),U=c(m,1,1),A=c(n,1,1),B=c(m,m,1),Q=c(m,m,1),R=c(n,n,1),x0=c(m,1,1),V0=c(m,m,1))
## Error-checking section that is specific to marss form
# Note most error checking happens in checkMARSSInputs, checkModelList, and is.marssMLE
#If a model elements passed in as a factors, make sure it is the correct length otherwise construction of marssMODEL object will break
problem = FALSE
msg=NULL
## Check model structures that are passed in as factors
correct.factor.len = list(Z=n,A=n,R=n,B=m,U=m,Q=m,x0=m,V0=m)
for (el in model.elem) {
#if a factor then it needs to have the correct length otherwise construction of marssMODEL object will break and no NAs
if( is.factor(model[[el]]) ) {
if(length(model[[el]]) != correct.factor.len[[el]]) {
problem=TRUE
msg = c(msg, paste(" The model$", el, " is being passed as a factor and should be length ", correct.factor.len[[el]], " based on data dims. It's not. See help file.\n", sep="")) }
if(NA.in.fac <- NA %in% model[[el]]) {
problem=TRUE
msg = c(msg, paste(" NAs are not allowed in model factor for ", el, ". See help file.\n", sep="")) }
} #is factor
} # end for (el in model.elem)
#if el == Z then factor needs to have m levels
if( is.factor(model$Z) ) {
if(length(levels(model$Z)) != m) {
problem=TRUE
msg=c(msg," When Z is a factor, the number of levels must equal the number of state processes (m).\n")
} }
#Check that if A is scaling, then Z spec must lead to a design matrix
if(identical(model$A,"scaling")){
if(is.array(model$Z) & length(dim(model$Z))==3)
if(dim(model$Z)[3]!=1 && !is.design(model$Z, zero.cols.ok = TRUE)) { #if it is a matrix
problem=TRUE
msg = c(msg, " If A is scaling(the default), then Z must be a time-constant design matrix:(0,1) and rowsums=1.\nYou can construct a scaling A matrix and pass that in.")
}
if((!is.array(model$Z) & !is.factor(model$Z))) #if it is a string
if(!(model$Z %in% c("onestate","identity"))){
problem=TRUE
msg = c(msg, " If A is scaling(the default), then Z must be a time-constant design matrix:(0,1) and rowsums=1.\n")
}
if(is.matrix(model$Z) && !is.design(model$Z, zero.cols.ok = TRUE)) { #if it is a matrix, won't be array due to first test
problem=TRUE
msg = c(msg, " If A is scaling(the default), then Z must be a time-constant design matrix:(0,1) and rowsums=1.\n")
}
}
if(is.array(model$x0) & length(dim(model$x0))==3)
if(dim(model$x0)[3]!=1){
problem=TRUE
msg = c(msg, " x0 cannot be time-varying thus if x0 in model arg is 3D, the 3rd dim must equal 1.\n")
}
if(is.array(model$V0) & length(dim(model$V0))==3)
if(dim(model$V0)[3]!=1){
problem=TRUE
msg = c(msg, " V0 cannot be time-varying thus if V0 in model arg is 3D, the 3rd dim must equal 1.\n")
}
#If there are problems
if(problem) {
cat("\n","Errors were caught in MARSS.marss \n", msg, sep="")
stop("Stopped in MARSS.marss() due to problem(s) with model specification.\n", call.=FALSE)
}
#end of error section
##Translate the text shortcuts into a marssMODEL object
## Translate the model structure names (shortcuts) into fixed and free
## fixed is a dim(1)*dim(2) X 1 vector of the fixed (intercepts) values
## free is a dim(1)*dim(2) X p vector of the free (betas) values for the p estimated elements
model.elem = c("Z","A","R","B","U","Q","x0","V0")
if(which(model.elem=="Z")>which(model.elem=="A")) model.elem=rev(model.elem) #Z must go first
tmp=list()
for(el in model.elem) {
tmp[[el]]="not assigned"
if(el=="Z" & is.factor(model$Z)) {
tmp[[el]] = matrix(0,model.dims$Z[1], model.dims$Z[2])
for(i in X.names) tmp[[el]][which(model$Z==i), which(as.vector(X.names)==i)] = 1
}
if(el=="Z" & identical(model$Z,"onestate")) { #m=1
tmp[[el]] = matrix(1, n, 0)
}
if( identical(model[[el]],"identity") ) {
tmp[[el]] = diag(1,model.dims[[el]][1])
}
if(identical(model[[el]],"diagonal and equal")) {
tmp[[el]] = array(list(0),dim=c(model.dims[[el]][1],model.dims[[el]][2]))
diag(tmp[[el]])="diag" #paste(el,"(diag)",sep="")
if(length(tmp[[el]])==1) tmp[[el]][1,1]=el
}
if(identical(model[[el]],"diagonal and unequal")) {
tmp[[el]] = array(list(0),dim=c(model.dims[[el]][1],model.dims[[el]][2]))
dim.mat = model.dims[[el]][1]
el.labs=as.character(1:dim.mat)
if(el %in% c("V0","Q","B")) el.labs=X.names
if(el %in% c("Z","R")) el.labs=Y.names
diag(tmp[[el]])=paste("(",el.labs,",", el.labs,")",sep="") #paste(el,"(",as.character(1:dim.mat),",",as.character(1:dim.mat),")",sep="")
if(length(tmp[[el]])==1) tmp[[el]][1,1]=el
}
if(identical(model[[el]],"unconstrained") | identical(model[[el]],"unequal")){
tmp[[el]]=array(NA,dim=c(model.dims[[el]][1],model.dims[[el]][2]))
if(el %in% c("Q","R","V0")){ #variance-covariance matrices
dim.mat = model.dims[[el]][1]
for(i in 1:dim.mat){
for(j in 1:dim.mat) tmp[[el]][i,j]=tmp[[el]][j,i]=paste("(",i,",",j,")",sep="") #paste(el,"(",i,",",j,")",sep="")
}
}else{ #not var-cov matrix
row.name=1:model.dims[[el]][1]
col.name=1:model.dims[[el]][2]
if(el %in% c("U","x0")) row.name=X.names
if(el == "A") row.name=Y.names
for(i in 1:model.dims[[el]][1]){
for(j in 1:model.dims[[el]][2]){
if(model.dims[[el]][2]>1) tmp[[el]][i,j]=paste("(",row.name[i],",",col.name[j],")",sep="") #paste(el,"(",row.name[i],",",col.name[j],")",sep="")
else tmp[[el]][i,j]=paste(row.name[i],sep=",") #paste(el,row.name[i],sep=",")
}
}
}
if(length(tmp[[el]])==1) tmp[[el]][1,1]=el
} #unconstrained
if(identical(model[[el]],"equalvarcov")) {
tmp[[el]]=array("offdiag",dim=c(model.dims[[el]][1],model.dims[[el]][2])) #array(paste(el,"(offdiag)",sep=""),dim=model.dims[[el]])
diag(tmp[[el]])="diag" #paste(el,"(diag)",sep="")
if(length(tmp[[el]])==1) tmp[[el]][1,1]=el
}
if(identical(model[[el]],"equal")) {
tmp[[el]]=array("1",dim=c(model.dims[[el]][1],model.dims[[el]][2])) #array(el,dim=model.dims[[el]])
}
if(identical(model[[el]],"zero")) {
tmp[[el]]=array(0,dim=c(model.dims[[el]][1],model.dims[[el]][2]))
}
if(is.array(model[[el]])) {
tmp[[el]]=model[[el]]
}
if(el=="A" & identical(model[[el]], "scaling")) { #check above ensures that Z is design and time-invariant
## Construct A from fixed Z matrix
tmp[[el]] = matrix(list(),model.dims$A[1],model.dims$A[2])
tmp[[el]][,1]=Y.names
for(i in 1:m) {
tmp[[el]][min(which(tmp$Z[,i]!=0)), 1] = 0
}
}
if(identical(tmp[[el]],"not assigned")) stop(paste("MARSS.marxss: tmp was not assigned for ",el,".\n",sep=""))
free[[el]] = convert.model.mat(tmp[[el]])$free
fixed[[el]] = convert.model.mat(tmp[[el]])$fixed
#set the last dim of the model.dims since it was at a temp value to start
model.dims[[el]][3]=max(dim(free[[el]])[3],dim(fixed[[el]])[3])
}
#Set the marssMODEL form marss
#This is the f+Dp form for the MARSS model used for user displays, printing and such
marss_object = list(fixed=fixed, free=free, data=dat, tinitx=model$tinitx, diffuse=model$diffuse)
#set the attributes
class(marss_object) = "marssMODEL"
attr(marss_object, "obj.elements")=c("fixed","free","data","tinitx","diffuse")
attr(marss_object, "form")="marss"
attr(marss_object, "model.dims")=model.dims
#par.names are what needs to be in fixed/free pair
attr(marss_object, "par.names")=c("Z","A","R","B","U","Q","x0","V0")
attr(marss_object, "X.names")=X.names
attr(marss_object, "Y.names")=Y.names
attr(marss_object, "equation")="x_{t}=B_{t}*x_{t-1}+U_{t}+w_{t}; w_{t}~MVN(0,Q_{t})\ny_{t}=Z_{t}*x_{t}+A_{t}+v_{t}; v_{t}~MVN(0,R_{t})"
#Put the marss model into model and marss
MARSS.call$model = marss_object
MARSS.call$marss=marss_object
## Return MARSS call list with $marss and $model added
MARSS.call
}
#the par element of a marssMLE object is already in form=marss. No change needed if form is marss
print_marss = function(x){ return(x) }
#the par element of a marssMLE object is already in form=marss. No change needed
coef_marss = function(x){ return(x) }
#MARSSinits needs a valid $par element with inits were marss form
#so no need to change inits
MARSSinits_marss = function(MLEobj, inits){ return(inits) }
predict_marss = function(x, newdata, n.ahead, t.start){
#the predict function takes the marssMLE obj along with newdata and returns a
#marssMODEL (form=marss) object constructed using info in newdata that is ready for use in prediction
#x means the x (marssMLE object) for prediction
#a marss marssMODEL object can be made into marxss, so convert to that and then use the predice_marxss function
x.marxss = marss_to_marxss(x, C.and.D.are.zero=TRUE)
return(predict_marxss(x.marxss, newdata, n.ahead, t.start))
}
#describe_marss works generally with marxss form models (of which marss is one type)
describe_marss = function(modelObj, model.elem=NULL){
#This returns the structure of a model using text strings; works for form=marss or marxss
#You can pass in model.elem to have is drop some elements of marxss (if you don't want these to print)
fixed = modelObj$fixed
free = modelObj$free
model.dims=attr(modelObj,"model.dims")
m=model.dims$x[1]; n=model.dims$y[1]
model.elem=attr(modelObj,"par.names")
#set up the constr type list with an element for each par name and init val of "none"
constr.type=vector("list",length(model.elem))
names(constr.type)=model.elem
constr.type[model.elem]="none"
for(elem in model.elem){
TT.max = max( dim(fixed[[elem]])[3], dim(free[[elem]])[3])
if(TT.max > 1) constr.type[[elem]]="time-varying"
}
# For everything below, elem is time constant however tests are written to take a 3D matrix and it is not req that t=1
############ Check the other matrices
# This part determines the form of free/fixed
for(elem in model.elem){
dimm = model.dims[[elem]][1]
dimc = model.dims[[elem]][2]
while(identical(constr.type[[elem]],"none")) {
#Z is not time varying and is a design matrix; note test above would have ensured Z is time constant
if( elem=="Z" & is.fixed(free[[elem]]) & dim(fixed[[elem]])[3]==1 ){
if( all(apply(fixed[[elem]],3,is.design,dim=c(dimm,dimc))) ) {
Z.names = c()
X.names=attr(modelObj,"X.names")
if(is.null(X.names)) X.names=as.character(1:m)
for(i in 1:n) Z.names = c(Z.names, X.names[as.logical(fixed$Z[i,,1])]) #fixed$Z will be time constant
constr.type[[elem]] = paste("design matrix with rows:",paste(Z.names,collapse=" "))
break }} #break out of the while for elem
if(is.fixed(free[[elem]])) { #it's fixed
if(all(fixed[[elem]]==0)) { constr.type[[elem]] = "fixed and zero"; break }
if(all(fixed[[elem]]==1)) { constr.type[[elem]] = "fixed and all one"; break }
if(length(unique(fixed[[elem]]))==1)
{ constr.type[[elem]] = paste("fixed and all",fixed[[elem]][1]); break }
if( all(apply(fixed[[elem]],3,is.identity,dim=c(dimm,dimc)))) constr.type[[elem]] = "identity"
else constr.type[[elem]] = "fixed"
break #break out of while
}
if(length(free[[elem]])==1) { constr.type[[elem]] = "scalar"; break }
if(length(unique(as.vector(free[[elem]])))==1 & dim(free[[elem]])[2]==1) { constr.type[[elem]] = "all equal"; break }
if(all(apply(free[[elem]],3,is.identity))) { constr.type[[elem]] = "unconstrained"; break }
tmp.free = free[[elem]]
#The next parts are only for parameters that are not column vectors, so dim 2 !=1
#fixed.free.to.formula requires that 3D mats have dim3=1, which will be true here since check above that time constant
#creates a list matrix version of model
if(model.dims[[elem]][2]!=1){
tmp.mat = fixed.free.to.formula( fixed[[elem]],free[[elem]],c(dimm,dimc) )
if(elem %in% c("Q","R","V0")){ #variance-covariance matrices
dimm=sqrt(dim(fixed[[elem]])[1])
if(
all(fixed[[elem]]==0) &
all(unlist(tmp.mat[upper.tri(tmp.mat)])==unlist(t(tmp.mat)[upper.tri(tmp.mat)])) &
is.design(free[[elem]]) &
length(unique(unlist(tmp.mat[upper.tri(tmp.mat)])))==dimm*(dimm-1)/2 &
length(unique(unlist(diag(tmp.mat))))==dimm )
{ constr.type[[elem]] = "unconstrained"; break }
}
if(is.diagonal(tmp.mat)) {
tmp.diag = unique(diag(tmp.mat)) #it's a list
tmp.diag.est=tmp.diag[sapply(tmp.diag,is.character)]
tmp.diag.fixed=tmp.diag[sapply(tmp.diag,is.numeric)]
if(length(tmp.diag.est)==dimm) { constr.type[[elem]] = "diagonal and unequal"; break }
if(length(tmp.diag.est)==1 & length(tmp.diag.fixed)==0) { constr.type[[elem]] = "diagonal and equal"; break }
if(length(tmp.diag.est)>1 & length(tmp.diag.fixed)==0) {
constr.type[[elem]] = paste("diagonal with ",length(tmp.diag.est)," groups",sep="") ; break }
if(length(tmp.diag.est)>=1 & length(tmp.diag.fixed)>0) {
constr.type[[elem]] = paste("diagonal with fixed elements and ",length(tmp.diag.est)," estimated values",sep="") ; break }
}
if(is.equaltri(tmp.mat)) {
if(elem %in% c("Q","R","V0")){ #variance-covariance matrices
constr.type[[elem]] = "one variance value and covariance value"; break
}else{ constr.type[[elem]] = "one diagonal value and one off-diagonal value"; break }
}
if(is.blockdiag(tmp.mat)) {
if(elem %in% c("Q","R","V0")){ #variance-covariance matrices
constr.type[[elem]] = "variance-covariance matrix with block diagonal structure"; break
}else{ constr.type[[elem]] = "block diagonal matrix"; break }
}
}
constr.type[[elem]]="see summary()" #not assigned to one of the above cases
}
} #model.elem
for(elem in model.elem){
constr.type[[elem]] = paste(constr.type[[elem]]," (",paste(model.dims[[elem]][1:2],collapse=" x "),")",sep="")
}
return(constr.type)
}
########################################################################
# is.marssMODEL_marss function
# Check that the marss object has all the parts it needs
# fixed, free, and par.names
# and that these have the proper size and form
# m is pulled from fixed$x0
########################################################################
is.marssMODEL_marss <- function(modelObj, method="kem"){
msg=NULL
## Set up par.names that should be in a marss model
en = c("Z", "A", "R", "B", "U", "Q", "x0", "V0")
#Check that par.names has these and only these names
par.names=attr(modelObj, "par.names")
if( !all(en %in% par.names ) ){
msg = c(msg, "Element ", en[!(en %in% par.names)], " is missing from the par.names attribute of the model object.\n")
}
if( !all( par.names %in% en ) ) {
msg = c(msg, "Only ", en, "should be in the par.names attribute of the model object.\n")
}
model.dims=attr(modelObj, "model.dims")
if( !all(en %in% names(model.dims) ) ){
msg = c(msg, "Element ", en[!(en %in% names(model.dims) )], " is missing from the model.dims attribute of the model object.\n")
}
if(!is.null(msg)){ #rest of the tests won't work so stop now
return(msg)
}
###########################
# Check model.dims are correct
###########################
n = dim(modelObj$data)[1]
TT = dim(modelObj$data)[2]
m = dim(modelObj$fixed$x0)[1]
en = c("Z", "A", "R", "B", "U", "Q", "x0", "V0", "data", "x", "y", "w", "v")
correct.dim1 = c(Z=n,A=n,R=n,B=m, U=m, Q=m, x0=m, V0=m, data=n, x=m, y=n, w=m, v=n)
correct.dim2 = c(Z=m,A=1,R=n,B=m, U=1, Q=m, x0=1, V0=m, data=TT, x=TT, y=TT, w=TT, v=TT)
for (elem in en) {
## Check for problems in the fixed/free pairs. Problems show up as TRUE
dim.flag1 = dim.flag2 = FALSE
# check dim
dim.flag1 = c(dim.flag1,!( model.dims[[elem]][1] == correct.dim1[[elem]] ))
dim.flag2 = c(dim.flag2,!( model.dims[[elem]][2] == correct.dim2[[elem]] ))
}
if (any(c(dim.flag1, dim.flag2))) { #There's a problem
if(any(dim.flag1)) {
msg = c(msg, paste("Dim 1 of ", en[dim.flag1], "is incorrect. Dims should be ", correct.dim1[dim.flag1], ", for a marss model.\n"))
}
if(any(dim.flag2)) {
msg = c(msg, paste("Dim 2 of ", en[dim.flag2], "is incorrect. Dims should be ", correct.dim2[dim.flag2], ", for a marss model.\n"))
}
msg=c("\nErrors were caught in is.marssMODEL_marss()\n", msg)
return(msg)
}
fixed=modelObj$fixed; free=modelObj$free
###########################
# Check that x0 and V0 are not time-varying
###########################
en = c("x0", "V0")
time.var = NULL
for (elem in en) {
time.var.flag = FALSE
time.var.flag = dim(fixed[[elem]])[3]!=1 | dim(free[[elem]])[3]!=1
time.var <- c(time.var, time.var.flag)
}
if(any(time.var)) { #There's a problem
msg = c(msg, paste(en[time.var], "cannot be time-varying. 3rd dim of fixed and free must equal 1.\n"))
msg=c("\nErrors were caught in is.marssMODEL_marss()\n", msg)
return(msg)
}
###########################
# Check that none of the var-cov matrices have negative values on the diagonal
# and that there are no f+Dq elements only f+0q or 0+Dq
# and D must be a design matrix, so no beta_1*q1 + beta_2*q2 elements
###########################
en = c("R", "Q", "V0")
neg = bad.var = not.design = NULL
for (elem in en) {
neg.flag = bad.var.flag = not.design.flag = FALSE
for(i in 1:max(dim(free[[elem]])[3],dim(fixed[[elem]])[3])){
if(dim(fixed[[elem]])[3]==1){i1=1}else{i1=i}
if(dim(free[[elem]])[3]==1){i2=1}else{i2=i}
if(is.fixed(free[[elem]][,,min(i,dim(free[[elem]])[3]),drop=FALSE])){ #this works on 3D mats
zero.free.rows = matrix(TRUE,correct.dim1[[elem]]*correct.dim2[[elem]],1)
}else{
zero.free.rows=apply(free[[elem]][,,i2,drop=FALSE]==0,1,all) #works on 3D mat
#the requirement is that each estimated element (in p) appears only in one place in the varcov mat, but fixed rows (0 rows) are ok
not.design.flag = !is.design(free[[elem]][,,i2,drop=FALSE],strict=FALSE,zero.rows.ok=TRUE,zero.cols.ok=TRUE) #works on 3D if dim3=1
}
zero.fixed.rows=apply(fixed[[elem]][,,i1,drop=FALSE]==0,1,all) #works on 3D
fixed.mat = unvec(fixed[[elem]][,,i1],dim=c(correct.dim1[[elem]],correct.dim2[[elem]]))
if( any(!zero.fixed.rows & !zero.free.rows) ) bad.var.flag = TRUE #no f+Dq rows
if(any(takediag(fixed.mat)<0,na.rm=TRUE)) neg.flag=TRUE #no negative diagonals
} #end the for loop over time
not.design = c(not.design, not.design.flag)
neg = c(neg, neg.flag)
bad.var = c(bad.var, bad.var.flag)
} #enf the for loop over elem
if(any(neg)) {
msg = c(msg, paste("Negative values are on the diagonal of ", en[neg], ". Neg values are illegal on the diag of a var-cov matrix.\n", sep=""))
}
if(any(bad.var)) {
msg = c(msg, paste("Fixed and estimated values are combined in some elements of ", en[bad.var], ". This is not allowed.\n", sep=""))
}
if(any(not.design)) {
msg = c(msg, paste("The D matrices of ", en[not.design], " must be design matrices.\n", sep=""))
}
###########################
# Check that V0, Q and R matrices are symmetric and positive-definite
###########################
en = c("R", "Q", "V0")
pos = symm = NULL
for (elem in en) {
varcov.flag = TRUE; varcov.msg=""
var.dim = c(correct.dim1[[elem]],correct.dim2[[elem]])
for(i in 1:model.dims[[elem]][3]){
if(dim(fixed[[elem]])[3]==1){i1=1}else{i1=i}
if(dim(free[[elem]])[3]==1){i2=1}else{i2=i}
#works on 3D if dim3=1
par.as.list = fixed.free.to.formula(fixed[[elem]][,,i1,drop=FALSE],free[[elem]][,,i2,drop=FALSE],var.dim) #coverts the fixed,free pair to a list matrix
tmp=is.validvarcov(par.as.list, method=method)
varcov.flag=varcov.flag & tmp$ok
if(!tmp$ok) varcov.msg = c(varcov.msg, paste(" ", tmp$error, "at t=", i, "\n",sep=""))
if(!varcov.flag) msg = c(msg, paste("The variance-covariance matrix ", elem, " is not properly constrained.\n", sep=""), varcov.msg)
} #end for loop over time
} #end for loop over elements
if(length(msg) == 0){ return(NULL)
}else {
msg=c("\nErrors were caught in is.marssMODEL_marss()\n", msg)
return(msg)
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.