Nothing
R/utilities.R
In MultOrd: Multivariate Ordinal Regression Models
## TODO if no intercept in the covGeneral model set default to fix1first
set.threshold.type <- function(rho){
#fixall
if (all(sapply(1:rho$ndim, function(j) all(!is.na(rho$threshold.values[[j]]))))) {#all thresholds are fixed in all dimensions
if (rho$intercept == FALSE) cat("We suggest to include an intercept in the model (formula = y ~ 1 + ...)")
type <- "fixall"
#fix2first
} else if (all(sapply(1:rho$ndim, function(j){
(all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1,2))) && all(which(!is.na(rho$threshold.values[[j]])) == c(1,2)) #all first two thresholds are not NA
) || ((length(rho$threshold.values[[j]]) == 1) && (which(!is.na(rho$threshold.values[[j]])) == 1))
}))){
if (rho$error.structure$type %in% c("corGeneral", "corAR1", "corEqui")) cat("We suggest to fix only one threshold or the intercept in a correlation model.\n")
if ((rho$error.structure$type %in% c("covGeneral"))&& (rho$intercept.type == "fixed")) cat("We suggest to fix either two thresholds or one threshold and the intercept in a covGeneral model.\n")
type <- "fix2first"
#fix2firstlast
} else if (all(sapply(1:rho$ndim, function(j){
(all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1,rho$ntheta[j]))) && all(which(!is.na(rho$threshold.values[[j]])) == c(1,rho$ntheta[j]))#all first and last two thresholds are not NA
) || ((length(rho$threshold.values[[j]]) == 1) && (which(!is.na(rho$threshold.values[[j]])) == 1))
}))){
if (rho$error.structure$type %in% c("corGeneral", "corAR1", "corEqui")) cat("We suggest to fix only one threshold or the intercept in a correlation model.\n")
if ((rho$error.structure$type %in% c("covGeneral"))&& (rho$intercept.type == "fixed")) cat("We suggest to fix either two thresholds or one threshold and the intercept in a covGeneral model.\n")
type <- "fix2firstlast"
#fix1first
} else if (all(sapply(1:rho$ndim, function(j) (length(which(!is.na(rho$threshold.values[[j]])) >= 1) && all(which(!is.na(rho$threshold.values[[j]])) == 1))))){ #all first thresholds are not NA (and no additional)
if ((rho$error.structure$type == "covGeneral") && (rho$intercept.type == "flexible")) stop("Model with covGeneral is not identifiable.
Please either fix two thresholds or one threshold and the intercept.\n", call. = FALSE)
if ((rho$error.structure$type %in% c("corGeneral", "corAR1", "corEqui"))&& (rho$intercept.type == "fixed")) cat("We suggest to fix only one threshold or the intercept in a correlation model.\n")
type <- "fix1first"
#flexible
} else if (all(sapply(1:rho$ndim, function(j) all(is.na(rho$threshold.values[[j]]))))){#all thresholds NA
if (rho$error.structure$type == "covGeneral") stop("Model with covGeneral is not identifiable.
Please either fix two thresholds or one threshold and the intercept.\n", call. = FALSE)
if ((rho$error.structure$type %in% c("corGeneral", "corAR1", "corEqui")) && (rho$intercept.type == "flexible")) stop("Model is not identifiable.
Please either fix one threshold or the intercept.", call. = FALSE)
type <- "flexible"
#ERRORS
} else stop("Either fix all thresholds in one or more outcome dimensions,
or consistently in all other outcome dimensions, all first thresholds or none.\n", call. = FALSE)
if((rho$error.structure$type == "covGeneral") && (rho$binary == TRUE)){
if(type != "fix1first") stop("In the presence of binary outcomes intercept and one threshold
have to be fixed to some value.\n", call. = FALSE)
}
type
}
# #CORRELATION
# } else if (rho$error.structure$type %in% c("corGeneral", "corAR1", "corEqui")){
# #NO INTERCEPT
# if (rho$intercept.type == "fixed"){
# if (all(sapply(1:rho$ndim, function(j) all(is.na(rho$threshold.values[[j]]))))){#all thresholds NA
# type <- "flexible"
# } else stop("If you want fix a threshold value please use intercept in formula (remove ~ 0 + ...)", call. = FALSE)
# #TODO? other parameterization are not allowed for intercept == FALSE?
#
# #INTERCEPT
# } else if (rho$intercept == TRUE){
# if (all(sapply(1:rho$ndim, function(j) all(is.na(rho$threshold.values[[j]]))))){ #all thresholds NA
# if (all(!is.na(rho$coef.values[,1]))){ #intercepts are fixed to constants
# type <- "flexible"
# } else stop("Remove intercept in formula (~ 0 + ...) or fix the first threshold in threshold.values", call. = FALSE)
# }
# if (all(sapply(1:rho$ndim, function(j) (length(which(!is.na(rho$threshold.values[[j]])) >= 1) && all(which(!is.na(rho$threshold.values[[j]])) == 1))))){ #all first thresholds are not NA (and no additional)
# type = "fix1first"
# } else if (any(sapply(1:rho$ndim, function(j) all(!is.na(rho$threshold.values[[j]]))))){#all thresholds are fixed in at least one dimension
# if (all(sapply(1:rho$ndim, function(j) all(!is.na(rho$threshold.values[[j]])) || all(
# (length(which(!is.na(rho$threshold.values[[j]])) >= 1) && (which(!is.na(rho$threshold.values[[j]])) == 1)))))){ #all not NA or only first fixed for each j
# type <- "fix1first"
# } else if (all(sapply(1:rho$ndim, function(j) all(!is.na(rho$threshold.values[[j]])) || all(is.na(rho$threshold.values[[j]]))))){ #all not NA or all NA for each j
# type <- "flexible"
# } else stop("Either fix all thresholds in one or more outcome dimensions,
# or consistently in all other outcome dimensions, all first thresholds or none.
# In addition, intercept must be included in the model.", call. = FALSE)
# } else{
# stop("Model is not identifiable. One of the following could help:
# 1. Remove intercept in formula (~ 0 + ...) or
# 2. Fix all first thresholds in threshold.values or
# 3. Fix all thresholds in threshold.values or
# 4. Either fix all thresholds in one or more outcome dimensions,
# or consistently in all other outcome dimensions, all first thresholds or none", call. = FALSE)
# }
# }
#
# #covGeneral
# } else if (rho$error.structure$type == "covGeneral"){
# if (rho$intercept == FALSE) stop("formula needs intercept if you want to fix a threshold in model CMORcov", call. = FALSE)#TODO?
# if (all(sapply(1:rho$ndim, function(j){
# (all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1,2))) && all(which(!is.na(rho$threshold.values[[j]])) == c(1,2)) #all first two thresholds are not NA
# ) || ((length(rho$threshold.values[[j]]) == 1) && (which(!is.na(rho$threshold.values[[j]])) == 1))
# }))){
# type == "fix2first"
# } else if (all(sapply(1:rho$ndim, function(j){
# (all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1,rho$ntheta[j]))) && all(which(!is.na(rho$threshold.values[[j]])) == c(1,rho$ntheta[j]))#all first and last two thresholds are not NA
# ) || ((length(rho$threshold.values[[j]]) == 1) && (which(!is.na(rho$threshold.values[[j]])) == 1))
# }))){
# type == "fix2firstlast"
# }else if (any(sapply(1:rho$ndim, function(j) all(!is.na(rho$threshold.values[[j]]))))){#all thresholds are fixed in at least one dimensions
# ##CHECK somewhere if all values are fixed in one dimensions + ?warning(sense)
# if (all(sapply(1:rho$ndim, function(j) #
# #(all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1))) && all((which(!is.na(rho$threshold.values[[j]])) == 1))) ||
# (all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1,2))) && all((which(!is.na(rho$threshold.values[[j]])) == c(1,2)))) ||
# (all(!is.na(rho$threshold.values[[j]])))))){
# type <- "fix2first"
# } else if (all(sapply(1:rho$ndim, function(j)
# #(all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1))) && all((which(!is.na(rho$threshold.values[[j]])) == 1))) ||
# (all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1,rho$ntheta[j]))) && all((which(!is.na(rho$threshold.values[[j]])) == c(1,rho$ntheta[j])))) ||
# (all(!is.na(rho$threshold.values[[j]])))))){
# type <- "fix2firstlast"
# } else stop("Either fix all thresholds in one or more outcome dimensions,
# or consistently in all other outcome dimensions, all first two thresholds or all first and last thresholds", call. = FALSE)
# } else stop("Either fix first two or first and last threshold in threshold.values if you use a model with covGeneral", call. = FALSE)
# }
# type
# }
#data.x <- data.multord$x
#ndim <- rho$ndim
set.error.structure <- function(error.structure, data.x, ndim){
if (error.structure$formula == ~1 ) {
if (error.structure$type %in% c("corEqui", "corAR1")) {
error.structure$x <- model.matrix(error.structure$formula, data.x[[1]])
} else {
error.structure$x <- as.factor(model.matrix(error.structure$formula, data.x[[1]]))
}
} else if (error.structure$type %in% c("corEqui", "corAR1")){
tmp <- lapply(1:ndim, function(j) {
tmp1 <- model.matrix(error.structure$formula, data.x[[j]])
attribute <- attr(tmp1, "assign")
tmp2 <- tmp1[match(rownames(data.x[[j]]),rownames(tmp1)),]
rownames(tmp2) <- rownames(data.x[[j]])
attr(tmp2, "assign") <- attribute
tmp2
}
)
error.x <- sapply(1:ncol(tmp[[1]]), function(k){
xtcol <- do.call(cbind,lapply(tmp, `[`,,k))
xtcol_final <- apply(xtcol,1,function(i) unique(i[!is.na(i)]))
if(is.list(xtcol_final)) stop("Covariate dependent variables need to be constant across multiple observations", call.=FALSE)
xtcol_final
})
colnames(error.x) <- colnames(tmp[[1]])
attr(error.x, "assign") <- attr(tmp[[1]], "assign")
error.structure$x <- error.x
# } else if (error.structure$type == "corAR1"){
# #depends on one factor (factor of first year)
# if (!is.factor(data.x[[1]][,all.vars(error.structure$formula)])) stop("variable must be of type factor in corAR1", call. = FALSE)
# tmp <- lapply(1:ndim,function(j) data.x[[j]][,all.vars(error.structure$formula)])
# if (!all(sapply(1:ndim, function(j) all(tmp[[1]]==tmp[[j]], na.rm=T)))) stop("Covariate dependent variables need to be constant across multiple observations", call. = FALSE)
# error.structure$x <- tmp[[1]]#data.x[[1]][,all.vars(error.structure$formula)]
} else if (error.structure$type %in% c("corGeneral", "covGeneral")){
#depends on one factor
if (length(all.vars(error.structure$formula[[2]])) > 1) stop("only one factor is allowed in covGeneral and corGeneral", call. = FALSE)
if (!is.factor(data.x[[1]][,all.vars(error.structure$formula)])) stop("variable must be of type factor in covGeneral and corGeneral", call. = FALSE)
tmp <- lapply(1:ndim,function(j) data.x[[j]][,all.vars(error.structure$formula)])
if (!all(sapply(1:ndim, function(j) all(tmp[[1]]==tmp[[j]], na.rm=T)))) stop("Covariate dependent variables need to be constant across multiple observations", call. = FALSE)
error.structure$x <- as.factor(apply(do.call("cbind.data.frame", tmp), 1, function(x) unique(x[!is.na(x)])))
} else {
#depends on one factor
if (length(all.vars(error.structure$formula[[2]])) > 1) stop("only one factor is allowed in covGeneral and corGeneral", call. = FALSE)
if (!is.factor(data[,all.vars(error.structure$formula)])) stop("variable must be of type factor in covGeneral and corGeneral", call. = FALSE)
error.structure$x <- data.x[,all.vars(error.structure$formula)]
}
error.structure
}
set.error.structure.CMOR <- function(error.structure, data){
if (error.structure$type == "corEqui"){
error.structure$x <- model.matrix(error.structure$formula, data)
} else if (error.structure$type == "corAR1"){
stop("corAR1 is not applicable in CMOR", call. = FALSE)
# if (error.structure$formula == ~1){
# error.structure$x <- as.factor(model.matrix(error.structure$formula, data[[1]]))
# } else{ #depends on one factor (factor of first year)
# if (!is.factor(data[[1]][,all.vars(error.structure$formula)])) stop("variable must be of type factor in covGeneral and corGeneral", call. = FALSE)
# error.structure$x <- data[[1]][,all.vars(error.structure$formula)]
# }
} else if (error.structure$type %in% c("corGeneral", "covGeneral")){
if (error.structure$formula == ~ 1){
error.structure$x <- as.factor(model.matrix(error.structure$formula, data))
} else{ #depends on one factor
if (length(all.vars(error.structure$formula[[2]])) > 1) stop("only one factor is allowed in covGeneral and corGeneral", call. = FALSE)
if (!is.factor(data[,all.vars(error.structure$formula)])) stop("variable must be of type factor in covGeneral and corGeneral", call. = FALSE)
error.structure$x <- data[,all.vars(error.structure$formula)]
}
}
error.structure
}
checkArgs <- function(rho){
#CHECK if treshold.values is in line with threshold.constraints
if (length(rho$threshold.constraints) != rho$ndim) stop("dimensions of threshold.values and number of thresholds do not match", call. = FALSE)
if (any(sapply(1:rho$ndim, function(j) length(rho$threshold.values[[j]]) != rho$ntheta[j])))
stop("dimensions of threshold.values and number of thresholds do not match", call. = FALSE)
for (j in unique(rho$threshold.constraints)){
ind <- which(rho$threshold.constraints == j)
if (length(unique(rho$threshold.values[ind]))!=1) stop("If constraints are set on thresholds (by threshold.constraints), threshold.values need to be specified accordingly for these outcome dimensions.
Maybe dimension do not have the same number of threshold parameters.", call. = FALSE)
}
if (nrow(rho$coef.constraints) != rho$ndim) stop("row dimension of coef.constraints and outcome dimension do not match (?factor)", call. = FALSE)
for (j in 1:ncol(rho$coef.constraints)){
indj <- unique(rho$coef.constraints[,j])
indj <- indj[!is.na(indj)]
lapply(seq_len(length(indj)), function(k) {
tmpind <- which(rho$coef.constraints[,j] == indj[k])
tmp <- rho$coef.values[tmpind,j]
if(length(unique(tmp)) != 1) stop("If constraints are set on the coefficients (by coef.constraints),
coef.values need to be specified accordingly for these outcome dimensions.", call. = FALSE)
})
}
}
##########################################
###### AUXILIARY FUNCTIONS ##############
###########################################################################
# rectbiv.norm.prob_TEST <- function(U, L, r) {
# # computes the rectangle probabilities for biv.normal-distribution
# if (is.null(dim(U)[2]) & is.null(dim(L)[2])) {
# dim(U) <- dim(L) <- c(1, length(U))
# }
# n <- nrow(U)
# x <- c(U[, 1],L[, 1],U[, 1],L[, 1])
# y <- c(U[, 2],U[, 2],L[, 2],L[, 2])
# r
# p <- pbivnorm(x,y,r)
# p[is.nan(p)] <- 0
# p[seq_len(n)]-p[n + seq_len(n)]-p[2*n + seq_len(n)] + p[3*n + seq_len(n)]
# }
rectbiv.norm.prob <- function(U, L, r) {
# computes the rectangle probabilities for biv.normal-distribution
if (is.null(dim(U)[2]) & is.null(dim(L)[2])) {
dim(U) <- dim(L) <- c(1, length(U))
}
p1 <- pbivnorm(U[, 1], U[, 2], r)
p2 <- pbivnorm(L[, 1], U[, 2], r)
p3 <- pbivnorm(U[, 1], L[, 2], r)
p4 <- pbivnorm(L[, 1], L[, 2], r)
p1[is.nan(p1)] <- 0
p2[is.nan(p2)] <- 0
p3[is.nan(p3)] <- 0
p4[is.nan(p4)] <- 0
p1-p2-p3+p4
}
biv.nt.prob2 <-function (df, lower, upper,
mean=c(0,0), r) {
if (df == Inf) {
nu <- 0
} else {
nu <- df
}
rho <- r
if (any(lower > upper)) stop("lower > upper integration limits")
if (any(lower == upper)) return(0)
infin <- as.integer(c(2, 2))
prob <- as.double(0)
a <- .Fortran("smvbvt", prob, nu, lower, upper, infin, rho,
#PACKAGE = "mnormt")
PACKAGE = "MultOrd")
return(a[[1]])
}
casesNA<- function(y) {
# returns a column vector v, assigning a different number/label to each unique combination of NA pattern
contrast <- as.matrix((!is.na(y)) + 0)
contrast.char <- apply(contrast, 1, paste, collapse = "") # make each contrast row as character i.e., 111 has all three response variables
factor(contrast.char, levels = unique(contrast.char))
}
split.NA.pattern <- function(y) {
# y should be a data.frame
# return a list, whose elements are groups of the matrix y according to each unique combination of response variables
if (!is.data.frame(y)) y <- as.data.frame(y)
yg <- split(y, casesNA(y)) # grouped y by cases of NA pattern
lapply(yg, function(x) which(rownames(y) %in% rownames(x)))
}
getStart.values <- function(rho){
# # set starting value for beta and the gammas thresholds
# #thetas <- qlogis((1:rho$ntheta) / (rho$ntheta + 1))
# #if (rho$link != "logistic")
# # thetas <- thetas/1.7
# mclm <- lapply(1:rho$ndim, function(j)
# ordinal:::clm(factor(rho$y[,j]) ~ rho$x[,which(rho$coef.constraints[j,] != 0)], link = "logit"))
# #thetas <- lapply(mclm, "[[", "alpha")
# betas <- lapply(mclm, "[[", "beta") #Reduce("+", lapply(mclm, "[[", "beta"))
# #gammas <- lapply(thetas, function(x) c(x[1L], log(diff(x))))
gammas <- sapply(1:rho$ndim, function(j) {
if (rho$npar.theta.opt[j] != 0){
theta <- if (rho$ntheta[j] >= 2) polr(rho$y[, j] ~1)$zeta else 0#qlogis((1:rho$ntheta[j])/(rho$ntheta[j] + 1))
if (!grepl("logit", rho$link)) theta <- theta/1.7
c(theta[1L], log(diff(theta)))[1:rho$npar.theta.opt[j]]
} else NULL
})
c(unlist(gammas))#, unlist(betas))# rep(, sum(rho$npar.betas))) }
}
transf.par.cor <- function(par, rho) {
sigmas <- rho$transf.sigmas(par[rho$npar.thetas + rho$npar.betas + seq_len(rho$npar.sigmas)],rho)
theta <- rho$transf.thresholds(par[seq_len(rho$npar.thetas)], rho)
par_beta <- par[rho$npar.thetas + seq(rho$npar.betas)]
beta <- sapply(1:ncol(rho$coef.constraints), function(j){
sapply(1:nrow(rho$coef.constraints), function(i,j) ifelse(is.na(rho$ind.coef[i,j]), rho$coef.values[i,j], par_beta[rho$ind.coef[i, j]]), j)
})
pred.fixed <- sapply(1:rho$ndim, function(j) rho$x[[j]] %*% beta[j, ])
theta.lower <- sapply(1:rho$ndim, function(j) c(-rho$inf.value, theta[[j]])[rho$y[, j]])
theta.upper <- sapply(1:rho$ndim, function(j) c(theta[[j]], rho$inf.value)[rho$y[, j]])
pred.lower <- (theta.lower - pred.fixed)/rho$sd.y
pred.upper <- (theta.upper - pred.fixed)/rho$sd.y
list(U = pred.upper,
L = pred.lower,
sigmas = sigmas)
}
transf.par.cov <- function(par, rho) {
exp.par.sd <- exp(par[rho$npar.thetas + rho$npar.betas +
rho$npar.cor * rho$ncor.levels +
seq_len(rho$npar.cor.sd * rho$ncor.levels)])
sigmas <- rho$transf.sigmas(par[rho$npar.thetas + rho$npar.betas + seq_len(rho$npar.sigmas)],rho,
exp.par.sd)
exp.par.sd <- lapply(1:rho$ncor.levels, function(l) exp.par.sd[ (l-1) * rho$npar.cor.sd + seq_len(rho$npar.cor.sd) ])
theta <- rho$transf.thresholds(par[seq_len(rho$npar.thetas)], rho)
par_beta <- par[rho$npar.thetas + seq_len(rho$npar.betas)]
beta <- sapply(1:ncol(rho$coef.constraints), function(j){
sapply(1:nrow(rho$coef.constraints), function(i,j)
ifelse(is.na(rho$ind.coef[i,j]), rho$coef.values[i,j], par_beta[rho$ind.coef[i, j]]), j)
})
pred.fixed <- sapply(1:rho$ndim, function(j) rho$x[[j]] %*% beta[j, ])
theta.lower <- sapply(1:rho$ndim, function(j) c(-rho$inf.value, theta[[j]])[rho$y[, j]])
theta.upper <- sapply(1:rho$ndim, function(j) c(theta[[j]], rho$inf.value)[rho$y[, j]])
## make matrix of exp.par.sd
lev <- match(rho$error.structure$x, rho$error.structure$levels)
#pred.lower <- (theta.lower - pred.fixed)/rho$sd.y
#pred.upper <- (theta.lower - pred.fixed)/rho$sd.y
#exp.par.sd.mat <- matrix(ncol = ncol(theta.upper), nrow = nrow(theta.upper))
#for (l in 1:rho$ncor.levels){
# exp.par.sd.mat[lev == l, ] <-
# exp.par.sd[(l - 1) * rho$npar.cor.sd + seq_len(rho$npar.cor.sd)]
#}
#for (i in 1:nrow(pred.lower)){
# pred.lower[i,] <- pred.lower[i,]/exp.par.sd[(lev[i] - 1) * rho$npar.cor.sd + seq_len(rho$npar.cor.sd)]
# pred.upper[i,] <- pred.upper[i,]/exp.par.sd[(lev[i] - 1) * rho$npar.cor.sd + seq_len(rho$npar.cor.sd)]
#}
pred.lower <- theta.lower - pred.fixed
pred.upper <- theta.upper - pred.fixed
for (l in 1:rho$ncor.levels){
pred.lower[lev==l,] <-t(t(pred.lower[lev==l,]/rho$sd.y)/exp.par.sd[[l]])
pred.upper[lev==l,] <-t(t((pred.upper[lev==l,])/rho$sd.y)/exp.par.sd[[l]])
}
#pred.lower <-t(t((theta.lower - pred.fixed)/rho$sd.y)/exp.par.sd)
#pred.upper <-t(t((theta.upper - pred.fixed)/rho$sd.y)/exp.par.sd)
# t(t(pred.upper)/exp.par.sd)
# test speed up
sigmas <- lapply(sigmas, cov2cor)#cov2cor(sigmas)
list(U = pred.upper,
L = pred.lower,
sigmas = sigmas,
std.dev = exp.par.sd)
}
#########################################################################
## transformation of the threshold parameters (to ensure monotonicity) ##
#########################################################################
transf.thresholds.fixall <- function(gamma, rho){
rho$threshold.values
}
transf.thresholds.flexible <- function(gamma, rho){
lapply(1:rho$ndim, function(j)
if (anyNA(rho$threshold.values[[j]])){
if (rho$ntheta[j] > 1) {
cumsum(c(gamma[rho$ind.thresholds[[j]][1]],
exp(gamma[rho$ind.thresholds[[j]][2:rho$ntheta[j]]])))
} else if (rho$ntheta[j] == 1) gamma[rho$ind.thresholds[[j]]] else NULL
} else rho$threshold.values[[j]]
)
}
transf.thresholds.fix1.first <- function(gamma, rho){
lapply(1:rho$ndim, function(j)
if (length(rho$ind.thresholds[[j]]) >= 1) {
cumsum(c(rho$threshold.values.fixed[[j]][1],
exp(gamma[rho$ind.thresholds[[j]][seq_len(length(rho$ind.thresholds[[j]]))]])))
} else if (length(rho$ind.thresholds[[j]]) == 0) rho$threshold.values.fixed[[j]] else NULL)
}
#gamma <- par[seq_len(rho$npar.thetas)]
transf.thresholds.fix2.first <- function(gamma, rho){
lapply(1:rho$ndim, function(j)
if (rho$npar.theta[j] >= 1) {
c(rho$threshold.values.fixed[[j]][1],
cumsum(c(rho$threshold.values.fixed[[j]][2] ,exp(gamma[rho$first.ind.theta[j] - 1 +
seq_len(rho$npar.theta[j])]))))
} else {
if (rho$npar.theta[j] == 0) rho$threshold.values.fixed[[j]]
}
)
}
# transf.thresholds.fix2.first <- function(gamma, rho){
# lapply(1:rho$ndim, function(j)
# if (rho$npar.theta[rho$threshold.constraints[j]] >= 1) {
# c(rho$threshold.values.fixed[[j]][1],
# cumsum(c(rho$threshold.values.fixed[[j]][2] ,exp(gamma[rho$first.ind.theta[rho$threshold.constraints[j]] - 1 +
# seq_len(rho$npar.theta[rho$threshold.constraints[j]])]))))
# } else {
# if (rho$npar.theta[rho$threshold.constraints[j]] == 0) rho$threshold.values.fixed[[j]]
# }
# )
# }
transf.thresholds.fix2.firstlast <- function(gamma, rho){
lapply(1:rho$ndim, function(j){
gamma1 <- gamma[rho$first.ind.theta[j] - 1 + seq_len(rho$npar.theta[j])]
recursive.theta <- function(i) {
if (i == 0) 0
else return ((exp(gamma1[i]) +
recursive.theta(i - 1))/(1 + exp(gamma1[i])))
}
if (rho$npar.theta[j] >= 1) {
theta <- sapply(1:rho$npar.theta[j], function(i)
recursive.theta(i))
c(0, theta, 1) * (rho$threshold.values.fixed[[j]][2] - rho$threshold.values.fixed[[j]][1]) +
rho$threshold.values.fixed[[j]][1]
} else {
rho$threshold.values.fixed[[j]]
}
})
}
# transf.thresholds.fix2.firstlast <- function(gamma, rho){
# lapply(1:rho$ndim, function(j){
# gamma1 <- gamma[rho$first.ind.theta[j] - 1 + seq_len(rho$npar.theta[rho$threshold.constraints[j]])]
# recursive.theta <- function(i) {
# if (i == 0) 0
# else return ((exp(gamma1[i]) +
# recursive.theta(i - 1))/(1 + exp(gamma1[i])))
# }
# if (rho$npar.theta[rho$threshold.constraints[j]] >= 1) {
# theta <- sapply(1:rho$npar.theta[rho$threshold.constraints[j]], function(i)
# recursive.theta(i))
# c(0, theta, 1) * (rho$threshold.values.fixed[[j]][2] - rho$threshold.values.fixed[[j]][1]) +
# rho$threshold.values.fixed[[j]][1]
# } else {
# rho$threshold.values.fixed[[j]]
# }
# })
# }
##############################################################################
## transformation of the parameters of the correlation/covariance structure ##
##############################################################################
transf.sigmas.spheric <- function(par.nu, rho, exp.par.sd = NULL) {
lapply(1:rho$ncor.levels, function(l) {
nu <- par.nu[(l - 1) * rho$npar.cor + seq_len(rho$npar.cor)]
angles <- pi * exp(nu)/(1 + exp(nu))
cosmat <- diag(rho$ndim)
cosmat[lower.tri(cosmat)] <- cos(angles)
S1 <- matrix(0, nrow = rho$ndim, ncol = rho$ndim)
S1[, 1L] <- 1
S1[-1L, -1L][lower.tri(S1[-1L, -1L], diag = T)] <- sin(angles)
tLmat <- sapply(1:rho$ndim,
function(j) cosmat[j, ] * cumprod(S1[j, ]))
sigma <- crossprod(tLmat)
if (is.null(exp.par.sd)) {
return(sigma)
} else {
stdev <- exp.par.sd[(l - 1) * rho$npar.cor.sd + seq_len(rho$npar.cor.sd)]
return(t(stdev * sigma) * stdev)
}
})
}
# transf.sigmas.corAR1 <- function(par.sigma, rho, exp.par.sd = NULL) {
# lapply(1:rho$ncor.levels, function(l){
# ###TODO include covariate dependent corAR1? (like corEqui)
# r <- z2r(par.sigma[seq_len(rho$npar.cor) + l -1])
# #r <- 0.9
# sigma <- diag(rho$ndim)
# sigma[lower.tri(sigma)] <- r^sequence((rho$ndim-1):1)
# ### CHECK: easier way?
# sigma <- sigma + t(sigma) - diag(diag(sigma))
# sigma
# })
# }
transf.sigmas.corAR1 <- function(par.sigma, rho, exp.par.sd = NULL) {
w <- par.sigma[seq_len(rho$npar.cor)]
z <- rho$error.structure$x %*% w
r <- z2r(z)
lapply(seq_len(length(r)), function(i){
sigma <- diag(rho$ndim)
sigma[lower.tri(sigma)] <- r[i]^sequence((rho$ndim-1):1)
### CHECK: easier way?
sigma <- sigma + t(sigma) - diag(diag(sigma))
sigma
})
}
transf.sigmas.corEqui <- function(par.sigma, rho, exp.par.sd = NULL) {
w <- par.sigma[seq_len(rho$npar.cor)]
z <- rho$error.structure$x %*% w
z2r(z)
}
z2r <- function (z) {
ifelse(z > 354, 1, (exp(2 * z) - 1)/(1 + exp(2 * z)))
}
getInd.coef <- function(coef.constraints, coef.values) {
ind <- matrix(NA,ncol = ncol(coef.constraints), nrow = nrow(coef.constraints))
k <- 1
i <- 1
for(j in 1:ncol(coef.constraints)){
if (is.na(coef.values[i,j])){
ind[i,j] <- k
k <- k + 1
}
}
for(i in 2:nrow(coef.constraints)){
for(j in 1:ncol(coef.constraints)){
if (is.na(coef.values[i,j]) && !(coef.constraints[i,j] %in% coef.constraints[1:(i-1),j])){
ind[i,j] <- k
k <- k + 1
}
if (is.na(coef.values[i,j]) && coef.constraints[i,j] %in% coef.constraints[1:(i-1),j]){
ind[i,j] <- ind[min(which(coef.constraints[i,j] == coef.constraints[1:(i-1),j])),j]
}
}
}
ind
}
getInd.thresholds.flexible <- function(threshold.constraints,rho){
pointer <- lapply(1:rho$ndim, function(j) rep(0,rho$npar.theta[j]))
pointer[[1]] <- seq_len(rho$npar.theta[1])
k <- length(pointer[[1]])
for(i in 2:rho$ndim){
if (rho$npar.theta[i] == 0){
pointer[[i]] <- integer(0)
} else {
if (threshold.constraints[i] %in% threshold.constraints[1:(i-1)]){
min.ind <- min(which(threshold.constraints[i] == threshold.constraints[1:(i-1)]))
if (rho$ntheta[i] != rho$ntheta[min.ind])
stop("Constraints on threshold parameters are not valid (different number of categories)", call. = FALSE)
pointer[[i]] <- pointer[[min.ind]]
} else{
pointer[[i]] <- (k+1):(k+rho$ntheta[i])
k <- k + rho$ntheta[i]
}
}
}
pointer
}
getInd.thresholds.fix2 <- function(threshold.constraints,rho) {
#pointer <- lapply((1:rho$ndim)[which(rho$npar.theta > 0)], function(j) 1:rho$npar.theta[j])
pointer <- lapply((1:rho$ndim), function(j) seq_len(rho$npar.theta[j]))
k <- length(pointer[[1]])
#for(i in (1:rho$ndim)[which(rho$npar.theta > 0)][-1]){
for(i in (1:rho$ndim)){
if (rho$npar.theta[i] == 0){
pointer[[i]] <- integer(0)
} else {
if (threshold.constraints[i] %in% threshold.constraints[1:(i - 1)]){
min.ind <- min(which(threshold.constraints[i] == threshold.constraints[1:(i-1)]))
if (rho$ntheta[i] != rho$ntheta[min.ind])
stop("Constraints on threshold parameters are not valid (different number of categories)", call. = FALSE)
pointer[[i]] <- pointer[[min.ind]]
} else{
pointer[[i]] <- (k+1):(k+rho$ntheta[i]-2)
k <- k + rho$ntheta[i]-2
}
}
}
pointer
}
getInd.thresholds.fix1 <- function(threshold.constraints,rho) {
pointer <- lapply(1:rho$ndim, function(j) rep(0,rho$npar.theta[j]))
pointer[[1]] <- seq_len(rho$npar.theta[1])#1:(rho$ntheta[1]-1)
k <- length(pointer[[1]])
for(i in 2:rho$ndim){
if (rho$npar.theta[i] == 0){
pointer[[i]] <- integer(0)
} else {
if (threshold.constraints[i] %in% threshold.constraints[1:(i - 1)]){
min.ind <- min(which(threshold.constraints[i] == threshold.constraints[1:(i - 1)]))
if (rho$ntheta[i] != rho$ntheta[min.ind])
stop("Constraints on threshold parameters are not valid (different number of categories)", call. = FALSE)
pointer[[i]] <- pointer[[min.ind]]
} else{
pointer[[i]] <- (k+1):(k+rho$ntheta[i]-1)
k <- k + rho$ntheta[i]-1
}
}
}
pointer
}
getInd.thresholds.fixall <- function(threshold.constraints,rho) {
lapply(1:rho$ndim, function(j) integer())
}
get.labels.theta <- function(rho,j) {
lev <- levels(rho$y[, j])
sapply(1:(rho$ntheta[j]), function(i){
paste(lev[i], lev[i + 1], sep = "|")
})
}
#' @title Error Structures in Multord
#' @description Different \code{error.structures} are available in \pkg{MultOrd}:
#' \itemize{
#' \item general correlation structure (default) \code{corGeneral(~1)},
#' \item general covariance structure \code{covGeneral(~1)},
#' \item factor dependent correlation structure \code{covGeneral(~f)},
#' \item factor dependent covariance structure \code{covGeneral(~f)},
#' \item covariate dependent equicorrelation structure \code{corEqui(~X)},
#' \item AR(1) correlation structure \code{corAR1(~1)}, or
#' \item factor dependent AR(1) correlation structure \code{corAR1(~f)}.
#' }
#' See \code{\link{error.structures}} or 'Details'.
#' @details General (symmetric) correlation structure (cross-sectional)
#' @param formula \code{\link{formula}} object
#' @export
#' @name error.structures
corGeneral <- function(formula) UseMethod("corGeneral")
#' @rdname error.structures
#' @export
corGeneral.formula <- function(formula){
return(list(type = "corGeneral", formula = formula))
}
#' @details Equicorelation structure (cross-sectional)
# #' @param formula \code{\link{formula}} object
#' @export
#' @rdname error.structures
corEqui <- function(formula) UseMethod("corEqui")
#' @rdname error.structures
#' @export
corEqui.formula <- function(formula){
return(list(type = "corEqui", formula = formula))
}
#' @details General (symmetric) covariance structure (cross-sectional)
# #' @param formula \code{\link{formula}} object
#' @export
#' @rdname error.structures
covGeneral <- function(formula) UseMethod("covGeneral")
#' @rdname error.structures
#' @export
covGeneral.formula <- function(formula){
return(list(type = "covGeneral", formula = formula))
}
#' @details AR(1) correlation structure (panel; only applicable in line with \code{\link{multord}})
# #' @param formula \code{\link{formula}} object
#' @export
#' @rdname error.structures
corAR1 <- function(formula) UseMethod("corAR1")
#' @rdname error.structures
#' @export
corAR1.formula <- function(formula){
return(list(type = "corAR1", formula = formula))
}
backtransf.sigmas <- function(R){
J <- nrow(R)
l <- t(chol(R))
angmat <- matrix(1,ncol=J,nrow=J)
angmat[-1,1] <- acos(l[-1,1])
for (j in 2:(J-1)){
sinprod <- apply(sin(angmat[, seq_len(j-1), drop=F]), 1, prod) ## denominator in division
angmat[-(1:j),j]<-acos((l/sinprod)[-(1:j),j])
}
angdivpi <- angmat[lower.tri(angmat)]/pi
log(angdivpi/(1-angdivpi))
}
get.sigma.i <- function(object) {
if (object$rho$error.structure$type == "corEqui") {
lapply(1:object$rho$n, function(i) {
tmp <- matrix(object$r[i], nrow = object$rho$ndim, ncol = object$rho$ndim)
diag(tmp) <- 1
tmp
})
} else {
lev <- match(object$rho$error.structure$x, object$rho$error.structure$levels)
lapply(1:object$rho$n, function(i) object$sigmas[lev[i]][[1]])
}
}
#' @title Data preparation for multord
#'
#' @description
#' This function is an (internally) used to transforms the \code{data}, into a "multivariate setting",
#' where all repeated measurements are matched accordingly to their ID. A matrix of all ordinal responses with \code{J} columns
#' as well as a list of length \code{J} of matrices with all the covariates are created.
#' @param data a \code{data.frame}, where each row corresponds to a single measurement.
#' @param index is an (optional) argument that specifies the index for the subjects and the response index of the multiple measurement.
#' This is usually performed
#' by a character vector of length two specifying the column names of the subject index and
#' the multiple response index in data. The default value of index is NULL assuming that the
#' first column of data contains the subject index and the second column the multiple response index.
#' @param y.names column name of \code{data} where the ordinal observations are stored.
#' @param x.names column names of all the covariates in {data}.
#' @param y.levels (optional) list of length \code{J} that specifies the levels of each repeated measurement. If the categories
#' differ across repeated measurements (either the number of categories or the category labels) it is recommended to set them.
#' @param response.names (optional) vector of names of the repeated measurements in \code{data}
#'which specifies the ordering of the repeated measurements.
#' @export
#y.names <- "yy"
#y.names <- rho$response.name
#x.names <- rho$x.names
# data <- data_multord
# y.names <- rho$response.name
# x.names <- unique(c(rho$x.names,weights))
# y.levels = response.levels
# response.names = response.names = c("rater1", "rater2", "rater3")
# index = c("firmID", "rater")
multord.data <- function(data, index, y.names, x.names,
y.levels = NULL, response.names = NULL) {
df <- list()
if (any(duplicated(data[,index]))) stop("duplicated observation(s) for one index", call. = FALSE)
index.levels <- levels(as.factor(data[, index[2]]))
data.split <- split(data[,c(y.names, index[1])], data[, index[2]])
data.split.y <- lapply(seq_len(length(data.split)), function(j) {
colnames(data.split[[j]]) <- c(index.levels[j],index[1])
data.split[[j]]})
df$y <- Reduce(function(...) merge(..., by = index[1], all = TRUE), data.split.y )[, -1]
response.names.NA <- c()
if (!is.null(response.names)) {
response.names.NA <- response.names[!response.names %in% colnames(df$y)]
df$y <- cbind(df$y, matrix(NA, ncol = length(response.names.NA), nrow = nrow(df$y), dimnames = list(c(),response.names.NA)))
df$y <- df$y[,as.character(response.names)]
}
colnames.y <- colnames(df$y)
if (is.null(y.levels)) {
df$y <- do.call(cbind.data.frame, lapply(1:ncol(df$y), function(j) ordered(df$y[, j])))
} else {
df$y <- do.call(cbind.data.frame, lapply(1:ncol(df$y), function(j) {
if (!all(unique(df$y[!is.na(df$y[, j]), j]) %in% y.levels[[j]])) stop("levels of response do not match with y.levels", call. = FALSE)
ordered(df$y[, j], levels = y.levels[[j]])}
))
}
colnames(df$y) <- colnames.y#paste0(y.names, 1:ncol(df$y))
#X
data.split.x <- split(data[, c(x.names, index[1])], data[, index[2]])
names.x <- names(data.split.x)
#REMOVE this
#names.x <- names(data.split.x)[names(data.split.x) %in% response.names]
#set colnames (otherwise warning due to identical colnames in reduce)
data.split.x <- lapply(seq_len(length(data.split.x)), function(j) {
colnames(data.split.x[[j]]) <- c(paste(x.names,j, sep = "."), index[1])
data.split.x[[j]]})
xdatadf <- Reduce(function(...) merge(...,by = index[1], all = TRUE), data.split.x)[, -1]
xdatadf <- cbind(xdatadf, matrix(NA, ncol = length(response.names.NA) * length(x.names),
nrow = nrow(xdatadf)))
#xdatadf <- Reduce(function(x,y) merge(x,y,by = index[1], all = TRUE), split(data[,c(x.names, index[1])], data[, index[2]]) )[, -1]
#REMOVE this
# df$x <- lapply(1:ncol(df$y), function(i) {
# tmp <- xdatadf[,(i - 1) * length(x.names) + seq_len(length(x.names))]
# names(tmp) <- x.names
# tmp
# })
df$x <- lapply(1:(length(index.levels) + length(response.names.NA)), function(i) {
tmp <- xdatadf[,(i - 1) * length(x.names) + seq_len(length(x.names))]
names(tmp) <- x.names
tmp
})
names(df$x) <- c(names.x, response.names.NA)
if (!is.null(response.names)) df$x <- df$x[as.character(response.names)]
df
}
check <- function(...){
stopifnot(...)
}
#threshold.values <- list(c(NA,NA,NA),c(2,3,NA),c(1,2,NA,NA,NA))
# set.model <- function(rho){
# if (rho$error.structure$type == "corAR1"){
# m <- "PMOR"
# } else{
# if (rho$error.structure$type == "corGeneral"){
# #if (length(all.vars(rho$error.structure$formula[[2]])) > 1) stop("only one factor is allowed in corGeneral", call. = FALSE)
# rho$model <- "CMORcor"
# }
# if (rho$error.structure$type == "corEqui"){
# rho$model <- "CMORcor"
# }
# if (rho$error.structure$type == "covGeneral"){
# #if (length(all.vars(rho$error.structure$formula[[2]])) > 1) stop("only one factor is allowed in covGeneral", call. = FALSE)
# rho$model <- "CMORcov"
# }
# if (is.null(rho$threshold.values)){
# if ((rho$model == "CMORcor") && (attr(terms.formula(rho$formula), "intercept") == 1)){
# stop("Remove intercept in formula (~ 0 + ...) or fix the first threshold in threshold.values", call. = FALSE)
# }
# if ((rho$model == "CMORcov") && (attr(terms.formula(rho$formula), "intercept") == 0)){
# stop("Include intercept in formula (~ 1 + ...)", call. = FALSE)
# }
# m <- rho$model
# } else{
# if (rho$model == "CMORcor"){
# if (all(sapply(1:rho$ndim, function(j) all(is.na(rho$threshold.values[[j]]))))){ #all thresholds NA
# m <- "CMORcor"
# } else if (all(sapply(1:rho$ndim, function(j) (length(which(!is.na(rho$threshold.values[[j]])) >= 1) && all(which(!is.na(rho$threshold.values[[j]])) == 1))))){ #all first thresholds are not NA (and no additional)
# #warning if no intercept
# if (attr(terms.formula(rho$formula), "intercept") == 0) stop("formula needs intercept if you want to fix a threshold in model with corGeneral", call. = FALSE)
# m <- "CMORcor2"
# } else if (all(sapply(1:rho$ndim, function(j) all(!is.na(rho$threshold.values[[j]]))))){#all thresholds are fixed in all dimensions
# if (attr(terms.formula(rho$formula), "intercept") == 0) stop("formula needs intercept if you want to fix all thresholds in model with corGeneral", call. = FALSE)
# m <- "CMORcor"
# } else if (any(sapply(1:rho$ndim, function(j) all(!is.na(rho$threshold.values[[j]]))))){#all thresholds are fixed in at least one dimension
# ##CHECK somewhere if all values are fixed in one dimension + ?warning(sense)
# if (all(sapply(1:rho$ndim, function(j) all(!is.na(rho$threshold.values[[j]])) || all(
# (length(which(!is.na(rho$threshold.values[[j]])) >= 1) && (which(!is.na(rho$threshold.values[[j]])) == 1)))))){ #all not NA or only first fixed for each j
# if (attr(terms.formula(rho$formula), "intercept") == 0) stop("formula needs intercept if you want to fix a threshold in model with corGeneral", call. = FALSE)
# m <- "CMORcor2"
# } else if (all(sapply(1:rho$ndim, function(j) all(!is.na(rho$threshold.values[[j]])) || all(is.na(rho$threshold.values[[j]]))))){ #all not NA or all NA for each j
# if (attr(terms.formula(rho$formula), "intercept") == 0) stop("formula needs intercept if you want to fix a threshold in model with corGeneral", call. = FALSE)
# m <- "CMORcor"
# } else stop("Either fix all thresholds in one or more outcome dimensions,
# or consistently in all other outcome dimensions, all first thresholds or none.
# In addition, intercept must be included in the model.", call. = FALSE)
# } else{
# stop("Model is not identifiable. One of the following could help:
# 1. Remove intercept in formula (~ 0 + ...) or
# 2. Fix all first thresholds in threshold.values or
# 3. Fix all thresholds in threshold.values or
# 4. Either fix all thresholds in one or more outcome dimensions,
# or consistently in all other outcome dimensions, all first thresholds or none", call. = FALSE)
# }
# }
# if (rho$model == "CMORcov"){
# if (attr(terms.formula(rho$formula), "intercept") == 0) stop("formula needs intercept if you want to fix a threshold in model CMORcov", call. = FALSE)
# if (all(sapply(1:rho$ndim, function(j){
# (all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1,2))) && all(which(!is.na(rho$threshold.values[[j]])) == c(1,2)) #all first two thresholds are not NA
# ) || ((length(rho$threshold.values[[j]]) == 1) && (which(!is.na(rho$threshold.values[[j]])) == 1))
# }))){
# m <- "CMORcov"
# } else if (all(sapply(1:rho$ndim, function(j){
# (all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1,rho$ntheta[j]))) && all(which(!is.na(rho$threshold.values[[j]])) == c(1,rho$ntheta[j]))#all first and last two thresholds are not NA
# ) || ((length(rho$threshold.values[[j]]) == 1) && (which(!is.na(rho$threshold.values[[j]])) == 1))
# }))){
# m <- "CMORcov2"
# } else if (any(sapply(1:rho$ndim, function(j) all(!is.na(rho$threshold.values[[j]]))))){#all thresholds are fixed in at least one dimensions
# ##CHECK somewhere if all values are fixed in one dimensions + ?warning(sense)
# if (all(sapply(1:rho$ndim, function(j) #
# #(all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1))) && all((which(!is.na(rho$threshold.values[[j]])) == 1))) ||
# (all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1,2))) && all((which(!is.na(rho$threshold.values[[j]])) == c(1,2)))) ||
# (all(!is.na(rho$threshold.values[[j]])))))){
# m <- "CMORcov"
# } else if (all(sapply(1:rho$ndim, function(j)
# #(all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1))) && all((which(!is.na(rho$threshold.values[[j]])) == 1))) ||
# (all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1,rho$ntheta[j]))) && all((which(!is.na(rho$threshold.values[[j]])) == c(1,rho$ntheta[j])))) ||
# (all(!is.na(rho$threshold.values[[j]])))))){
# m <- "CMORcov2"
# } else stop("Either fix all thresholds in one or more outcome dimensions,
# or consistently in all other outcome dimensions, all first two thresholds or all first and last thresholds", call. = FALSE)
# } else stop("Either fix first two or first and last threshold in threshold.values if you use a model with covGeneral", call. = FALSE)
# }
# }
# }
# m
# }
Try the MultOrd package in your browser
Any scripts or data that you put into this service are public.
MultOrd documentation built on May 2, 2019, 4:49 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.
## TODO if no intercept in the covGeneral model set default to fix1first
set.threshold.type <- function(rho){
#fixall
if (all(sapply(1:rho$ndim, function(j) all(!is.na(rho$threshold.values[[j]]))))) {#all thresholds are fixed in all dimensions
if (rho$intercept == FALSE) cat("We suggest to include an intercept in the model (formula = y ~ 1 + ...)")
type <- "fixall"
#fix2first
} else if (all(sapply(1:rho$ndim, function(j){
(all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1,2))) && all(which(!is.na(rho$threshold.values[[j]])) == c(1,2)) #all first two thresholds are not NA
) || ((length(rho$threshold.values[[j]]) == 1) && (which(!is.na(rho$threshold.values[[j]])) == 1))
}))){
if (rho$error.structure$type %in% c("corGeneral", "corAR1", "corEqui")) cat("We suggest to fix only one threshold or the intercept in a correlation model.\n")
if ((rho$error.structure$type %in% c("covGeneral"))&& (rho$intercept.type == "fixed")) cat("We suggest to fix either two thresholds or one threshold and the intercept in a covGeneral model.\n")
type <- "fix2first"
#fix2firstlast
} else if (all(sapply(1:rho$ndim, function(j){
(all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1,rho$ntheta[j]))) && all(which(!is.na(rho$threshold.values[[j]])) == c(1,rho$ntheta[j]))#all first and last two thresholds are not NA
) || ((length(rho$threshold.values[[j]]) == 1) && (which(!is.na(rho$threshold.values[[j]])) == 1))
}))){
if (rho$error.structure$type %in% c("corGeneral", "corAR1", "corEqui")) cat("We suggest to fix only one threshold or the intercept in a correlation model.\n")
if ((rho$error.structure$type %in% c("covGeneral"))&& (rho$intercept.type == "fixed")) cat("We suggest to fix either two thresholds or one threshold and the intercept in a covGeneral model.\n")
type <- "fix2firstlast"
#fix1first
} else if (all(sapply(1:rho$ndim, function(j) (length(which(!is.na(rho$threshold.values[[j]])) >= 1) && all(which(!is.na(rho$threshold.values[[j]])) == 1))))){ #all first thresholds are not NA (and no additional)
if ((rho$error.structure$type == "covGeneral") && (rho$intercept.type == "flexible")) stop("Model with covGeneral is not identifiable.
Please either fix two thresholds or one threshold and the intercept.\n", call. = FALSE)
if ((rho$error.structure$type %in% c("corGeneral", "corAR1", "corEqui"))&& (rho$intercept.type == "fixed")) cat("We suggest to fix only one threshold or the intercept in a correlation model.\n")
type <- "fix1first"
#flexible
} else if (all(sapply(1:rho$ndim, function(j) all(is.na(rho$threshold.values[[j]]))))){#all thresholds NA
if (rho$error.structure$type == "covGeneral") stop("Model with covGeneral is not identifiable.
Please either fix two thresholds or one threshold and the intercept.\n", call. = FALSE)
if ((rho$error.structure$type %in% c("corGeneral", "corAR1", "corEqui")) && (rho$intercept.type == "flexible")) stop("Model is not identifiable.
Please either fix one threshold or the intercept.", call. = FALSE)
type <- "flexible"
#ERRORS
} else stop("Either fix all thresholds in one or more outcome dimensions,
or consistently in all other outcome dimensions, all first thresholds or none.\n", call. = FALSE)
if((rho$error.structure$type == "covGeneral") && (rho$binary == TRUE)){
if(type != "fix1first") stop("In the presence of binary outcomes intercept and one threshold
have to be fixed to some value.\n", call. = FALSE)
}
type
}
# #CORRELATION
# } else if (rho$error.structure$type %in% c("corGeneral", "corAR1", "corEqui")){
# #NO INTERCEPT
# if (rho$intercept.type == "fixed"){
# if (all(sapply(1:rho$ndim, function(j) all(is.na(rho$threshold.values[[j]]))))){#all thresholds NA
# type <- "flexible"
# } else stop("If you want fix a threshold value please use intercept in formula (remove ~ 0 + ...)", call. = FALSE)
# #TODO? other parameterization are not allowed for intercept == FALSE?
#
# #INTERCEPT
# } else if (rho$intercept == TRUE){
# if (all(sapply(1:rho$ndim, function(j) all(is.na(rho$threshold.values[[j]]))))){ #all thresholds NA
# if (all(!is.na(rho$coef.values[,1]))){ #intercepts are fixed to constants
# type <- "flexible"
# } else stop("Remove intercept in formula (~ 0 + ...) or fix the first threshold in threshold.values", call. = FALSE)
# }
# if (all(sapply(1:rho$ndim, function(j) (length(which(!is.na(rho$threshold.values[[j]])) >= 1) && all(which(!is.na(rho$threshold.values[[j]])) == 1))))){ #all first thresholds are not NA (and no additional)
# type = "fix1first"
# } else if (any(sapply(1:rho$ndim, function(j) all(!is.na(rho$threshold.values[[j]]))))){#all thresholds are fixed in at least one dimension
# if (all(sapply(1:rho$ndim, function(j) all(!is.na(rho$threshold.values[[j]])) || all(
# (length(which(!is.na(rho$threshold.values[[j]])) >= 1) && (which(!is.na(rho$threshold.values[[j]])) == 1)))))){ #all not NA or only first fixed for each j
# type <- "fix1first"
# } else if (all(sapply(1:rho$ndim, function(j) all(!is.na(rho$threshold.values[[j]])) || all(is.na(rho$threshold.values[[j]]))))){ #all not NA or all NA for each j
# type <- "flexible"
# } else stop("Either fix all thresholds in one or more outcome dimensions,
# or consistently in all other outcome dimensions, all first thresholds or none.
# In addition, intercept must be included in the model.", call. = FALSE)
# } else{
# stop("Model is not identifiable. One of the following could help:
# 1. Remove intercept in formula (~ 0 + ...) or
# 2. Fix all first thresholds in threshold.values or
# 3. Fix all thresholds in threshold.values or
# 4. Either fix all thresholds in one or more outcome dimensions,
# or consistently in all other outcome dimensions, all first thresholds or none", call. = FALSE)
# }
# }
#
# #covGeneral
# } else if (rho$error.structure$type == "covGeneral"){
# if (rho$intercept == FALSE) stop("formula needs intercept if you want to fix a threshold in model CMORcov", call. = FALSE)#TODO?
# if (all(sapply(1:rho$ndim, function(j){
# (all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1,2))) && all(which(!is.na(rho$threshold.values[[j]])) == c(1,2)) #all first two thresholds are not NA
# ) || ((length(rho$threshold.values[[j]]) == 1) && (which(!is.na(rho$threshold.values[[j]])) == 1))
# }))){
# type == "fix2first"
# } else if (all(sapply(1:rho$ndim, function(j){
# (all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1,rho$ntheta[j]))) && all(which(!is.na(rho$threshold.values[[j]])) == c(1,rho$ntheta[j]))#all first and last two thresholds are not NA
# ) || ((length(rho$threshold.values[[j]]) == 1) && (which(!is.na(rho$threshold.values[[j]])) == 1))
# }))){
# type == "fix2firstlast"
# }else if (any(sapply(1:rho$ndim, function(j) all(!is.na(rho$threshold.values[[j]]))))){#all thresholds are fixed in at least one dimensions
# ##CHECK somewhere if all values are fixed in one dimensions + ?warning(sense)
# if (all(sapply(1:rho$ndim, function(j) #
# #(all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1))) && all((which(!is.na(rho$threshold.values[[j]])) == 1))) ||
# (all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1,2))) && all((which(!is.na(rho$threshold.values[[j]])) == c(1,2)))) ||
# (all(!is.na(rho$threshold.values[[j]])))))){
# type <- "fix2first"
# } else if (all(sapply(1:rho$ndim, function(j)
# #(all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1))) && all((which(!is.na(rho$threshold.values[[j]])) == 1))) ||
# (all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1,rho$ntheta[j]))) && all((which(!is.na(rho$threshold.values[[j]])) == c(1,rho$ntheta[j])))) ||
# (all(!is.na(rho$threshold.values[[j]])))))){
# type <- "fix2firstlast"
# } else stop("Either fix all thresholds in one or more outcome dimensions,
# or consistently in all other outcome dimensions, all first two thresholds or all first and last thresholds", call. = FALSE)
# } else stop("Either fix first two or first and last threshold in threshold.values if you use a model with covGeneral", call. = FALSE)
# }
# type
# }
#data.x <- data.multord$x
#ndim <- rho$ndim
set.error.structure <- function(error.structure, data.x, ndim){
if (error.structure$formula == ~1 ) {
if (error.structure$type %in% c("corEqui", "corAR1")) {
error.structure$x <- model.matrix(error.structure$formula, data.x[[1]])
} else {
error.structure$x <- as.factor(model.matrix(error.structure$formula, data.x[[1]]))
}
} else if (error.structure$type %in% c("corEqui", "corAR1")){
tmp <- lapply(1:ndim, function(j) {
tmp1 <- model.matrix(error.structure$formula, data.x[[j]])
attribute <- attr(tmp1, "assign")
tmp2 <- tmp1[match(rownames(data.x[[j]]),rownames(tmp1)),]
rownames(tmp2) <- rownames(data.x[[j]])
attr(tmp2, "assign") <- attribute
tmp2
}
)
error.x <- sapply(1:ncol(tmp[[1]]), function(k){
xtcol <- do.call(cbind,lapply(tmp, `[`,,k))
xtcol_final <- apply(xtcol,1,function(i) unique(i[!is.na(i)]))
if(is.list(xtcol_final)) stop("Covariate dependent variables need to be constant across multiple observations", call.=FALSE)
xtcol_final
})
colnames(error.x) <- colnames(tmp[[1]])
attr(error.x, "assign") <- attr(tmp[[1]], "assign")
error.structure$x <- error.x
# } else if (error.structure$type == "corAR1"){
# #depends on one factor (factor of first year)
# if (!is.factor(data.x[[1]][,all.vars(error.structure$formula)])) stop("variable must be of type factor in corAR1", call. = FALSE)
# tmp <- lapply(1:ndim,function(j) data.x[[j]][,all.vars(error.structure$formula)])
# if (!all(sapply(1:ndim, function(j) all(tmp[[1]]==tmp[[j]], na.rm=T)))) stop("Covariate dependent variables need to be constant across multiple observations", call. = FALSE)
# error.structure$x <- tmp[[1]]#data.x[[1]][,all.vars(error.structure$formula)]
} else if (error.structure$type %in% c("corGeneral", "covGeneral")){
#depends on one factor
if (length(all.vars(error.structure$formula[[2]])) > 1) stop("only one factor is allowed in covGeneral and corGeneral", call. = FALSE)
if (!is.factor(data.x[[1]][,all.vars(error.structure$formula)])) stop("variable must be of type factor in covGeneral and corGeneral", call. = FALSE)
tmp <- lapply(1:ndim,function(j) data.x[[j]][,all.vars(error.structure$formula)])
if (!all(sapply(1:ndim, function(j) all(tmp[[1]]==tmp[[j]], na.rm=T)))) stop("Covariate dependent variables need to be constant across multiple observations", call. = FALSE)
error.structure$x <- as.factor(apply(do.call("cbind.data.frame", tmp), 1, function(x) unique(x[!is.na(x)])))
} else {
#depends on one factor
if (length(all.vars(error.structure$formula[[2]])) > 1) stop("only one factor is allowed in covGeneral and corGeneral", call. = FALSE)
if (!is.factor(data[,all.vars(error.structure$formula)])) stop("variable must be of type factor in covGeneral and corGeneral", call. = FALSE)
error.structure$x <- data.x[,all.vars(error.structure$formula)]
}
error.structure
}
set.error.structure.CMOR <- function(error.structure, data){
if (error.structure$type == "corEqui"){
error.structure$x <- model.matrix(error.structure$formula, data)
} else if (error.structure$type == "corAR1"){
stop("corAR1 is not applicable in CMOR", call. = FALSE)
# if (error.structure$formula == ~1){
# error.structure$x <- as.factor(model.matrix(error.structure$formula, data[[1]]))
# } else{ #depends on one factor (factor of first year)
# if (!is.factor(data[[1]][,all.vars(error.structure$formula)])) stop("variable must be of type factor in covGeneral and corGeneral", call. = FALSE)
# error.structure$x <- data[[1]][,all.vars(error.structure$formula)]
# }
} else if (error.structure$type %in% c("corGeneral", "covGeneral")){
if (error.structure$formula == ~ 1){
error.structure$x <- as.factor(model.matrix(error.structure$formula, data))
} else{ #depends on one factor
if (length(all.vars(error.structure$formula[[2]])) > 1) stop("only one factor is allowed in covGeneral and corGeneral", call. = FALSE)
if (!is.factor(data[,all.vars(error.structure$formula)])) stop("variable must be of type factor in covGeneral and corGeneral", call. = FALSE)
error.structure$x <- data[,all.vars(error.structure$formula)]
}
}
error.structure
}
checkArgs <- function(rho){
#CHECK if treshold.values is in line with threshold.constraints
if (length(rho$threshold.constraints) != rho$ndim) stop("dimensions of threshold.values and number of thresholds do not match", call. = FALSE)
if (any(sapply(1:rho$ndim, function(j) length(rho$threshold.values[[j]]) != rho$ntheta[j])))
stop("dimensions of threshold.values and number of thresholds do not match", call. = FALSE)
for (j in unique(rho$threshold.constraints)){
ind <- which(rho$threshold.constraints == j)
if (length(unique(rho$threshold.values[ind]))!=1) stop("If constraints are set on thresholds (by threshold.constraints), threshold.values need to be specified accordingly for these outcome dimensions.
Maybe dimension do not have the same number of threshold parameters.", call. = FALSE)
}
if (nrow(rho$coef.constraints) != rho$ndim) stop("row dimension of coef.constraints and outcome dimension do not match (?factor)", call. = FALSE)
for (j in 1:ncol(rho$coef.constraints)){
indj <- unique(rho$coef.constraints[,j])
indj <- indj[!is.na(indj)]
lapply(seq_len(length(indj)), function(k) {
tmpind <- which(rho$coef.constraints[,j] == indj[k])
tmp <- rho$coef.values[tmpind,j]
if(length(unique(tmp)) != 1) stop("If constraints are set on the coefficients (by coef.constraints),
coef.values need to be specified accordingly for these outcome dimensions.", call. = FALSE)
})
}
}
##########################################
###### AUXILIARY FUNCTIONS ##############
###########################################################################
# rectbiv.norm.prob_TEST <- function(U, L, r) {
# # computes the rectangle probabilities for biv.normal-distribution
# if (is.null(dim(U)[2]) & is.null(dim(L)[2])) {
# dim(U) <- dim(L) <- c(1, length(U))
# }
# n <- nrow(U)
# x <- c(U[, 1],L[, 1],U[, 1],L[, 1])
# y <- c(U[, 2],U[, 2],L[, 2],L[, 2])
# r
# p <- pbivnorm(x,y,r)
# p[is.nan(p)] <- 0
# p[seq_len(n)]-p[n + seq_len(n)]-p[2*n + seq_len(n)] + p[3*n + seq_len(n)]
# }
rectbiv.norm.prob <- function(U, L, r) {
# computes the rectangle probabilities for biv.normal-distribution
if (is.null(dim(U)[2]) & is.null(dim(L)[2])) {
dim(U) <- dim(L) <- c(1, length(U))
}
p1 <- pbivnorm(U[, 1], U[, 2], r)
p2 <- pbivnorm(L[, 1], U[, 2], r)
p3 <- pbivnorm(U[, 1], L[, 2], r)
p4 <- pbivnorm(L[, 1], L[, 2], r)
p1[is.nan(p1)] <- 0
p2[is.nan(p2)] <- 0
p3[is.nan(p3)] <- 0
p4[is.nan(p4)] <- 0
p1-p2-p3+p4
}
biv.nt.prob2 <-function (df, lower, upper,
mean=c(0,0), r) {
if (df == Inf) {
nu <- 0
} else {
nu <- df
}
rho <- r
if (any(lower > upper)) stop("lower > upper integration limits")
if (any(lower == upper)) return(0)
infin <- as.integer(c(2, 2))
prob <- as.double(0)
a <- .Fortran("smvbvt", prob, nu, lower, upper, infin, rho,
#PACKAGE = "mnormt")
PACKAGE = "MultOrd")
return(a[[1]])
}
casesNA<- function(y) {
# returns a column vector v, assigning a different number/label to each unique combination of NA pattern
contrast <- as.matrix((!is.na(y)) + 0)
contrast.char <- apply(contrast, 1, paste, collapse = "") # make each contrast row as character i.e., 111 has all three response variables
factor(contrast.char, levels = unique(contrast.char))
}
split.NA.pattern <- function(y) {
# y should be a data.frame
# return a list, whose elements are groups of the matrix y according to each unique combination of response variables
if (!is.data.frame(y)) y <- as.data.frame(y)
yg <- split(y, casesNA(y)) # grouped y by cases of NA pattern
lapply(yg, function(x) which(rownames(y) %in% rownames(x)))
}
getStart.values <- function(rho){
# # set starting value for beta and the gammas thresholds
# #thetas <- qlogis((1:rho$ntheta) / (rho$ntheta + 1))
# #if (rho$link != "logistic")
# # thetas <- thetas/1.7
# mclm <- lapply(1:rho$ndim, function(j)
# ordinal:::clm(factor(rho$y[,j]) ~ rho$x[,which(rho$coef.constraints[j,] != 0)], link = "logit"))
# #thetas <- lapply(mclm, "[[", "alpha")
# betas <- lapply(mclm, "[[", "beta") #Reduce("+", lapply(mclm, "[[", "beta"))
# #gammas <- lapply(thetas, function(x) c(x[1L], log(diff(x))))
gammas <- sapply(1:rho$ndim, function(j) {
if (rho$npar.theta.opt[j] != 0){
theta <- if (rho$ntheta[j] >= 2) polr(rho$y[, j] ~1)$zeta else 0#qlogis((1:rho$ntheta[j])/(rho$ntheta[j] + 1))
if (!grepl("logit", rho$link)) theta <- theta/1.7
c(theta[1L], log(diff(theta)))[1:rho$npar.theta.opt[j]]
} else NULL
})
c(unlist(gammas))#, unlist(betas))# rep(, sum(rho$npar.betas))) }
}
transf.par.cor <- function(par, rho) {
sigmas <- rho$transf.sigmas(par[rho$npar.thetas + rho$npar.betas + seq_len(rho$npar.sigmas)],rho)
theta <- rho$transf.thresholds(par[seq_len(rho$npar.thetas)], rho)
par_beta <- par[rho$npar.thetas + seq(rho$npar.betas)]
beta <- sapply(1:ncol(rho$coef.constraints), function(j){
sapply(1:nrow(rho$coef.constraints), function(i,j) ifelse(is.na(rho$ind.coef[i,j]), rho$coef.values[i,j], par_beta[rho$ind.coef[i, j]]), j)
})
pred.fixed <- sapply(1:rho$ndim, function(j) rho$x[[j]] %*% beta[j, ])
theta.lower <- sapply(1:rho$ndim, function(j) c(-rho$inf.value, theta[[j]])[rho$y[, j]])
theta.upper <- sapply(1:rho$ndim, function(j) c(theta[[j]], rho$inf.value)[rho$y[, j]])
pred.lower <- (theta.lower - pred.fixed)/rho$sd.y
pred.upper <- (theta.upper - pred.fixed)/rho$sd.y
list(U = pred.upper,
L = pred.lower,
sigmas = sigmas)
}
transf.par.cov <- function(par, rho) {
exp.par.sd <- exp(par[rho$npar.thetas + rho$npar.betas +
rho$npar.cor * rho$ncor.levels +
seq_len(rho$npar.cor.sd * rho$ncor.levels)])
sigmas <- rho$transf.sigmas(par[rho$npar.thetas + rho$npar.betas + seq_len(rho$npar.sigmas)],rho,
exp.par.sd)
exp.par.sd <- lapply(1:rho$ncor.levels, function(l) exp.par.sd[ (l-1) * rho$npar.cor.sd + seq_len(rho$npar.cor.sd) ])
theta <- rho$transf.thresholds(par[seq_len(rho$npar.thetas)], rho)
par_beta <- par[rho$npar.thetas + seq_len(rho$npar.betas)]
beta <- sapply(1:ncol(rho$coef.constraints), function(j){
sapply(1:nrow(rho$coef.constraints), function(i,j)
ifelse(is.na(rho$ind.coef[i,j]), rho$coef.values[i,j], par_beta[rho$ind.coef[i, j]]), j)
})
pred.fixed <- sapply(1:rho$ndim, function(j) rho$x[[j]] %*% beta[j, ])
theta.lower <- sapply(1:rho$ndim, function(j) c(-rho$inf.value, theta[[j]])[rho$y[, j]])
theta.upper <- sapply(1:rho$ndim, function(j) c(theta[[j]], rho$inf.value)[rho$y[, j]])
## make matrix of exp.par.sd
lev <- match(rho$error.structure$x, rho$error.structure$levels)
#pred.lower <- (theta.lower - pred.fixed)/rho$sd.y
#pred.upper <- (theta.lower - pred.fixed)/rho$sd.y
#exp.par.sd.mat <- matrix(ncol = ncol(theta.upper), nrow = nrow(theta.upper))
#for (l in 1:rho$ncor.levels){
# exp.par.sd.mat[lev == l, ] <-
# exp.par.sd[(l - 1) * rho$npar.cor.sd + seq_len(rho$npar.cor.sd)]
#}
#for (i in 1:nrow(pred.lower)){
# pred.lower[i,] <- pred.lower[i,]/exp.par.sd[(lev[i] - 1) * rho$npar.cor.sd + seq_len(rho$npar.cor.sd)]
# pred.upper[i,] <- pred.upper[i,]/exp.par.sd[(lev[i] - 1) * rho$npar.cor.sd + seq_len(rho$npar.cor.sd)]
#}
pred.lower <- theta.lower - pred.fixed
pred.upper <- theta.upper - pred.fixed
for (l in 1:rho$ncor.levels){
pred.lower[lev==l,] <-t(t(pred.lower[lev==l,]/rho$sd.y)/exp.par.sd[[l]])
pred.upper[lev==l,] <-t(t((pred.upper[lev==l,])/rho$sd.y)/exp.par.sd[[l]])
}
#pred.lower <-t(t((theta.lower - pred.fixed)/rho$sd.y)/exp.par.sd)
#pred.upper <-t(t((theta.upper - pred.fixed)/rho$sd.y)/exp.par.sd)
# t(t(pred.upper)/exp.par.sd)
# test speed up
sigmas <- lapply(sigmas, cov2cor)#cov2cor(sigmas)
list(U = pred.upper,
L = pred.lower,
sigmas = sigmas,
std.dev = exp.par.sd)
}
#########################################################################
## transformation of the threshold parameters (to ensure monotonicity) ##
#########################################################################
transf.thresholds.fixall <- function(gamma, rho){
rho$threshold.values
}
transf.thresholds.flexible <- function(gamma, rho){
lapply(1:rho$ndim, function(j)
if (anyNA(rho$threshold.values[[j]])){
if (rho$ntheta[j] > 1) {
cumsum(c(gamma[rho$ind.thresholds[[j]][1]],
exp(gamma[rho$ind.thresholds[[j]][2:rho$ntheta[j]]])))
} else if (rho$ntheta[j] == 1) gamma[rho$ind.thresholds[[j]]] else NULL
} else rho$threshold.values[[j]]
)
}
transf.thresholds.fix1.first <- function(gamma, rho){
lapply(1:rho$ndim, function(j)
if (length(rho$ind.thresholds[[j]]) >= 1) {
cumsum(c(rho$threshold.values.fixed[[j]][1],
exp(gamma[rho$ind.thresholds[[j]][seq_len(length(rho$ind.thresholds[[j]]))]])))
} else if (length(rho$ind.thresholds[[j]]) == 0) rho$threshold.values.fixed[[j]] else NULL)
}
#gamma <- par[seq_len(rho$npar.thetas)]
transf.thresholds.fix2.first <- function(gamma, rho){
lapply(1:rho$ndim, function(j)
if (rho$npar.theta[j] >= 1) {
c(rho$threshold.values.fixed[[j]][1],
cumsum(c(rho$threshold.values.fixed[[j]][2] ,exp(gamma[rho$first.ind.theta[j] - 1 +
seq_len(rho$npar.theta[j])]))))
} else {
if (rho$npar.theta[j] == 0) rho$threshold.values.fixed[[j]]
}
)
}
# transf.thresholds.fix2.first <- function(gamma, rho){
# lapply(1:rho$ndim, function(j)
# if (rho$npar.theta[rho$threshold.constraints[j]] >= 1) {
# c(rho$threshold.values.fixed[[j]][1],
# cumsum(c(rho$threshold.values.fixed[[j]][2] ,exp(gamma[rho$first.ind.theta[rho$threshold.constraints[j]] - 1 +
# seq_len(rho$npar.theta[rho$threshold.constraints[j]])]))))
# } else {
# if (rho$npar.theta[rho$threshold.constraints[j]] == 0) rho$threshold.values.fixed[[j]]
# }
# )
# }
transf.thresholds.fix2.firstlast <- function(gamma, rho){
lapply(1:rho$ndim, function(j){
gamma1 <- gamma[rho$first.ind.theta[j] - 1 + seq_len(rho$npar.theta[j])]
recursive.theta <- function(i) {
if (i == 0) 0
else return ((exp(gamma1[i]) +
recursive.theta(i - 1))/(1 + exp(gamma1[i])))
}
if (rho$npar.theta[j] >= 1) {
theta <- sapply(1:rho$npar.theta[j], function(i)
recursive.theta(i))
c(0, theta, 1) * (rho$threshold.values.fixed[[j]][2] - rho$threshold.values.fixed[[j]][1]) +
rho$threshold.values.fixed[[j]][1]
} else {
rho$threshold.values.fixed[[j]]
}
})
}
# transf.thresholds.fix2.firstlast <- function(gamma, rho){
# lapply(1:rho$ndim, function(j){
# gamma1 <- gamma[rho$first.ind.theta[j] - 1 + seq_len(rho$npar.theta[rho$threshold.constraints[j]])]
# recursive.theta <- function(i) {
# if (i == 0) 0
# else return ((exp(gamma1[i]) +
# recursive.theta(i - 1))/(1 + exp(gamma1[i])))
# }
# if (rho$npar.theta[rho$threshold.constraints[j]] >= 1) {
# theta <- sapply(1:rho$npar.theta[rho$threshold.constraints[j]], function(i)
# recursive.theta(i))
# c(0, theta, 1) * (rho$threshold.values.fixed[[j]][2] - rho$threshold.values.fixed[[j]][1]) +
# rho$threshold.values.fixed[[j]][1]
# } else {
# rho$threshold.values.fixed[[j]]
# }
# })
# }
##############################################################################
## transformation of the parameters of the correlation/covariance structure ##
##############################################################################
transf.sigmas.spheric <- function(par.nu, rho, exp.par.sd = NULL) {
lapply(1:rho$ncor.levels, function(l) {
nu <- par.nu[(l - 1) * rho$npar.cor + seq_len(rho$npar.cor)]
angles <- pi * exp(nu)/(1 + exp(nu))
cosmat <- diag(rho$ndim)
cosmat[lower.tri(cosmat)] <- cos(angles)
S1 <- matrix(0, nrow = rho$ndim, ncol = rho$ndim)
S1[, 1L] <- 1
S1[-1L, -1L][lower.tri(S1[-1L, -1L], diag = T)] <- sin(angles)
tLmat <- sapply(1:rho$ndim,
function(j) cosmat[j, ] * cumprod(S1[j, ]))
sigma <- crossprod(tLmat)
if (is.null(exp.par.sd)) {
return(sigma)
} else {
stdev <- exp.par.sd[(l - 1) * rho$npar.cor.sd + seq_len(rho$npar.cor.sd)]
return(t(stdev * sigma) * stdev)
}
})
}
# transf.sigmas.corAR1 <- function(par.sigma, rho, exp.par.sd = NULL) {
# lapply(1:rho$ncor.levels, function(l){
# ###TODO include covariate dependent corAR1? (like corEqui)
# r <- z2r(par.sigma[seq_len(rho$npar.cor) + l -1])
# #r <- 0.9
# sigma <- diag(rho$ndim)
# sigma[lower.tri(sigma)] <- r^sequence((rho$ndim-1):1)
# ### CHECK: easier way?
# sigma <- sigma + t(sigma) - diag(diag(sigma))
# sigma
# })
# }
transf.sigmas.corAR1 <- function(par.sigma, rho, exp.par.sd = NULL) {
w <- par.sigma[seq_len(rho$npar.cor)]
z <- rho$error.structure$x %*% w
r <- z2r(z)
lapply(seq_len(length(r)), function(i){
sigma <- diag(rho$ndim)
sigma[lower.tri(sigma)] <- r[i]^sequence((rho$ndim-1):1)
### CHECK: easier way?
sigma <- sigma + t(sigma) - diag(diag(sigma))
sigma
})
}
transf.sigmas.corEqui <- function(par.sigma, rho, exp.par.sd = NULL) {
w <- par.sigma[seq_len(rho$npar.cor)]
z <- rho$error.structure$x %*% w
z2r(z)
}
z2r <- function (z) {
ifelse(z > 354, 1, (exp(2 * z) - 1)/(1 + exp(2 * z)))
}
getInd.coef <- function(coef.constraints, coef.values) {
ind <- matrix(NA,ncol = ncol(coef.constraints), nrow = nrow(coef.constraints))
k <- 1
i <- 1
for(j in 1:ncol(coef.constraints)){
if (is.na(coef.values[i,j])){
ind[i,j] <- k
k <- k + 1
}
}
for(i in 2:nrow(coef.constraints)){
for(j in 1:ncol(coef.constraints)){
if (is.na(coef.values[i,j]) && !(coef.constraints[i,j] %in% coef.constraints[1:(i-1),j])){
ind[i,j] <- k
k <- k + 1
}
if (is.na(coef.values[i,j]) && coef.constraints[i,j] %in% coef.constraints[1:(i-1),j]){
ind[i,j] <- ind[min(which(coef.constraints[i,j] == coef.constraints[1:(i-1),j])),j]
}
}
}
ind
}
getInd.thresholds.flexible <- function(threshold.constraints,rho){
pointer <- lapply(1:rho$ndim, function(j) rep(0,rho$npar.theta[j]))
pointer[[1]] <- seq_len(rho$npar.theta[1])
k <- length(pointer[[1]])
for(i in 2:rho$ndim){
if (rho$npar.theta[i] == 0){
pointer[[i]] <- integer(0)
} else {
if (threshold.constraints[i] %in% threshold.constraints[1:(i-1)]){
min.ind <- min(which(threshold.constraints[i] == threshold.constraints[1:(i-1)]))
if (rho$ntheta[i] != rho$ntheta[min.ind])
stop("Constraints on threshold parameters are not valid (different number of categories)", call. = FALSE)
pointer[[i]] <- pointer[[min.ind]]
} else{
pointer[[i]] <- (k+1):(k+rho$ntheta[i])
k <- k + rho$ntheta[i]
}
}
}
pointer
}
getInd.thresholds.fix2 <- function(threshold.constraints,rho) {
#pointer <- lapply((1:rho$ndim)[which(rho$npar.theta > 0)], function(j) 1:rho$npar.theta[j])
pointer <- lapply((1:rho$ndim), function(j) seq_len(rho$npar.theta[j]))
k <- length(pointer[[1]])
#for(i in (1:rho$ndim)[which(rho$npar.theta > 0)][-1]){
for(i in (1:rho$ndim)){
if (rho$npar.theta[i] == 0){
pointer[[i]] <- integer(0)
} else {
if (threshold.constraints[i] %in% threshold.constraints[1:(i - 1)]){
min.ind <- min(which(threshold.constraints[i] == threshold.constraints[1:(i-1)]))
if (rho$ntheta[i] != rho$ntheta[min.ind])
stop("Constraints on threshold parameters are not valid (different number of categories)", call. = FALSE)
pointer[[i]] <- pointer[[min.ind]]
} else{
pointer[[i]] <- (k+1):(k+rho$ntheta[i]-2)
k <- k + rho$ntheta[i]-2
}
}
}
pointer
}
getInd.thresholds.fix1 <- function(threshold.constraints,rho) {
pointer <- lapply(1:rho$ndim, function(j) rep(0,rho$npar.theta[j]))
pointer[[1]] <- seq_len(rho$npar.theta[1])#1:(rho$ntheta[1]-1)
k <- length(pointer[[1]])
for(i in 2:rho$ndim){
if (rho$npar.theta[i] == 0){
pointer[[i]] <- integer(0)
} else {
if (threshold.constraints[i] %in% threshold.constraints[1:(i - 1)]){
min.ind <- min(which(threshold.constraints[i] == threshold.constraints[1:(i - 1)]))
if (rho$ntheta[i] != rho$ntheta[min.ind])
stop("Constraints on threshold parameters are not valid (different number of categories)", call. = FALSE)
pointer[[i]] <- pointer[[min.ind]]
} else{
pointer[[i]] <- (k+1):(k+rho$ntheta[i]-1)
k <- k + rho$ntheta[i]-1
}
}
}
pointer
}
getInd.thresholds.fixall <- function(threshold.constraints,rho) {
lapply(1:rho$ndim, function(j) integer())
}
get.labels.theta <- function(rho,j) {
lev <- levels(rho$y[, j])
sapply(1:(rho$ntheta[j]), function(i){
paste(lev[i], lev[i + 1], sep = "|")
})
}
#' @title Error Structures in Multord
#' @description Different \code{error.structures} are available in \pkg{MultOrd}:
#' \itemize{
#' \item general correlation structure (default) \code{corGeneral(~1)},
#' \item general covariance structure \code{covGeneral(~1)},
#' \item factor dependent correlation structure \code{covGeneral(~f)},
#' \item factor dependent covariance structure \code{covGeneral(~f)},
#' \item covariate dependent equicorrelation structure \code{corEqui(~X)},
#' \item AR(1) correlation structure \code{corAR1(~1)}, or
#' \item factor dependent AR(1) correlation structure \code{corAR1(~f)}.
#' }
#' See \code{\link{error.structures}} or 'Details'.
#' @details General (symmetric) correlation structure (cross-sectional)
#' @param formula \code{\link{formula}} object
#' @export
#' @name error.structures
corGeneral <- function(formula) UseMethod("corGeneral")
#' @rdname error.structures
#' @export
corGeneral.formula <- function(formula){
return(list(type = "corGeneral", formula = formula))
}
#' @details Equicorelation structure (cross-sectional)
# #' @param formula \code{\link{formula}} object
#' @export
#' @rdname error.structures
corEqui <- function(formula) UseMethod("corEqui")
#' @rdname error.structures
#' @export
corEqui.formula <- function(formula){
return(list(type = "corEqui", formula = formula))
}
#' @details General (symmetric) covariance structure (cross-sectional)
# #' @param formula \code{\link{formula}} object
#' @export
#' @rdname error.structures
covGeneral <- function(formula) UseMethod("covGeneral")
#' @rdname error.structures
#' @export
covGeneral.formula <- function(formula){
return(list(type = "covGeneral", formula = formula))
}
#' @details AR(1) correlation structure (panel; only applicable in line with \code{\link{multord}})
# #' @param formula \code{\link{formula}} object
#' @export
#' @rdname error.structures
corAR1 <- function(formula) UseMethod("corAR1")
#' @rdname error.structures
#' @export
corAR1.formula <- function(formula){
return(list(type = "corAR1", formula = formula))
}
backtransf.sigmas <- function(R){
J <- nrow(R)
l <- t(chol(R))
angmat <- matrix(1,ncol=J,nrow=J)
angmat[-1,1] <- acos(l[-1,1])
for (j in 2:(J-1)){
sinprod <- apply(sin(angmat[, seq_len(j-1), drop=F]), 1, prod) ## denominator in division
angmat[-(1:j),j]<-acos((l/sinprod)[-(1:j),j])
}
angdivpi <- angmat[lower.tri(angmat)]/pi
log(angdivpi/(1-angdivpi))
}
get.sigma.i <- function(object) {
if (object$rho$error.structure$type == "corEqui") {
lapply(1:object$rho$n, function(i) {
tmp <- matrix(object$r[i], nrow = object$rho$ndim, ncol = object$rho$ndim)
diag(tmp) <- 1
tmp
})
} else {
lev <- match(object$rho$error.structure$x, object$rho$error.structure$levels)
lapply(1:object$rho$n, function(i) object$sigmas[lev[i]][[1]])
}
}
#' @title Data preparation for multord
#'
#' @description
#' This function is an (internally) used to transforms the \code{data}, into a "multivariate setting",
#' where all repeated measurements are matched accordingly to their ID. A matrix of all ordinal responses with \code{J} columns
#' as well as a list of length \code{J} of matrices with all the covariates are created.
#' @param data a \code{data.frame}, where each row corresponds to a single measurement.
#' @param index is an (optional) argument that specifies the index for the subjects and the response index of the multiple measurement.
#' This is usually performed
#' by a character vector of length two specifying the column names of the subject index and
#' the multiple response index in data. The default value of index is NULL assuming that the
#' first column of data contains the subject index and the second column the multiple response index.
#' @param y.names column name of \code{data} where the ordinal observations are stored.
#' @param x.names column names of all the covariates in {data}.
#' @param y.levels (optional) list of length \code{J} that specifies the levels of each repeated measurement. If the categories
#' differ across repeated measurements (either the number of categories or the category labels) it is recommended to set them.
#' @param response.names (optional) vector of names of the repeated measurements in \code{data}
#'which specifies the ordering of the repeated measurements.
#' @export
#y.names <- "yy"
#y.names <- rho$response.name
#x.names <- rho$x.names
# data <- data_multord
# y.names <- rho$response.name
# x.names <- unique(c(rho$x.names,weights))
# y.levels = response.levels
# response.names = response.names = c("rater1", "rater2", "rater3")
# index = c("firmID", "rater")
multord.data <- function(data, index, y.names, x.names,
y.levels = NULL, response.names = NULL) {
df <- list()
if (any(duplicated(data[,index]))) stop("duplicated observation(s) for one index", call. = FALSE)
index.levels <- levels(as.factor(data[, index[2]]))
data.split <- split(data[,c(y.names, index[1])], data[, index[2]])
data.split.y <- lapply(seq_len(length(data.split)), function(j) {
colnames(data.split[[j]]) <- c(index.levels[j],index[1])
data.split[[j]]})
df$y <- Reduce(function(...) merge(..., by = index[1], all = TRUE), data.split.y )[, -1]
response.names.NA <- c()
if (!is.null(response.names)) {
response.names.NA <- response.names[!response.names %in% colnames(df$y)]
df$y <- cbind(df$y, matrix(NA, ncol = length(response.names.NA), nrow = nrow(df$y), dimnames = list(c(),response.names.NA)))
df$y <- df$y[,as.character(response.names)]
}
colnames.y <- colnames(df$y)
if (is.null(y.levels)) {
df$y <- do.call(cbind.data.frame, lapply(1:ncol(df$y), function(j) ordered(df$y[, j])))
} else {
df$y <- do.call(cbind.data.frame, lapply(1:ncol(df$y), function(j) {
if (!all(unique(df$y[!is.na(df$y[, j]), j]) %in% y.levels[[j]])) stop("levels of response do not match with y.levels", call. = FALSE)
ordered(df$y[, j], levels = y.levels[[j]])}
))
}
colnames(df$y) <- colnames.y#paste0(y.names, 1:ncol(df$y))
#X
data.split.x <- split(data[, c(x.names, index[1])], data[, index[2]])
names.x <- names(data.split.x)
#REMOVE this
#names.x <- names(data.split.x)[names(data.split.x) %in% response.names]
#set colnames (otherwise warning due to identical colnames in reduce)
data.split.x <- lapply(seq_len(length(data.split.x)), function(j) {
colnames(data.split.x[[j]]) <- c(paste(x.names,j, sep = "."), index[1])
data.split.x[[j]]})
xdatadf <- Reduce(function(...) merge(...,by = index[1], all = TRUE), data.split.x)[, -1]
xdatadf <- cbind(xdatadf, matrix(NA, ncol = length(response.names.NA) * length(x.names),
nrow = nrow(xdatadf)))
#xdatadf <- Reduce(function(x,y) merge(x,y,by = index[1], all = TRUE), split(data[,c(x.names, index[1])], data[, index[2]]) )[, -1]
#REMOVE this
# df$x <- lapply(1:ncol(df$y), function(i) {
# tmp <- xdatadf[,(i - 1) * length(x.names) + seq_len(length(x.names))]
# names(tmp) <- x.names
# tmp
# })
df$x <- lapply(1:(length(index.levels) + length(response.names.NA)), function(i) {
tmp <- xdatadf[,(i - 1) * length(x.names) + seq_len(length(x.names))]
names(tmp) <- x.names
tmp
})
names(df$x) <- c(names.x, response.names.NA)
if (!is.null(response.names)) df$x <- df$x[as.character(response.names)]
df
}
check <- function(...){
stopifnot(...)
}
#threshold.values <- list(c(NA,NA,NA),c(2,3,NA),c(1,2,NA,NA,NA))
# set.model <- function(rho){
# if (rho$error.structure$type == "corAR1"){
# m <- "PMOR"
# } else{
# if (rho$error.structure$type == "corGeneral"){
# #if (length(all.vars(rho$error.structure$formula[[2]])) > 1) stop("only one factor is allowed in corGeneral", call. = FALSE)
# rho$model <- "CMORcor"
# }
# if (rho$error.structure$type == "corEqui"){
# rho$model <- "CMORcor"
# }
# if (rho$error.structure$type == "covGeneral"){
# #if (length(all.vars(rho$error.structure$formula[[2]])) > 1) stop("only one factor is allowed in covGeneral", call. = FALSE)
# rho$model <- "CMORcov"
# }
# if (is.null(rho$threshold.values)){
# if ((rho$model == "CMORcor") && (attr(terms.formula(rho$formula), "intercept") == 1)){
# stop("Remove intercept in formula (~ 0 + ...) or fix the first threshold in threshold.values", call. = FALSE)
# }
# if ((rho$model == "CMORcov") && (attr(terms.formula(rho$formula), "intercept") == 0)){
# stop("Include intercept in formula (~ 1 + ...)", call. = FALSE)
# }
# m <- rho$model
# } else{
# if (rho$model == "CMORcor"){
# if (all(sapply(1:rho$ndim, function(j) all(is.na(rho$threshold.values[[j]]))))){ #all thresholds NA
# m <- "CMORcor"
# } else if (all(sapply(1:rho$ndim, function(j) (length(which(!is.na(rho$threshold.values[[j]])) >= 1) && all(which(!is.na(rho$threshold.values[[j]])) == 1))))){ #all first thresholds are not NA (and no additional)
# #warning if no intercept
# if (attr(terms.formula(rho$formula), "intercept") == 0) stop("formula needs intercept if you want to fix a threshold in model with corGeneral", call. = FALSE)
# m <- "CMORcor2"
# } else if (all(sapply(1:rho$ndim, function(j) all(!is.na(rho$threshold.values[[j]]))))){#all thresholds are fixed in all dimensions
# if (attr(terms.formula(rho$formula), "intercept") == 0) stop("formula needs intercept if you want to fix all thresholds in model with corGeneral", call. = FALSE)
# m <- "CMORcor"
# } else if (any(sapply(1:rho$ndim, function(j) all(!is.na(rho$threshold.values[[j]]))))){#all thresholds are fixed in at least one dimension
# ##CHECK somewhere if all values are fixed in one dimension + ?warning(sense)
# if (all(sapply(1:rho$ndim, function(j) all(!is.na(rho$threshold.values[[j]])) || all(
# (length(which(!is.na(rho$threshold.values[[j]])) >= 1) && (which(!is.na(rho$threshold.values[[j]])) == 1)))))){ #all not NA or only first fixed for each j
# if (attr(terms.formula(rho$formula), "intercept") == 0) stop("formula needs intercept if you want to fix a threshold in model with corGeneral", call. = FALSE)
# m <- "CMORcor2"
# } else if (all(sapply(1:rho$ndim, function(j) all(!is.na(rho$threshold.values[[j]])) || all(is.na(rho$threshold.values[[j]]))))){ #all not NA or all NA for each j
# if (attr(terms.formula(rho$formula), "intercept") == 0) stop("formula needs intercept if you want to fix a threshold in model with corGeneral", call. = FALSE)
# m <- "CMORcor"
# } else stop("Either fix all thresholds in one or more outcome dimensions,
# or consistently in all other outcome dimensions, all first thresholds or none.
# In addition, intercept must be included in the model.", call. = FALSE)
# } else{
# stop("Model is not identifiable. One of the following could help:
# 1. Remove intercept in formula (~ 0 + ...) or
# 2. Fix all first thresholds in threshold.values or
# 3. Fix all thresholds in threshold.values or
# 4. Either fix all thresholds in one or more outcome dimensions,
# or consistently in all other outcome dimensions, all first thresholds or none", call. = FALSE)
# }
# }
# if (rho$model == "CMORcov"){
# if (attr(terms.formula(rho$formula), "intercept") == 0) stop("formula needs intercept if you want to fix a threshold in model CMORcov", call. = FALSE)
# if (all(sapply(1:rho$ndim, function(j){
# (all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1,2))) && all(which(!is.na(rho$threshold.values[[j]])) == c(1,2)) #all first two thresholds are not NA
# ) || ((length(rho$threshold.values[[j]]) == 1) && (which(!is.na(rho$threshold.values[[j]])) == 1))
# }))){
# m <- "CMORcov"
# } else if (all(sapply(1:rho$ndim, function(j){
# (all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1,rho$ntheta[j]))) && all(which(!is.na(rho$threshold.values[[j]])) == c(1,rho$ntheta[j]))#all first and last two thresholds are not NA
# ) || ((length(rho$threshold.values[[j]]) == 1) && (which(!is.na(rho$threshold.values[[j]])) == 1))
# }))){
# m <- "CMORcov2"
# } else if (any(sapply(1:rho$ndim, function(j) all(!is.na(rho$threshold.values[[j]]))))){#all thresholds are fixed in at least one dimensions
# ##CHECK somewhere if all values are fixed in one dimensions + ?warning(sense)
# if (all(sapply(1:rho$ndim, function(j) #
# #(all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1))) && all((which(!is.na(rho$threshold.values[[j]])) == 1))) ||
# (all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1,2))) && all((which(!is.na(rho$threshold.values[[j]])) == c(1,2)))) ||
# (all(!is.na(rho$threshold.values[[j]])))))){
# m <- "CMORcov"
# } else if (all(sapply(1:rho$ndim, function(j)
# #(all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1))) && all((which(!is.na(rho$threshold.values[[j]])) == 1))) ||
# (all(length(which(!is.na(rho$threshold.values[[j]])))==length(c(1,rho$ntheta[j]))) && all((which(!is.na(rho$threshold.values[[j]])) == c(1,rho$ntheta[j])))) ||
# (all(!is.na(rho$threshold.values[[j]])))))){
# m <- "CMORcov2"
# } else stop("Either fix all thresholds in one or more outcome dimensions,
# or consistently in all other outcome dimensions, all first two thresholds or all first and last thresholds", call. = FALSE)
# } else stop("Either fix first two or first and last threshold in threshold.values if you use a model with covGeneral", call. = FALSE)
# }
# }
# }
# m
# }
Try the MultOrd package in your browser
Any scripts or data that you put into this service are public.
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.