R/MARSSkemcheck.r
In nwfsc-timeseries/MARSS: Multivariate Autoregressive State-Space Modeling
Defines functions
MARSSkemcheck
Documented in
MARSSkemcheck
MARSSkemcheck <- function(MLEobj) {
# This checks that the model can be handled by the MARSSkem algorithm
# Most of this is implementing the restrictions in Summary of Requirements for Degenerate Models in derivation
MODELobj <- MLEobj$marss
fixed <- MODELobj$fixed
free <- MODELobj$free
par.dims <- attr(MODELobj, "model.dims")
m <- par.dims[["x"]][1]
n <- par.dims[["y"]][1]
TT <- par.dims[["data"]][2]
pseudolim <- 1E-8
ok <- TRUE
msg <- NULL
# If TT=2 then kf will break
if (TT <= 2) {
msg <- c(msg, "The number of time steps is <=2.\nMore than 2 data points are needed to estimate parameters.\n")
ok <- FALSE
}
############ Check that if fixed, B is within the unit circle
for (t in 1:max(dim(free$B)[3], dim(fixed$B)[3])) {
ifixed <- min(t, dim(fixed$B)[3])
if (is.fixed(free$B[, , ifixed, drop = FALSE])) { # works on 3D matrices
# parmat needs MODELobj and par list
if (is.null(MLEobj$par$B)) tmpparB <- MLEobj$start$B else tmpparB <- MLEobj$par$B
tmp.MLEobj <- list(marss = MODELobj, par = list(B = tmpparB)) # B is fixed but par might have cols from other times
parB <- parmat(tmp.MLEobj, "B", t = t)$B
if (!is.unitcircle(parB)) {
msg <- c(msg, " All the eigenvalues of B must be within the unit circle: all(abs(eigen(B)$values)<=1)\n")
ok <- FALSE
}
}
} # end for over time to check B
############ Check that if R has 0s, then the corresponding row of A and Z are fixed
# See Summary of Requirements for Degenenate Models in EMDerivation.pdf
el <- "R"
# extracts the r=0 rows (or cols)
diag.rows <- 1 + 0:(par.dims[[el]][1] - 1) * (par.dims[[el]][1] + 1)
Tmax <- 0
for (par.test in c("R", "Z", "A")) {
Tmax <- max(Tmax, dim(fixed[[par.test]])[3], dim(free[[par.test]])[3])
}
if (is.null(MLEobj[["par"]])) MLEobj$par <- MLEobj$start
# Run with MARSSkfss because it has more internal checks and error messages
MLEobj$kf <- MARSSkfss(MLEobj)
if (!MLEobj$kf$ok) {
return(list(ok = FALSE, msg = c(MLEobj$kf$errors, msg)))
}
MLEobj$Ey <- MARSShatyt(MLEobj)
msg.tmp <- NULL
for (i in 1:TT) {
ifixed <- min(i, dim(fixed[[el]])[3])
ifree <- min(i, dim(free[[el]])[3])
zero.diags <- is.fixed(free[[el]][diag.rows, , ifree, drop = FALSE], by.row = TRUE) & apply(fixed[[el]][diag.rows, , ifixed, drop = FALSE] == 0, 1, all) # fixed rows and 0
if (any(zero.diags)) {
II0 <- makediag(as.numeric(zero.diags))
if (i <= Tmax) {
for (par.test in c("Z", "A")) {
II <- diag(1, par.dims[[par.test]][2])
ifree.par <- min(i, dim(free[[par.test]])[3])
dpart <- sub3D(free[[par.test]], t = ifree.par)
par.not.fixed <- any(((II %x% II0) %*% dpart) != 0)
if (par.not.fixed) {
if (par.test == "Z") msg <- c(msg, paste("t=", i, ": For method=kem (EM), if an element of the diagonal of R is 0, the corresponding row of Z must be fixed. See MARSSinfo('AZR0').\n", sep = ""))
if (par.test == "A") msg <- c(msg, paste("t=", i, ": For method=kem (EM), if an element of the diagonal of R is 0, the corresponding row of both A and D must be fixed. See MARSSinfo('AZR0').\n", sep = ""))
}
} # for par.test
} # if par needs to be tested; don't need to do over all TT just up to Tmax
Z <- parmat(MLEobj, "Z", t = i)$Z
A <- parmat(MLEobj, "A", t = i)$A
Z.R0 <- Z[zero.diags, , drop = FALSE]
A.R0 <- A[zero.diags, , drop = FALSE]
y.R0 <- MLEobj$Ey$ytT[zero.diags, i, drop = FALSE]
y.resid <- y.R0 - Z.R0 %*% MLEobj$kf$xtT[, i, drop = FALSE] - A.R0
if (any(y.resid > pseudolim)) {
msg.tmp <- c(msg.tmp, paste(" Z, A, and y for the R=0 rows at time=", i, " do not agree. For the R=0 rows, E(y) must equal Z * E(x) + A.\n", sep = ""))
}
} # any zero.diags
} # end for loop over time
if (!is.null(msg)) {
msg <- c(msg, msg.tmp[1:min(10, length(msg.tmp))]) # error msgs could be really long, only print first 10
ok <- FALSE
}
############ Check that II_q^{0} is time constant
# See Summary of Requirements for Degenenate Models in EMDerivation.pdf
Tmax <- 0
for (par.test in c("Q")) {
Tmax <- max(Tmax, dim(fixed[[par.test]])[3], dim(free[[par.test]])[3])
}
II0 <- list()
for (i in 1:Tmax) {
for (el in c("Q")) {
ifixed <- min(i, dim(fixed[[el]])[3])
ifree <- min(i, dim(free[[el]])[3])
diag.rows <- 1 + 0:(par.dims[[el]][1] - 1) * (par.dims[[el]][1] + 1)
zero.diags <- is.fixed(free[[el]][diag.rows, , ifree, drop = FALSE], by.row = TRUE) & apply(fixed[[el]][diag.rows, , ifixed, drop = FALSE] == 0, 1, all) # fixed rows
II0Q <- makediag(as.numeric(zero.diags))
II0Q <- unname(II0[[el]])
}
if (i == 1) {
II0Q.1 <- II0Q
} else {
if (!isTRUE(all.equal(II0Q, II0Q.1))) {
msg <- c(msg, paste("t=", i, ": The placement of 0 variances in Q must be time constant.\n", sep = ""))
ok <- FALSE
}
}
} # end for loop over time
############ Check that if R and Q both have 0s, then the U and Bs are appropriately fixed
# See Summary of Requirements for Degenenate Models in EMDerivation.pdf
Tmax <- 0
for (par.test in c("R", "Q", "U", "B", "Z")) {
Tmax <- max(Tmax, dim(fixed[[par.test]])[3], dim(free[[par.test]])[3])
}
II0 <- list()
for (i in 1:Tmax) {
for (el in c("R", "Q")) {
ifixed <- min(i, dim(fixed[[el]])[3])
ifree <- min(i, dim(free[[el]])[3])
diag.rows <- 1 + 0:(par.dims[[el]][1] - 1) * (par.dims[[el]][1] + 1)
zero.diags <- is.fixed(free[[el]][diag.rows, , ifree, drop = FALSE], by.row = TRUE) & apply(fixed[[el]][diag.rows, , ifixed, drop = FALSE] == 0, 1, all) # fixed rows
II0[[el]] <- makediag(as.numeric(zero.diags))
}
for (el in c("B", "U")) {
ifree <- min(i, dim(free[[el]])[3])
# I don't know what par$Z will be. I want to create a par$Z where there is a non-zero value for any potentially non-zero Z's
tmp.MLEobj <- list(marss = MODELobj, par = list(Z = matrix(1, dim(free$Z)[2], 1)))
tmp.MLEobj$marss$fixed$Z[tmp.MLEobj$marss$fixed$Z != 0] <- 1
tmp.MLEobj$marss$free$Z[tmp.MLEobj$marss$free$Z != 0] <- 1
parZ <- parmat(tmp.MLEobj, "Z", t = i)$Z
II <- diag(1, par.dims[[el]][2])
dpart <- sub3D(free[[el]], t = ifree)
par.not.fixed <- any(((II %x% (II0$R %*% parZ %*% II0$Q)) %*% dpart) != 0)
if (par.not.fixed) {
msg <- c(msg, paste("t=", i, ": For method=kem (EM), if an element of the diagonal of R & Q is 0, the corresponding row of ", par.test, " must be fixed.\n", sep = ""))
ok <- FALSE
}
} # for par.test
} # end for loop over time
############ Check that B^{0} is fixed
# See Summary of Requirements for Degenenate Models in EMDerivation.pdf
Tmax <- 0
for (par.test in c("B")) {
Tmax <- max(Tmax, dim(fixed[[par.test]])[3], dim(free[[par.test]])[3])
}
II0 <- list()
for (i in 1:Tmax) {
for (el in c("Q")) {
ifixed <- min(i, dim(fixed[[el]])[3])
ifree <- min(i, dim(free[[el]])[3])
diag.rows <- 1 + 0:(par.dims[[el]][1] - 1) * (par.dims[[el]][1] + 1)
zero.diags <- is.fixed(free[[el]][diag.rows, , ifree, drop = FALSE], by.row = TRUE) & apply(fixed[[el]][diag.rows, , ifixed, drop = FALSE] == 0, 1, all) # fixed rows
II0[[el]] <- makediag(as.numeric(zero.diags))
}
for (el in c("B")) {
ifree <- min(i, dim(free[[el]])[3])
II <- diag(1, par.dims[[el]][2])
dpart <- sub3D(free[[el]], t = ifree)
par.not.fixed <- any(((II %x% II0$Q) %*% dpart) != 0)
if (par.not.fixed) {
msg <- c(msg, paste("t=", i, ": If an element of the diagonal of Q is 0, the corresponding row and col of ", par.test, " must be fixed.\n", sep = ""))
ok <- FALSE
}
} # for par.test
} # end for loop over time
############ Check that if u^{0} or xi^{0} are estimated, B adjacency matrix is time invariant
# See Summary of Requirements for Degenenate Models in EMDerivation.pdf
Tmax <- 0
for (par.test in c("U")) {
Tmax <- max(Tmax, dim(fixed[[par.test]])[3], dim(free[[par.test]])[3])
}
II0 <- list()
el <- "Q"
ifixed <- min(i, dim(fixed[[el]])[3])
ifree <- min(i, dim(free[[el]])[3])
diag.rows <- 1 + 0:(par.dims[[el]][1] - 1) * (par.dims[[el]][1] + 1)
zero.diags <- is.fixed(free[[el]][diag.rows, , ifree, drop = FALSE], by.row = TRUE) & apply(fixed[[el]][diag.rows, , 1, drop = FALSE] == 0, 1, all) # fixed rows
II0[[el]] <- makediag(as.numeric(zero.diags))
dpart <- sub3D(free$x0, t = 1)
test.adj <- any((II0$Q %*% dpart) != 0)
for (i in 1:Tmax) {
dpart <- sub3D(free$U, t = 1)
test.adj <- test.adj & any((II0$Q %*% dpart) != 0) # II0$Q required to be time constant above
}
if (test.adj) { # means x0^{0} or u^{0} being estimated
Tmax <- 0
for (par.test in c("U", "B")) {
Tmax <- max(Tmax, dim(fixed[[par.test]])[3], dim(free[[par.test]])[3])
}
for (i in 1:Tmax) {
el <- "B"
ifree <- min(i, dim(free[[el]])[3])
# I don't know what par$B will be. I want to create a par$B where there is a non-zero value for any potentially non-zero B's
tmp.MLEobj <- list(marss = MODELobj, par = list(B = matrix(1, dim(free$B)[2], 1)))
tmp.MLEobj$marss$fixed[[el]][tmp.MLEobj$marss$fixed[[el]] != 0] <- 1
tmp.MLEobj$marss$free[[el]][tmp.MLEobj$marss$free[[el]] != 0] <- 1
adjB <- parmat(tmp.MLEobj, el, t = i)[[el]]
adjB[adjB != 0] <- 1
adjB <- unname(adjB)
if (i == 1) {
adjB.1 <- adjB
} else {
if (!isTRUE(all.equal(adjB, adjB.1))) {
msg <- c(msg, paste("t=", i, ": If u^{0} or xi^{0} are estimated, the adjacency matrix specified by B must be time constant.\n", sep = ""))
ok <- FALSE
}
}
} # end for loop over time
} # if the adj test needs to be done
############ Check that II_q^{d} and II_q^{is} are time constant
# See Summary of Requirements for Degenenate Models in EMDerivation.pdf
Tmax <- 0
for (par.test in c("Q")) {
Tmax <- max(Tmax, dim(fixed[[par.test]])[3], dim(free[[par.test]])[3])
}
IIz <- IIp <- OMGz <- OMGp <- IId <- IIis <- list()
for (i in 1:Tmax) {
el <- "Q"
ifixed <- min(i, dim(fixed[[el]])[3])
ifree <- min(i, dim(free[[el]])[3])
diag.rows <- 1 + 0:(par.dims[[el]][1] - 1) * (par.dims[[el]][1] + 1)
zero.diags <- is.fixed(free[[el]][diag.rows, , ifree, drop = FALSE], by.row = TRUE) & apply(fixed[[el]][diag.rows, , ifixed, drop = FALSE] == 0, 1, all) # fixed rows
IIz[[el]] <- makediag(as.numeric(zero.diags))
IIp[[el]] <- diag(1, par.dims[[el]][1]) - IIz[[el]]
OMGz[[el]] <- diag(1, par.dims[[el]][1])[diag(IIz[[el]]) == 1, , drop = FALSE]
OMGp[[el]] <- diag(1, par.dims[[el]][1])[diag(IIp[[el]]) == 1, , drop = FALSE]
# I don't know what par$B will be. I want to create a par$B where there is a non-zero value for any potentially non-zero B's
tmp.MLEobj <- list(marss = MODELobj, par = list(B = matrix(1, dim(free$B)[2], 1)))
tmp.MLEobj$marss$fixed$B[tmp.MLEobj$marss$fixed$B != 0] <- 1
tmp.MLEobj$marss$free$B[tmp.MLEobj$marss$free$B != 0] <- 1
Adj.mat <- parmat(tmp.MLEobj, "B", t = i)$B
Adj.mat[Adj.mat != 0] <- 1
Adj.mat <- unname(Adj.mat)
Adj.mat.pow.m <- matrix.power(Adj.mat, m) # to find all the linkages
Q.0.rows.of.Adj.mat <- OMGz[[el]] %*% Adj.mat.pow.m
# These are the columns corresponding to the directly stochastic bit: Q.0.rows.of.Adj.mat%*%t(OMGp$Q)
if (dim(OMGp[[el]])[1] != 0) {
tmp <- Q.0.rows.of.Adj.mat
# which rows of the + columns are all zero
if (dim(Q.0.rows.of.Adj.mat)[1] != 0) tmp <- apply(Q.0.rows.of.Adj.mat %*% t(OMGp[[el]]) == 0, 1, all)
tmp <- t(OMGz[[el]]) %*% matrix(as.numeric(tmp), ncol = 1) # expand back outl 1s where the deterministic x's are
IId[[el]] <- makediag(tmp)
IId[[el]] <- unname(IId[[el]])
IIis[[el]] <- diag(1, m) - IId[[el]] - IIp[[el]]
IIis[[el]] <- unname(IIis[[el]])
} else { # Q all 0
IId[[el]] <- diag(1, par.dims[[el]][1])
IIis[[el]] <- diag(0, par.dims[[el]][1])
}
if (i == 1) {
IId.1 <- IId[[el]]
IIis.1 <- IIis[[el]]
} else {
if (!isTRUE(all.equal(IId[[el]], IId.1))) {
msg <- c(msg, paste("t=", i, ": The location of the deterministic x's must be time constant.\n", sep = ""))
ok <- FALSE
}
if (!isTRUE(all.equal(IIis[[el]], IIis.1))) {
msg <- c(msg, paste("t=", i, ": The location of the indirectly stochastic x's must be time constant.\n", sep = ""))
ok <- FALSE
}
}
} # end for loop over time
return(list(ok = ok, msg = msg))
}
nwfsc-timeseries/MARSS documentation built on June 3, 2023, 1:32 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions MARSSkemcheck
Documented in MARSSkemcheck
MARSSkemcheck <- function(MLEobj) {
# This checks that the model can be handled by the MARSSkem algorithm
# Most of this is implementing the restrictions in Summary of Requirements for Degenerate Models in derivation
MODELobj <- MLEobj$marss
fixed <- MODELobj$fixed
free <- MODELobj$free
par.dims <- attr(MODELobj, "model.dims")
m <- par.dims[["x"]][1]
n <- par.dims[["y"]][1]
TT <- par.dims[["data"]][2]
pseudolim <- 1E-8
ok <- TRUE
msg <- NULL
# If TT=2 then kf will break
if (TT <= 2) {
msg <- c(msg, "The number of time steps is <=2.\nMore than 2 data points are needed to estimate parameters.\n")
ok <- FALSE
}
############ Check that if fixed, B is within the unit circle
for (t in 1:max(dim(free$B)[3], dim(fixed$B)[3])) {
ifixed <- min(t, dim(fixed$B)[3])
if (is.fixed(free$B[, , ifixed, drop = FALSE])) { # works on 3D matrices
# parmat needs MODELobj and par list
if (is.null(MLEobj$par$B)) tmpparB <- MLEobj$start$B else tmpparB <- MLEobj$par$B
tmp.MLEobj <- list(marss = MODELobj, par = list(B = tmpparB)) # B is fixed but par might have cols from other times
parB <- parmat(tmp.MLEobj, "B", t = t)$B
if (!is.unitcircle(parB)) {
msg <- c(msg, " All the eigenvalues of B must be within the unit circle: all(abs(eigen(B)$values)<=1)\n")
ok <- FALSE
}
}
} # end for over time to check B
############ Check that if R has 0s, then the corresponding row of A and Z are fixed
# See Summary of Requirements for Degenenate Models in EMDerivation.pdf
el <- "R"
# extracts the r=0 rows (or cols)
diag.rows <- 1 + 0:(par.dims[[el]][1] - 1) * (par.dims[[el]][1] + 1)
Tmax <- 0
for (par.test in c("R", "Z", "A")) {
Tmax <- max(Tmax, dim(fixed[[par.test]])[3], dim(free[[par.test]])[3])
}
if (is.null(MLEobj[["par"]])) MLEobj$par <- MLEobj$start
# Run with MARSSkfss because it has more internal checks and error messages
MLEobj$kf <- MARSSkfss(MLEobj)
if (!MLEobj$kf$ok) {
return(list(ok = FALSE, msg = c(MLEobj$kf$errors, msg)))
}
MLEobj$Ey <- MARSShatyt(MLEobj)
msg.tmp <- NULL
for (i in 1:TT) {
ifixed <- min(i, dim(fixed[[el]])[3])
ifree <- min(i, dim(free[[el]])[3])
zero.diags <- is.fixed(free[[el]][diag.rows, , ifree, drop = FALSE], by.row = TRUE) & apply(fixed[[el]][diag.rows, , ifixed, drop = FALSE] == 0, 1, all) # fixed rows and 0
if (any(zero.diags)) {
II0 <- makediag(as.numeric(zero.diags))
if (i <= Tmax) {
for (par.test in c("Z", "A")) {
II <- diag(1, par.dims[[par.test]][2])
ifree.par <- min(i, dim(free[[par.test]])[3])
dpart <- sub3D(free[[par.test]], t = ifree.par)
par.not.fixed <- any(((II %x% II0) %*% dpart) != 0)
if (par.not.fixed) {
if (par.test == "Z") msg <- c(msg, paste("t=", i, ": For method=kem (EM), if an element of the diagonal of R is 0, the corresponding row of Z must be fixed. See MARSSinfo('AZR0').\n", sep = ""))
if (par.test == "A") msg <- c(msg, paste("t=", i, ": For method=kem (EM), if an element of the diagonal of R is 0, the corresponding row of both A and D must be fixed. See MARSSinfo('AZR0').\n", sep = ""))
}
} # for par.test
} # if par needs to be tested; don't need to do over all TT just up to Tmax
Z <- parmat(MLEobj, "Z", t = i)$Z
A <- parmat(MLEobj, "A", t = i)$A
Z.R0 <- Z[zero.diags, , drop = FALSE]
A.R0 <- A[zero.diags, , drop = FALSE]
y.R0 <- MLEobj$Ey$ytT[zero.diags, i, drop = FALSE]
y.resid <- y.R0 - Z.R0 %*% MLEobj$kf$xtT[, i, drop = FALSE] - A.R0
if (any(y.resid > pseudolim)) {
msg.tmp <- c(msg.tmp, paste(" Z, A, and y for the R=0 rows at time=", i, " do not agree. For the R=0 rows, E(y) must equal Z * E(x) + A.\n", sep = ""))
}
} # any zero.diags
} # end for loop over time
if (!is.null(msg)) {
msg <- c(msg, msg.tmp[1:min(10, length(msg.tmp))]) # error msgs could be really long, only print first 10
ok <- FALSE
}
############ Check that II_q^{0} is time constant
# See Summary of Requirements for Degenenate Models in EMDerivation.pdf
Tmax <- 0
for (par.test in c("Q")) {
Tmax <- max(Tmax, dim(fixed[[par.test]])[3], dim(free[[par.test]])[3])
}
II0 <- list()
for (i in 1:Tmax) {
for (el in c("Q")) {
ifixed <- min(i, dim(fixed[[el]])[3])
ifree <- min(i, dim(free[[el]])[3])
diag.rows <- 1 + 0:(par.dims[[el]][1] - 1) * (par.dims[[el]][1] + 1)
zero.diags <- is.fixed(free[[el]][diag.rows, , ifree, drop = FALSE], by.row = TRUE) & apply(fixed[[el]][diag.rows, , ifixed, drop = FALSE] == 0, 1, all) # fixed rows
II0Q <- makediag(as.numeric(zero.diags))
II0Q <- unname(II0[[el]])
}
if (i == 1) {
II0Q.1 <- II0Q
} else {
if (!isTRUE(all.equal(II0Q, II0Q.1))) {
msg <- c(msg, paste("t=", i, ": The placement of 0 variances in Q must be time constant.\n", sep = ""))
ok <- FALSE
}
}
} # end for loop over time
############ Check that if R and Q both have 0s, then the U and Bs are appropriately fixed
# See Summary of Requirements for Degenenate Models in EMDerivation.pdf
Tmax <- 0
for (par.test in c("R", "Q", "U", "B", "Z")) {
Tmax <- max(Tmax, dim(fixed[[par.test]])[3], dim(free[[par.test]])[3])
}
II0 <- list()
for (i in 1:Tmax) {
for (el in c("R", "Q")) {
ifixed <- min(i, dim(fixed[[el]])[3])
ifree <- min(i, dim(free[[el]])[3])
diag.rows <- 1 + 0:(par.dims[[el]][1] - 1) * (par.dims[[el]][1] + 1)
zero.diags <- is.fixed(free[[el]][diag.rows, , ifree, drop = FALSE], by.row = TRUE) & apply(fixed[[el]][diag.rows, , ifixed, drop = FALSE] == 0, 1, all) # fixed rows
II0[[el]] <- makediag(as.numeric(zero.diags))
}
for (el in c("B", "U")) {
ifree <- min(i, dim(free[[el]])[3])
# I don't know what par$Z will be. I want to create a par$Z where there is a non-zero value for any potentially non-zero Z's
tmp.MLEobj <- list(marss = MODELobj, par = list(Z = matrix(1, dim(free$Z)[2], 1)))
tmp.MLEobj$marss$fixed$Z[tmp.MLEobj$marss$fixed$Z != 0] <- 1
tmp.MLEobj$marss$free$Z[tmp.MLEobj$marss$free$Z != 0] <- 1
parZ <- parmat(tmp.MLEobj, "Z", t = i)$Z
II <- diag(1, par.dims[[el]][2])
dpart <- sub3D(free[[el]], t = ifree)
par.not.fixed <- any(((II %x% (II0$R %*% parZ %*% II0$Q)) %*% dpart) != 0)
if (par.not.fixed) {
msg <- c(msg, paste("t=", i, ": For method=kem (EM), if an element of the diagonal of R & Q is 0, the corresponding row of ", par.test, " must be fixed.\n", sep = ""))
ok <- FALSE
}
} # for par.test
} # end for loop over time
############ Check that B^{0} is fixed
# See Summary of Requirements for Degenenate Models in EMDerivation.pdf
Tmax <- 0
for (par.test in c("B")) {
Tmax <- max(Tmax, dim(fixed[[par.test]])[3], dim(free[[par.test]])[3])
}
II0 <- list()
for (i in 1:Tmax) {
for (el in c("Q")) {
ifixed <- min(i, dim(fixed[[el]])[3])
ifree <- min(i, dim(free[[el]])[3])
diag.rows <- 1 + 0:(par.dims[[el]][1] - 1) * (par.dims[[el]][1] + 1)
zero.diags <- is.fixed(free[[el]][diag.rows, , ifree, drop = FALSE], by.row = TRUE) & apply(fixed[[el]][diag.rows, , ifixed, drop = FALSE] == 0, 1, all) # fixed rows
II0[[el]] <- makediag(as.numeric(zero.diags))
}
for (el in c("B")) {
ifree <- min(i, dim(free[[el]])[3])
II <- diag(1, par.dims[[el]][2])
dpart <- sub3D(free[[el]], t = ifree)
par.not.fixed <- any(((II %x% II0$Q) %*% dpart) != 0)
if (par.not.fixed) {
msg <- c(msg, paste("t=", i, ": If an element of the diagonal of Q is 0, the corresponding row and col of ", par.test, " must be fixed.\n", sep = ""))
ok <- FALSE
}
} # for par.test
} # end for loop over time
############ Check that if u^{0} or xi^{0} are estimated, B adjacency matrix is time invariant
# See Summary of Requirements for Degenenate Models in EMDerivation.pdf
Tmax <- 0
for (par.test in c("U")) {
Tmax <- max(Tmax, dim(fixed[[par.test]])[3], dim(free[[par.test]])[3])
}
II0 <- list()
el <- "Q"
ifixed <- min(i, dim(fixed[[el]])[3])
ifree <- min(i, dim(free[[el]])[3])
diag.rows <- 1 + 0:(par.dims[[el]][1] - 1) * (par.dims[[el]][1] + 1)
zero.diags <- is.fixed(free[[el]][diag.rows, , ifree, drop = FALSE], by.row = TRUE) & apply(fixed[[el]][diag.rows, , 1, drop = FALSE] == 0, 1, all) # fixed rows
II0[[el]] <- makediag(as.numeric(zero.diags))
dpart <- sub3D(free$x0, t = 1)
test.adj <- any((II0$Q %*% dpart) != 0)
for (i in 1:Tmax) {
dpart <- sub3D(free$U, t = 1)
test.adj <- test.adj & any((II0$Q %*% dpart) != 0) # II0$Q required to be time constant above
}
if (test.adj) { # means x0^{0} or u^{0} being estimated
Tmax <- 0
for (par.test in c("U", "B")) {
Tmax <- max(Tmax, dim(fixed[[par.test]])[3], dim(free[[par.test]])[3])
}
for (i in 1:Tmax) {
el <- "B"
ifree <- min(i, dim(free[[el]])[3])
# I don't know what par$B will be. I want to create a par$B where there is a non-zero value for any potentially non-zero B's
tmp.MLEobj <- list(marss = MODELobj, par = list(B = matrix(1, dim(free$B)[2], 1)))
tmp.MLEobj$marss$fixed[[el]][tmp.MLEobj$marss$fixed[[el]] != 0] <- 1
tmp.MLEobj$marss$free[[el]][tmp.MLEobj$marss$free[[el]] != 0] <- 1
adjB <- parmat(tmp.MLEobj, el, t = i)[[el]]
adjB[adjB != 0] <- 1
adjB <- unname(adjB)
if (i == 1) {
adjB.1 <- adjB
} else {
if (!isTRUE(all.equal(adjB, adjB.1))) {
msg <- c(msg, paste("t=", i, ": If u^{0} or xi^{0} are estimated, the adjacency matrix specified by B must be time constant.\n", sep = ""))
ok <- FALSE
}
}
} # end for loop over time
} # if the adj test needs to be done
############ Check that II_q^{d} and II_q^{is} are time constant
# See Summary of Requirements for Degenenate Models in EMDerivation.pdf
Tmax <- 0
for (par.test in c("Q")) {
Tmax <- max(Tmax, dim(fixed[[par.test]])[3], dim(free[[par.test]])[3])
}
IIz <- IIp <- OMGz <- OMGp <- IId <- IIis <- list()
for (i in 1:Tmax) {
el <- "Q"
ifixed <- min(i, dim(fixed[[el]])[3])
ifree <- min(i, dim(free[[el]])[3])
diag.rows <- 1 + 0:(par.dims[[el]][1] - 1) * (par.dims[[el]][1] + 1)
zero.diags <- is.fixed(free[[el]][diag.rows, , ifree, drop = FALSE], by.row = TRUE) & apply(fixed[[el]][diag.rows, , ifixed, drop = FALSE] == 0, 1, all) # fixed rows
IIz[[el]] <- makediag(as.numeric(zero.diags))
IIp[[el]] <- diag(1, par.dims[[el]][1]) - IIz[[el]]
OMGz[[el]] <- diag(1, par.dims[[el]][1])[diag(IIz[[el]]) == 1, , drop = FALSE]
OMGp[[el]] <- diag(1, par.dims[[el]][1])[diag(IIp[[el]]) == 1, , drop = FALSE]
# I don't know what par$B will be. I want to create a par$B where there is a non-zero value for any potentially non-zero B's
tmp.MLEobj <- list(marss = MODELobj, par = list(B = matrix(1, dim(free$B)[2], 1)))
tmp.MLEobj$marss$fixed$B[tmp.MLEobj$marss$fixed$B != 0] <- 1
tmp.MLEobj$marss$free$B[tmp.MLEobj$marss$free$B != 0] <- 1
Adj.mat <- parmat(tmp.MLEobj, "B", t = i)$B
Adj.mat[Adj.mat != 0] <- 1
Adj.mat <- unname(Adj.mat)
Adj.mat.pow.m <- matrix.power(Adj.mat, m) # to find all the linkages
Q.0.rows.of.Adj.mat <- OMGz[[el]] %*% Adj.mat.pow.m
# These are the columns corresponding to the directly stochastic bit: Q.0.rows.of.Adj.mat%*%t(OMGp$Q)
if (dim(OMGp[[el]])[1] != 0) {
tmp <- Q.0.rows.of.Adj.mat
# which rows of the + columns are all zero
if (dim(Q.0.rows.of.Adj.mat)[1] != 0) tmp <- apply(Q.0.rows.of.Adj.mat %*% t(OMGp[[el]]) == 0, 1, all)
tmp <- t(OMGz[[el]]) %*% matrix(as.numeric(tmp), ncol = 1) # expand back outl 1s where the deterministic x's are
IId[[el]] <- makediag(tmp)
IId[[el]] <- unname(IId[[el]])
IIis[[el]] <- diag(1, m) - IId[[el]] - IIp[[el]]
IIis[[el]] <- unname(IIis[[el]])
} else { # Q all 0
IId[[el]] <- diag(1, par.dims[[el]][1])
IIis[[el]] <- diag(0, par.dims[[el]][1])
}
if (i == 1) {
IId.1 <- IId[[el]]
IIis.1 <- IIis[[el]]
} else {
if (!isTRUE(all.equal(IId[[el]], IId.1))) {
msg <- c(msg, paste("t=", i, ": The location of the deterministic x's must be time constant.\n", sep = ""))
ok <- FALSE
}
if (!isTRUE(all.equal(IIis[[el]], IIis.1))) {
msg <- c(msg, paste("t=", i, ": The location of the indirectly stochastic x's must be time constant.\n", sep = ""))
ok <- FALSE
}
}
} # end for loop over time
return(list(ok = ok, msg = msg))
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.