Nothing
R/se.R
In MultOrd: Multivariate Ordinal Regression Models
##########################################################
###### PL Standard Errors - no parallel ###############
##########################################################
transf.thresholds.fix2.firstlast.jac <- function(rho, j, gamma_j, i){
recursive.theta <- function(i) {
if (i == 0) 0
else return ((exp(gamma_j[i]) + recursive.theta(i - 1))/(1 + exp(gamma_j[i])))
}
theta <- sapply(1:length(gamma_j), function(i)
recursive.theta(i))
theta[i]
}
transf.thresholds.fix2.first.jac <- function(rho,j,gamma_j,i){
c(0, cumsum(c(1 ,exp(gamma_j))))[i+2]
}
corr.jac.num.fct <- function(rho, nu, i){
# i is the ith correlation parameter
L <- diag(rho$ndim)
angles <- pi * exp(nu)/(1 + exp(nu))
L[lower.tri(L)] <- cos(angles)
S <- matrix(0, nrow = rho$ndim - 1, ncol = rho$ndim - 1)
S[lower.tri(S,diag=T)] <- sin(angles)
S <- apply(cbind(1, rbind(0, S)), 1, cumprod)
L <- L * t(S)
sigma <- tcrossprod(L)
sigma[lower.tri(sigma)][i]
}
PL_se <- function(rho){
par <- rho$optpar
rho$transf.thresholds.jac <- switch(rho$threshold,
fix2firstlast = transf.thresholds.fix2.firstlast.jac,
fix2first = transf.thresholds.fix2.first.jac)
gamma <- par[1:rho$npar.thetas]
if (rho$threshold == "flexible") {
jac <- lapply((1:rho$ndim)[which(rho$npar.theta.opt > 0)], function(j){ #rho$npar.theta
emat <- diag(rho$ntheta[j])
if (ncol(emat) >= 2) {
emat[,1] <- 1
for (k in 2:ncol(emat))
emat[(k:nrow(emat)), k] <-
exp(gamma[(rho$first.ind.theta[j]) + seq_len(rho$ntheta[j]-1)])[k - 1]
}
emat
})
} else {
if (rho$threshold == "fix1first") {
jac <- lapply((1:rho$ndim)[which(rho$npar.theta.opt > 0)], function(j){ #rho$npar.theta
emat <- diag(rho$ntheta[j])
if (ncol(emat) >= 2) {
emat[,1] <- 1
for (k in 2:ncol(emat))
emat[(k:nrow(emat)), k] <-
exp(gamma[(rho$first.ind.theta[j]) + seq_len(rho$npar.theta.opt[j])-1])[k - 1] #rho$npar.theta
}
emat[-1,-1]
})
} else {
jac <- lapply((1:rho$ndim)[which(rho$npar.theta.opt > 0)], function(j){ #rho$npar.theta
gamma_j <- gamma[rho$first.ind.theta[j] + seq_len(rho$npar.theta.opt[j]) - 1] #rho$npar.theta
t(sapply(1:length(gamma_j),
function(i) grad(function(x) rho$transf.thresholds.jac(rho, j, x, i), x=gamma_j)))
})
}
}
## no transform for betas
jac[sum(rho$npar.theta.opt > 0) + seq_len(rho$npar.betas)] <- 1 #rho$npar.theta
## jacobian for spherical transf
nu <- par[rho$npar.thetas + rho$npar.betas + seq_len(rho$npar.cor * rho$ncor.levels)]
if (rho$error.structure$type %in% c("corEqui", "corAR1")){
corr.jac <- list(diag(rho$npar.sigmas))
} else {
corr.jac <- lapply(1:rho$ncor.levels, function(l) t(sapply(1:rho$npar.cor, function(i) grad(function(x) corr.jac.num.fct(rho, x, i),
x=nu[(l - 1) * rho$npar.cor + seq_len(rho$npar.cor)]))))
}
jac[sum(rho$npar.theta.opt > 0) + rho$npar.betas + seq_len(rho$ncor.levels)] <- corr.jac #rho$npar.theta
if (rho$error.structure$type %in% c("covGeneral")) jac[sum(rho$npar.theta.opt > 0) + rho$npar.betas + rho$ncor.levels + seq_len(rho$npar.cor.sd * rho$ncor.levels)] <-
exp(par[rho$npar.thetas + rho$npar.betas + rho$npar.cor * rho$ncor.levels + seq_len(rho$npar.cor.sd * rho$ncor.levels)]) #rho$npar.theta
# exp(par[rho$npar.thetas + rho$npar.betas + rho$npar.cor * rho$ncor.levels + seq_len(rho$npar.cor.sd * rho$ncor.levels)])
J <- as.matrix(bdiag(jac))
J.inv <- solve(J)
## analytic only for vector beta constraints, no fixall and with no PL.lag restriction TODO
if (#all(rho$coef.constraints == matrix(1:rho$ndim, ncol = rho$p, nrow = rho$ndim)) &&
(ncol(unique(rho$coef.constraints, MARGIN = 2)) == 1) &&
# all(rho$threshold.constraints == 1:rho$ndim) &&
(rho$PL.lag == rho$ndim) &&
# all(is.na(rho$coef.values)) &&
all(sapply(rho$threshold.values.fixed,length) != rho$ntheta )) {
cat("Computing variability matrix analytically ... \n")
rho$Vi <- Vi_ana(rho)
} else {
rho$neglog_comp_i <- switch(rho$link,
probit = neglogPL_comp_i.probit,
logit = neglogPL_comp_i.logit)
rho$transf.par_i <- switch(rho$error.structure$type,
corGeneral = transf.par.cor_i,
corEqui = transf.par.cor_i,
covGeneral = transf.par.cov_i,
corAR1 = transf.par.cor_i)
cat("Computing variability matrix numerically ... \n")
Vi <- matrix(0, ncol = length(par), nrow = rho$n)
for (i in 1:rho$n) {
if (i %% 100 == 0) cat('Computed gradient for', i, 'out of', rho$n,'observations\n')
Vi[i, ] <- grad(function(par) rho$neglog_comp_i(par, rho, i), par, method = "Richardson")
}
cat("\n")
# rho$Vi1 <- Vi %*% J.inv
# all.equal(rho$Vi1, rho$Vi)
rho$Vi <- Vi %*% J.inv
}
rho$V <- crossprod(rho$Vi) # original variability matrix
cat("\nComputing Hessian ... \n")
Ht <- hessian(function(par) rho$ObjFun(par, rho), par,
method = "Richardson", method.args=list(eps=1e-6)) # Fisher matrix H(Gamma transf)
rho$V <- rho$n/(rho$n - length(par)) * rho$V ## see maop code, correct for degrees of freedom
rho$H.inv <- J %*% solve(Ht) %*% t(J)
rho$varGamma <- rho$H.inv %*% rho$V %*% rho$H.inv ## inverse godambe
rho$seGamma <- sqrt(diag(rho$varGamma))
rho$claic <- 2 * rho$objective + 2 * sum(diag(rho$V %*% rho$H.inv))
rho$clbic <- 2 * rho$objective + log(rho$n) * sum(diag(rho$V %*% rho$H.inv))
rho
}
###########################################################################
###### neg loglikelihood component for each subject i (gradient) ##############
###########################################################################
transf.par.cor_i <- function(par, rho, i) {
sigmas <- rho$transf.sigmas(par[rho$npar.thetas + rho$npar.betas + seq_len(rho$npar.sigmas)],rho)
#transform thresholds due to monotonicity
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]][i, ] %*% beta[j, ])
theta.lower <- sapply(1:rho$ndim, function(j) c(-rho$inf.value, theta[[j]])[rho$y[i, j]])
theta.upper <- sapply(1:rho$ndim, function(j) c(theta[[j]], rho$inf.value)[rho$y[i, 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_i <- function(par, rho, i) {
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) ])
lev <- match(rho$error.structure$x[i], rho$error.structure$levels)
#transform thresholds due to monotonicity
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]][i, ] %*% beta[j, ])
theta.lower <- sapply(1:rho$ndim, function(j) c(-rho$inf.value, theta[[j]])[rho$y[i, j]])
theta.upper <- sapply(1:rho$ndim, function(j) c(theta[[j]], rho$inf.value)[rho$y[i, j]])
pred.lower <- ((theta.lower - pred.fixed)/rho$sd.y)/exp.par.sd[[lev]]
pred.upper <- ((theta.upper - pred.fixed)/rho$sd.y)/exp.par.sd[[lev]]
sigmas <- lapply(sigmas, cov2cor)
list(U = pred.upper, L = pred.lower,
sigmas = sigmas)
}
neglogPL_comp_i.probit <- function(par, rho, i) {
# transform parameters and get upper and lower bounds
tmp <- rho$transf.par_i(par, rho, i)
U <- tmp$U; L <- tmp$L; sigmas <- tmp$sigmas
q <- which(!is.na(rho$y[i, ]))
if (length(q) == 1){
pr <- pnorm(U[q]) - pnorm(L[q])
logPLi <- rho$weights[i] * log(max(pr, .Machine$double.eps))
} else {
combis <- combn(q, 2)
combis <- combis[,which((combis[2,] - combis[1,]) <= 2 * rho$PL.lag), drop = FALSE]
logPLi <- 0
for (h in 1:ncol(combis)){
if(rho$error.structure$type %in% c("corGeneral", "covGeneral")){ #, "corAR1"
lev <- match(rho$error.structure$x[i], rho$error.structure$levels)
r <- sigmas[[lev]][combis[1, h], combis[2, h]]
} else if(rho$error.structure$type %in% c("corAR1")){
r <- sigmas[[i]][combis[1, h], combis[2, h]]
} else r <- sigmas[i]
pr <- rectbiv.norm.prob(U[combis[, h]], L[combis[, h]], r)
pr[pr < .Machine$double.eps] <- .Machine$double.eps
logPLi <- logPLi + rho$weights[i] * log(max(pr, .Machine$double.eps))
}
}
-logPLi
}
neglogPL_comp_i.logit <- function(par, rho, i) {
# transform parameters and get upper and lower bounds
tmp <- rho$transf.par_i(par, rho, i)
U <- tmp$U; L <- tmp$L; sigmas <- tmp$sigmas
q <- which(!is.na(rho$y[i, ]))
if (length(q) == 1){
pr <- pt(U[q], df = rho$df.t) - pt(L[q], df = rho$df.t)
logPLi <- rho$weights[i] * log(pr)
} else {
combis <- combn(q, 2)
combis <- combis[,which((combis[2,] - combis[1,]) <= 2 * rho$PL.lag), drop = FALSE]
logPLi <- 0
for (h in 1:ncol(combis)){
if(rho$error.structure$type %in% c("corGeneral", "covGeneral")){ #, "corAR1"
lev <- match(rho$error.structure$x[i], rho$error.structure$levels)
r <- sigmas[[lev]][combis[1, h], combis[2, h]]
} else if(rho$error.structure$type %in% c("corAR1")){
r <- sigmas[[i]][combis[1, h], combis[2, h]]
} else r <- sigmas[i]
pr <- biv.nt.prob2(rho$df.t,
lower = L[combis[,h]],
upper = U[combis[,h]],
r = r)
pr[pr < .Machine$double.eps] <- .Machine$double.eps
logPLi <- logPLi + rho$weights[i] * log(max(pr, .Machine$double.eps))
}
}
-logPLi
}
biv.prob <- function(link, U, L, r, df){
## U and L must have dimension 2 and are upper and lower bounds
## r is the correlation parameter
## df are the degrees of freedom for the t-distribution
if (is.null(dim(U)[2]) & is.null(dim(L)[2])) {
dim(U) <- dim(L) <- c(1, length(U))
}
if (length(r) == 1) r <- rep(r, nrow(U))
pr <- switch(link,
probit = rectbiv.norm.prob(U, L, r),
logit = sapply(seq_len(nrow(U)), function(i)
biv.nt.prob2(df = df,
lower = L[i, ],
upper = U[i, ],
r = r[i])))
pr[pr < .Machine$double.eps] <- .Machine$double.eps
pr
}
deriv_biv_norm <- function(A, B, r){
dnorm(A) * pnorm((B - r * A)/sqrt(1 - r^2))
}
deriv_corr_norm <- function(A, B, r){
1/(2 * pi * sqrt(1 - r^2)) *
exp(-(A^2 - 2 * r * A * B + B^2)/(2 * (1 - r^2)))
}
deriv_biv_t <- function(A, B, r, df){
mu_c <- r * A
sigma_c <- sqrt((df + A^2)/(df + 1) * (1 - r^2))
df_c <- df + 1
dt(A, df = df) * pt((B - mu_c)/sigma_c, df = df_c)
}
deriv_corr_t <- function(A, B, r, df){
1/(2 * pi * sqrt(1 - r^2)) *
(1 + (A^2 - 2 * r * A * B + B^2)/(df * (1 - r^2)))^(- df/2)
}
Vi_ana <- function(rho){
pfun <- switch(rho$link,
probit = pnorm,
logit = function(q) pt(q, df = rho$df.t))
dfun <- switch(rho$link,
probit = dnorm,
logit = function(q) dt(q, df = rho$df.t))
deriv_biv_fun <- switch(rho$link,
probit = deriv_biv_norm,
logit = function(A, B, r) deriv_biv_t(A, B, r, df = rho$df.t))
deriv_corr_fun <- switch(rho$link,
probit = deriv_corr_norm,
logit = function(A, B, r) deriv_corr_t(A, B, r, df = rho$df.t))
par <- rho$optpar
tmp <- rho$transf.par(par, rho)
U <- tmp$U; L <- tmp$L; sigmas <- tmp$sigmas;
rho$B2 <- lapply(1:rho$ndim, function(j)
1 * (col(matrix(0, rho$n, rho$ntheta[j] + 1)) ==
c(unclass(rho$y[, j]))))
rho$B1 <- lapply(1:rho$ndim, function(i) as.matrix(rho$B2[[i]][,-(rho$ntheta[i] + 1), drop = FALSE]))
rho$B2 <- lapply(1:rho$ndim, function(i) as.matrix(rho$B2[[i]][,-1, drop = FALSE]))
dbeta <- matrix(0, nrow = rho$n, ncol = rho$npar.betas)
dtheta <- matrix(0, nrow = rho$n, ncol = sum(rho$npar.theta.opt)) #rho$npar.theta
dcorr <- matrix(0, nrow = rho$n, ncol = rho$npar.cor * rho$ncor.levels)
if (rho$error.structure$type == "covGeneral") {
dstddev <- matrix(0, nrow = rho$n, ncol = rho$npar.cor.sd * rho$ncor.levels)
std.dev <- tmp$std.dev;
} else {
dstddev <- NULL
std.dev <- rep(list(rep(1, rho$ndim)), rho$ncor.levels)
}
npar.beta.opt <- apply(rho$ind.coef, 1, function(x) sum(!is.na(x)))
pick.col.beta <- lapply(seq_len(nrow(rho$ind.coef)), function(i)
which(!is.na(rho$ind.coef[i, ])))
for (k in seq_len(length(rho$y.NA.ind))) {
q <- as.numeric(strsplit(names(rho$y.NA.ind[k]), "")[[1]])# turn "101"to 1 0 1)
ind <- rho$y.NA.ind[[k]]
if(! rho$error.structure$type %in% c("corEqui", "corAR1")){
lev <- match(rho$error.structure$x[ind], rho$error.structure$levels)
} else {
lev <- rep(1, length(ind))
}
Uq <- U[ind, , drop = F]
Lq <- L[ind, , drop = F]
std.dev.mat <- do.call("rbind", std.dev[lev])
if (sum(q) == 1){
idq <- which(q == 1)
pr <- pfun(Uq[, idq]) - pfun(Lq[, idq])
pick.col.theta <- switch(rho$threshold,
flexible = seq_len(rho$ntheta[idq]),
fix1first = seq_len(rho$ntheta[idq])[-1],
fix2first = seq_len(rho$ntheta[idq])[-c(1,2)],
fix2firstlast = seq_len(rho$ntheta[idq] - 1)[-1])
if (length(pick.col.theta) != 0) {
dtheta[ind, rho$first.ind.theta[idq] + seq_len(rho$npar.theta[idq]) - 1] <- 1/pr *
(- (dfun(Uq[, idq]) * rho$B1[[idq]][ind, pick.col.theta] -
dfun(Lq[, idq]) * rho$B2[[idq]][ind, pick.col.theta]) *
1/rho$sd.y * 1/std.dev.mat[, idq])
}
dbeta[ind, (rho$first.ind.beta[idq] - rho$npar.thetas) + seq_len(npar.beta.opt[idq]) - 1] <-
1/pr * 1/rho$sd.y * 1/std.dev.mat[, idq] *
(dfun(Uq[, idq]) - dfun(Lq[, idq])) * rho$x[[idq]][ind, pick.col.beta[[idq]]]
if (rho$error.structure$type == "covGeneral"){
for (l in 1:rho$ncor.levels) {
idl <- which(lev == l)
dsd <- rep(0, length(ind))
dsd[idl] <- 1/rho$sd.y * 1/std.dev.mat[idl, idq] *
(dfun(Uq[idl, idq]) * Uq[idl, idq] - dfun(Lq[idl, idq]) * Lq[idl, idq])
dstddev[ind, (l - 1) * rho$npar.cor.sd + idq] <- dsd/pr
}
}
} else {
combis <- combn(which(q == 1), 2, simplify = F)
# combis <- combis[sapply(combis, function(x) (x[2] - x[1]) <= 2 * rho$PL.lag)] # TODO if length(combis) == 0!!!
deriv_theta_rect <- function(k, l, U, L, r){
## derivatives of the rectangle probabilities wrt to thresholds
UU <- deriv_biv_fun(U[, k], U[, l], r)
UL <- deriv_biv_fun(U[, k], L[, l], r)
LU <- deriv_biv_fun(L[, k], U[, l], r)
LL <- deriv_biv_fun(L[, k], L[, l], r)
pick.col.theta <- switch(rho$threshold,
flexible = seq_len(rho$ntheta[k]),
fix1first = seq_len(rho$ntheta[k])[-1],
fix2first = seq_len(rho$ntheta[k])[-c(1,2)],
fix2firstlast = seq_len(rho$ntheta[k] - 1)[-1])
if (length(pick.col.theta) != 0) {
dtheta <- - 1/rho$sd.y * 1/std.dev.mat[, k] *
(UU * rho$B1[[k]][ind, pick.col.theta, drop = F] - UL * rho$B1[[k]][ind, pick.col.theta, drop = F] -
LU * rho$B2[[k]][ind, pick.col.theta, drop = F] + LL * rho$B2[[k]][ind, pick.col.theta, drop = F])
} else {
dtheta <- NULL
}
return(dtheta)
}
deriv_beta_rect <- function(k, l, U, L, r){
## derivatives of the rectangle probabilities wrt to regr.coeff
UU <- deriv_biv_fun(U[, k], U[, l], r)
UL <- deriv_biv_fun(U[, k], L[, l], r)
LU <- deriv_biv_fun(L[, k], U[, l], r)
LL <- deriv_biv_fun(L[, k], L[, l], r)
1/rho$sd.y * 1/std.dev.mat[, k] * (UU - UL - LU + LL) * rho$x[[k]][ind, pick.col.beta[[k]]]
}
deriv_corr_rect <- function(k, l, U, L, r) {
-(deriv_corr_fun(A = U[, k], B = U[, l], r) -
deriv_corr_fun(A = U[, k], B = L[, l], r) -
deriv_corr_fun(A = L[, k], B = U[, l], r) +
deriv_corr_fun(A = L[, k], B = L[, l], r))
}
## gradient wrt threshold parameters
if (rho$npar.thetas != 0) {
dtheta[ind, ] <- do.call("cbind", lapply(unique(rho$threshold.constraints), function(j){
indj <- (1:rho$ndim)[rho$threshold.constraints == j]
combisj <- combis[sapply(combis, function(x) any(indj %in% x))]
## if (length(combisj) != 0) combisj <- combisj[sapply(combisj, function(x) (x[2] - x[1]) <= 2 * rho$PL.lag)]
if (length(combisj) != 0) {
Reduce("+", lapply(combisj, function(comb){
r <- switch(rho$error.structure$type,
corGeneral = sapply(sigmas[lev],'[', comb[1], comb[2]),
covGeneral = sapply(sigmas[lev],'[', comb[1], comb[2]),
corEqui = sigmas[ind, 1],
corAR1 = sapply(sigmas[ind],'[', comb[1], comb[2]))
if (sum(comb %in% indj) == 2){
ind1 <- 1; ind2 <- 2;
} else {
ind1 <- which(comb %in% indj)
ind2 <- which(! comb %in% indj)
}
pr <- biv.prob(link = rho$link,
U = Uq[, comb[c(ind1, ind2)], drop = F],
L = Lq[, comb[c(ind1, ind2)], drop = F],
r = r,
df = rho$df.t)
term1 <- 1/pr * deriv_theta_rect(k = comb[ind1], l = comb[ind2],
U = Uq, L = Lq, r = r)
if (rho$threshold.constraints[comb[1]] == rho$threshold.constraints[comb[2]]) {
term1 <- term1 + 1/pr * deriv_theta_rect(k = comb[ind2], l = comb[ind1],
U = Uq, L = Lq, r = r)
}
if (length(term1) != 0 & is.null(dim(term1)[2])) dim(term1) <- c(1, length(term1))
term1
}))
} else {
matrix(0, nrow = length(ind), ncol = rho$npar.theta.opt[indj[1]])
}}))
}
## gradient wrt regression coefficients
beta.constraints <- c(unique(rho$coef.constraints, MARGIN = 2))
dbeta[ind, ]<- do.call("cbind", lapply(unique(beta.constraints), function(j){
indj <- (1:rho$ndim)[beta.constraints == j]
combisj <- combis[sapply(combis, function(x) any(indj %in% x))]
## if (length(combisj) != 0) combisj <- combisj[sapply(combisj, function(x) (x[2] - x[1]) <= 2 * rho$PL.lag)]
if (length(combisj) != 0) {
Reduce("+", lapply(combisj, function(comb){
r <- switch(rho$error.structure$type,
corGeneral = sapply(sigmas[lev],'[', comb[1], comb[2]),
covGeneral = sapply(sigmas[lev],'[', comb[1], comb[2]),
corEqui = sigmas[ind, 1],
corAR1 = sapply(sigmas[ind],'[', comb[1], comb[2]))
if (sum(comb %in% indj) == 2){
ind1 <- 1; ind2 <- 2;
} else {
ind1 <- which(comb %in% indj)
ind2 <- which(! comb %in% indj)
}
pr <- biv.prob(link = rho$link,
U = Uq[, comb[c(ind1, ind2)], drop = F],
L = Lq[, comb[c(ind1, ind2)], drop = F],
r = r,
df = rho$df.t)
term1 <- 1/pr * deriv_beta_rect(k = comb[ind1], l = comb[ind2],
U = Uq, L = Lq, r = r)
if (beta.constraints[comb[1]] == beta.constraints[comb[2]]) {
term1 <- term1 + 1/pr * deriv_beta_rect(k = comb[ind2], l = comb[ind1],
U = Uq, L = Lq, r = r)
}
if (is.null(dim(term1)[2])) dim(term1) <- c(1, length(term1))
term1
}))
} else {
matrix(0, nrow = length(ind), ncol = npar.beta.opt[indj[1]])
}}))
# dbeta[ind, ] <- do.call("cbind", lapply(1:rho$ndim, function(j){
# combisj <- combis[sapply(combis, function(x) j %in% x)]
# ## if (length(combisj) != 0) combisj <- combisj[sapply(combisj, function(x) (x[2] - x[1]) <= 2 * rho$PL.lag)]
# if (length(combisj) != 0) {
# Reduce("+", lapply(combisj, function(comb){
# r <- switch(rho$error.structure$type,
# corGeneral = sapply(sigmas[lev],'[', comb[1], comb[2]),
# covGeneral = sapply(sigmas[lev],'[', comb[1], comb[2]),
# corEqui = sigmas[ind, 1],
# corAR1 = sapply(sigmas[ind],'[', comb[1], comb[2]))
# pr <- biv.prob(link = rho$link,
# U = Uq[, comb, drop = F],
# L = Lq[, comb, drop = F],
# r = r,
# df = rho$df.t)
# term1 <- 1/pr * deriv_beta_rect(j, comb[comb!=j], Uq, Lq, r)
# if (is.null(dim(term1)[2])) dim(term1) <- c(1, length(term1))
# term1
# }))
# } else {
# matrix(0, nrow = length(ind), ncol = rho$p)
# }
# }))
## gradient wrt correlation parameters
if (rho$error.structure$type == "corAR1") {
dcorr[ind, ] <-
Reduce("+",
lapply(combis, function(comb) {
rpowlag <- sapply(sigmas[ind],'[', comb[1], comb[2])
r <- rpowlag^(1/abs(comb[1] - comb[2]))
pr <- biv.prob(link = rho$link,
U = Uq[, comb, drop = F],
L = Lq[, comb, drop = F],
r = rpowlag,
df = rho$df.t)
dLdr <- abs(comb[1] - comb[2]) * r^(abs(comb[1] - comb[2]) - 1) * 1/pr *
deriv_corr_rect(comb[1], comb[2], Uq, Lq, r = rpowlag)
xbeta <- 0.5 * (log(1 + r) - log(1 - r))
dLdr * exp(2 * xbeta)/(exp(2 * xbeta) + 1)^2 * 4 * rho$error.structure$x[ind, ]
}))
} else {
if (rho$error.structure$type == "corEqui") {
dcorr[ind, ] <- Reduce("+",
lapply(combis, function(comb) {
r <- sigmas[ind, 1]
pr <- biv.prob(link = rho$link,
U = Uq[, comb, drop = F],
L = Lq[, comb, drop = F],
r = r, df = rho$df.t)
dLdr <- 1/pr * deriv_corr_rect(comb[1], comb[2], Uq, Lq, r = r)
xbeta <- 0.5 * (log(1 + r) - log(1 - r))
dLdr * exp(2 * xbeta)/(exp(2 * xbeta) + 1)^2 * 4 * rho$error.structure$x[ind, ]
}))
} else {
poslev <- which(combn(rho$ndim, 2, simplify=F) %in% combis)
for (l in 1:rho$ncor.levels){
dcorr[ind, (l-1) * rho$npar.cor + poslev] <- sapply(combis, function(comb) {
r <- sigmas[[l]][comb[1], comb[2]]
pr <- biv.prob(link = rho$link,
U = Uq[lev == l, comb, drop = F],
L = Lq[lev == l, comb, drop = F],
r = r, df = rho$df.t)
dc <- rep(0, length(ind))
dc[lev == l] <- 1/pr * deriv_corr_rect(k = comb[1], l = comb[2],
U = Uq[lev == l, , drop = F],
L = Lq[lev == l, , drop = F],
r = r)
dc
})
}
}
}
## gradient wrt standard deviation parameters
if (rho$error.structure$type == "covGeneral") {
deriv_stddev_rect <- function(k, l, U, L, r){
UU <- deriv_biv_fun(U[, k], U[, l], r)
UL <- deriv_biv_fun(U[, k], L[, l], r)
LU <- deriv_biv_fun(L[, k], U[, l], r)
LL <- deriv_biv_fun(L[, k], L[, l], r)
(UU * U[, k] - UL * U[, k] - LU * L[, k] + LL * L[, k])
}
for (l in 1:rho$ncor.levels) {
dstddev[ind, (l - 1) * rho$npar.cor.sd + seq_len(rho$npar.cor.sd)] <-
do.call("cbind", lapply(1:rho$ndim, function(j){
combisj <- combis[sapply(combis, function(x) j %in% x)]
## if (length(combisj) != 0) combisj <- combisj[sapply(combisj, function(x) (x[2] - x[1]) <= 2 * rho$PL.lag)]
if (length(combisj) != 0) {
Reduce("+", lapply(combisj, function(comb){
ind1 <- which(comb == j)
ind2 <- which(comb != j)
r <- sigmas[[l]][comb[1], comb[2]]
pr <- biv.prob(link = rho$link,
U = Uq[lev == l, comb, drop = F],
L = Lq[lev == l, comb, drop = F],
r = r, df = rho$df.t)
dsd <- rep(0, length(ind))
dsd[lev == l] <- 1/rho$sd.y * 1/std.dev.mat[lev==l,j] * 1/pr *
deriv_stddev_rect(j, comb[comb!=j],
U = Uq[lev==l,, drop = F],
L = Lq[lev == l, , drop = F], r)
dsd
}))
} else {
matrix(0, nrow = length(ind), ncol = 1)
}
}))
}
}
}
}
Vii <- cbind(dtheta, dbeta, dcorr, dstddev)
return(Vii)
}
# all.equal(Vii, rho$Vi)
#all.equal(c(dcorr), c(rho$Vi[, c(21:23)]))
#all.equal(c(dbeta), c(rho$Vi[, c((rho$npar.theta.opts+1):17)]))
#summary(abs(dcorr - rho$Vi[,21:23]))
#head(rho$Vi[, c(1:rho$npar.theta.opts)])
#all.equal(c(Vi[, c(1:rho$npar.theta.opts)]), c(rho$Vi[, c(1:rho$npar.theta.opts)]))
#summary(abs(Vi[ind, c(1:8)] - rho$Vi[ind, c(1:8)]))
# Hesse_ana_probit <- function(rho){
# par <- rho$optpar
# tmp <- rho$transf.par(par, rho)
# U <- tmp$U; L <- tmp$L; sigmas <- tmp$sigmas
# rho$B2 <- lapply(1:rho$ndim, function(j)
# 1 * (col(matrix(0, rho$n, rho$ntheta[j] + 1)) ==
# c(unclass(rho$y[, j]))))
# rho$B1 <- lapply(1:rho$ndim, function(i) as.matrix(rho$B2[[i]][,-(rho$ntheta[i] + 1), drop = FALSE]))
# rho$B2 <- lapply(1:rho$ndim, function(i) as.matrix(rho$B2[[i]][,-1, drop = FALSE]))
#
# H_comb<-V_comb <- list()
# for (i in 1:length(combis)) {
# comb <- combis[[i]]
# dbeta <- matrix(0, nrow = rho$n, ncol = rho$npar.betas)
# dtheta <- matrix(0, nrow = rho$n, ncol = sum(rho$npar.theta.opt))
# dcorr <- matrix(0, nrow = rho$n, ncol = rho$npar.cor * rho$ncor.levels)
# for (k in 1:length(rho$y.NA.list)){
# q <- as.numeric(strsplit(names(rho$y.NA.list[k]), "")[[1]])# turn "101"to 1 0 1)
# ind <- rho$y.NA.list[[k]]$ind
# #combis <- combn(which(q == 1), 2, simplify = F)
# if(rho$error.structure$type != "corEqui"){
# lev <- match(rho$error.structure$x[ind], rho$error.structure$levels)
# }
# Uq <- U[ind, , drop = F]
# Lq <- L[ind, , drop = F]
# if (sum(q) == 1){
# pr <- pnorm(Uq[, q == 1]) - pnorm(Lq[, q == 1])
# if (rho$model == "CMORcor") dtheta[ind, ] <- - (dnorm(Uq[, q == 1]) * rho$B1[[q == 1]][ind, ] - dnorm(Lq[, q == 1]) * rho$B2[[q == 1]][ind, ])/pr
# if (rho$model == "CMORcor2") dtheta[ind, ] <- - (dnorm(Uq[, q == 1]) * rho$B1[[q == 1]][ind, -1] - dnorm(Lq[, q == 1]) * rho$B2[[q == 1]][ind, -1])/pr
# dbeta[ind, ] <- (dnorm(Uq[, q == 1]) - dnorm(Lq[, q == 1])) * rho$x[[q == 1]]/pr
# } else {
# derivtheta2norm <- function(k,l, U, L,r){
# UU <- derivs2norm(U[, k], U[, l], r)
# UL <- derivs2norm(U[, k], L[, l], r)
# LU <- derivs2norm(L[, k], U[, l], r)
# LL <- derivs2norm(L[, k], L[, l], r)
# if (rho$model == "CMORcor") dth <- -(UU * rho$B1[[k]][ind, ] - UL * rho$B1[[k]][ind, ] - LU * rho$B2[[k]][ind, ] + LL * rho$B2[[k]][ind, ])
# if (rho$model == "CMORcor2") dth <- - (UU * rho$B1[[k]][ind, -1] - UL * rho$B1[[k]][ind, -1] - LU * rho$B2[[k]][ind, -1] + LL * rho$B2[[k]][ind, -1])
# dth
# }
#
# derivbeta2norm <- function(k,l, U, L,r){
# UU <- derivs2norm(U[, k], U[, l], r)
# UL <- derivs2norm(U[, k], L[, l], r)
# LU <- derivs2norm(L[, k], U[, l], r)
# LL <- derivs2norm(L[, k], L[, l], r)
# (UU - UL - LU + LL) * rho$x[[k]][ind, ]
# }
# derivcornorm <- function(k, l, r, U, L) {
# -(corrderivnorm(r, x = U[, k], y = U[, l]) -
# corrderivnorm(r, x = U[, k], y = L[, l]) -
# corrderivnorm(r, x = L[, k], y = U[, l]) +
# corrderivnorm(r, x = L[, k], y = L[, l]))
# }
#
# dtheta[ind,] <- do.call("cbind", lapply(1:rho$ndim, function(j){
# if (j %in% comb){
# r <- sapply(sigmas[lev],'[', comb[1], comb[2])
# pr <- rectbiv.norm.prob(Uq[, c(j, comb[comb!=j])], Lq[, c(j, comb[comb!=j])], r)
# pr[pr < .Machine$double.eps] <- .Machine$double.eps
# derivtheta2norm(j, comb[comb!=j], U=Uq, L=Lq, r=r)/pr
# } else {
# matrix(0, nrow = length(ind), ncol = rho$npar.theta.opt[j])
# } }))
# dbeta[ind, ] <- do.call("cbind", lapply(1:rho$ndim, function(j){
# if (j %in% comb){
# r <- sapply(sigmas[lev],'[', comb[1], comb[2])
# pr <- rectbiv.norm.prob(Uq[, c(j, comb[comb!=j])], Lq[, c(j, comb[comb!=j])], r)
# pr[pr < .Machine$double.eps] <- .Machine$double.eps
# derivbeta2norm(j, comb[comb!=j], Uq, Lq, r)/pr
# } else {
# matrix(0, nrow = length(ind), ncol = rho$p)
# }
# }))
# poslev <- which(combn(rho$npar.cor, 2, simplify=F) %in% list(comb))
# for (l in 1:rho$ncor.levels){
# r <- sigmas[[l]][comb[1], comb[2]]
# pr <- rectbiv.norm.prob(Uq[lev == l, comb], Lq[lev == l, comb], r)
# pr[pr < .Machine$double.eps] <- .Machine$double.eps
# dc <- rep(0, length(ind))
# dc[lev == l] <- derivcornorm(comb[1], comb[2], r = r, Uq[lev == l,], Lq[lev == l, ])/pr
# dcorr[ind, (l-1) * rho$npar.cor + poslev] <- dc
# }
# }
# }
# H_comb[[i]] <- crossprod(cbind(dtheta, dbeta, dcorr))
# }
# i<-1
# Hi_comb <- list()
# aaa <- Reduce("+", H_comb)
# HH_comb <- Reduce("+", Hi_comb)
# head(HH_comb)
# str(H_comb)
# H_comb[[1]][1:3,1:3]
# H_comb[[2]][1:3,1:3]
# H_comb[[3]][1:3,1:3]
# Hnum[1:3,1:3]
# aaa[1:3,1:3]
# H_comb[[1]][4:6,4:6]
# H_comb[[2]][4:6,4:6]
# H_comb[[3]][4:6,4:6]
# Hnum[4:6,4:6]
# aaa[4:6,4:6]
# H_comb[[1]][7:11,7:11]
# H_comb[[2]][7:11,7:11]
# H_comb[[3]][7:11,7:11]
# round(Hnum[7:11,7:11],3)
# aaa[7:11,7:11]
#
# aaa[1:11, 1:11]
# H_comb[[1]][12:15, 12:15]
# H_comb[[2]][12:15, 12:15]
# H_comb[[3]][12:15, 12:15]
# round(Hnum[12:15, 12:15],3)
# aaa[12:15, 12:15]
#
#
# Hnum <- t(J.inv) %*% Ht %*% J.inv
#
# head(aaa)
# all.equal(aaa,HH_comb)
# sum(solve(aaa) == rho$H.inv)
# Hnum <- solve(rho$H.inv)
# head(round(Hnum,3))
# head(round(aaa,3))
# Hnum[1:11, 1:11]
# aaa[1:11, 1:11]
# Hnum[12:20, 12:20]
# aaa[12:20, 12:20]
# Hnum[21:23, 21:23]
# aaa[21:23, 21:23]
# ## no constraints, nada !!!
#
# r <- sapply(sigmas[lev],'[', comb[1], comb[2])
# pr <- rectbiv.norm.prob(Uq[, c(j, comb[comb!=j])], Lq[, c(j, comb[comb!=j])], r)
# dtheta[ind, ] <- derivtheta2norm(j, comb[comb!=j], U=Uq, L=Lq, r=r)/pr
#
#
#
# dbeta[ind, ] <- do.call("cbind", lapply(1:rho$ndim, function(j){
# combisj <- combis[sapply(combis, function(x) j %in% x)]
# if (length(combisj) != 0) {
# Reduce("+", lapply(combisj, function(comb){
# r <- sapply(sigmas[lev],'[', comb[1], comb[2])
# pr <- rectbiv.norm.prob(Uq[, c(j, comb[comb!=j])], Lq[, c(j, comb[comb!=j])], r)
# pr[pr < .Machine$double.eps] <- .Machine$double.eps
# derivbeta2norm(j, comb[comb!=j], Uq, Lq, r)/pr}))
# } else {
# matrix(0, nrow = length(ind), ncol = rho$p)
# }
# }))
# poslev <- which(combn(rho$npar.cor, 2, simplify=F) %in% combis)
# for (l in 1:rho$ncor.levels){
# dcorr[ind, (l-1) * rho$npar.cor + poslev] <- sapply(combis, function(comb) {
# r <- sigmas[[l]][comb[1], comb[2]]
# pr <- rectbiv.norm.prob(Uq[lev == l, comb], Lq[lev == l, comb], r)
# pr[pr < .Machine$double.eps] <- .Machine$double.eps
# dc <- rep(0, length(ind))
# dc[lev == l] <- derivcornorm(comb[1], comb[2], r = r, Uq[lev == l,], Lq[lev == l, ])/pr
# dc
# })
# }
# }
# }
# Vii <- cbind(dtheta, dbeta, dcorr)
# return(Vii)
# }
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.
##########################################################
###### PL Standard Errors - no parallel ###############
##########################################################
transf.thresholds.fix2.firstlast.jac <- function(rho, j, gamma_j, i){
recursive.theta <- function(i) {
if (i == 0) 0
else return ((exp(gamma_j[i]) + recursive.theta(i - 1))/(1 + exp(gamma_j[i])))
}
theta <- sapply(1:length(gamma_j), function(i)
recursive.theta(i))
theta[i]
}
transf.thresholds.fix2.first.jac <- function(rho,j,gamma_j,i){
c(0, cumsum(c(1 ,exp(gamma_j))))[i+2]
}
corr.jac.num.fct <- function(rho, nu, i){
# i is the ith correlation parameter
L <- diag(rho$ndim)
angles <- pi * exp(nu)/(1 + exp(nu))
L[lower.tri(L)] <- cos(angles)
S <- matrix(0, nrow = rho$ndim - 1, ncol = rho$ndim - 1)
S[lower.tri(S,diag=T)] <- sin(angles)
S <- apply(cbind(1, rbind(0, S)), 1, cumprod)
L <- L * t(S)
sigma <- tcrossprod(L)
sigma[lower.tri(sigma)][i]
}
PL_se <- function(rho){
par <- rho$optpar
rho$transf.thresholds.jac <- switch(rho$threshold,
fix2firstlast = transf.thresholds.fix2.firstlast.jac,
fix2first = transf.thresholds.fix2.first.jac)
gamma <- par[1:rho$npar.thetas]
if (rho$threshold == "flexible") {
jac <- lapply((1:rho$ndim)[which(rho$npar.theta.opt > 0)], function(j){ #rho$npar.theta
emat <- diag(rho$ntheta[j])
if (ncol(emat) >= 2) {
emat[,1] <- 1
for (k in 2:ncol(emat))
emat[(k:nrow(emat)), k] <-
exp(gamma[(rho$first.ind.theta[j]) + seq_len(rho$ntheta[j]-1)])[k - 1]
}
emat
})
} else {
if (rho$threshold == "fix1first") {
jac <- lapply((1:rho$ndim)[which(rho$npar.theta.opt > 0)], function(j){ #rho$npar.theta
emat <- diag(rho$ntheta[j])
if (ncol(emat) >= 2) {
emat[,1] <- 1
for (k in 2:ncol(emat))
emat[(k:nrow(emat)), k] <-
exp(gamma[(rho$first.ind.theta[j]) + seq_len(rho$npar.theta.opt[j])-1])[k - 1] #rho$npar.theta
}
emat[-1,-1]
})
} else {
jac <- lapply((1:rho$ndim)[which(rho$npar.theta.opt > 0)], function(j){ #rho$npar.theta
gamma_j <- gamma[rho$first.ind.theta[j] + seq_len(rho$npar.theta.opt[j]) - 1] #rho$npar.theta
t(sapply(1:length(gamma_j),
function(i) grad(function(x) rho$transf.thresholds.jac(rho, j, x, i), x=gamma_j)))
})
}
}
## no transform for betas
jac[sum(rho$npar.theta.opt > 0) + seq_len(rho$npar.betas)] <- 1 #rho$npar.theta
## jacobian for spherical transf
nu <- par[rho$npar.thetas + rho$npar.betas + seq_len(rho$npar.cor * rho$ncor.levels)]
if (rho$error.structure$type %in% c("corEqui", "corAR1")){
corr.jac <- list(diag(rho$npar.sigmas))
} else {
corr.jac <- lapply(1:rho$ncor.levels, function(l) t(sapply(1:rho$npar.cor, function(i) grad(function(x) corr.jac.num.fct(rho, x, i),
x=nu[(l - 1) * rho$npar.cor + seq_len(rho$npar.cor)]))))
}
jac[sum(rho$npar.theta.opt > 0) + rho$npar.betas + seq_len(rho$ncor.levels)] <- corr.jac #rho$npar.theta
if (rho$error.structure$type %in% c("covGeneral")) jac[sum(rho$npar.theta.opt > 0) + rho$npar.betas + rho$ncor.levels + seq_len(rho$npar.cor.sd * rho$ncor.levels)] <-
exp(par[rho$npar.thetas + rho$npar.betas + rho$npar.cor * rho$ncor.levels + seq_len(rho$npar.cor.sd * rho$ncor.levels)]) #rho$npar.theta
# exp(par[rho$npar.thetas + rho$npar.betas + rho$npar.cor * rho$ncor.levels + seq_len(rho$npar.cor.sd * rho$ncor.levels)])
J <- as.matrix(bdiag(jac))
J.inv <- solve(J)
## analytic only for vector beta constraints, no fixall and with no PL.lag restriction TODO
if (#all(rho$coef.constraints == matrix(1:rho$ndim, ncol = rho$p, nrow = rho$ndim)) &&
(ncol(unique(rho$coef.constraints, MARGIN = 2)) == 1) &&
# all(rho$threshold.constraints == 1:rho$ndim) &&
(rho$PL.lag == rho$ndim) &&
# all(is.na(rho$coef.values)) &&
all(sapply(rho$threshold.values.fixed,length) != rho$ntheta )) {
cat("Computing variability matrix analytically ... \n")
rho$Vi <- Vi_ana(rho)
} else {
rho$neglog_comp_i <- switch(rho$link,
probit = neglogPL_comp_i.probit,
logit = neglogPL_comp_i.logit)
rho$transf.par_i <- switch(rho$error.structure$type,
corGeneral = transf.par.cor_i,
corEqui = transf.par.cor_i,
covGeneral = transf.par.cov_i,
corAR1 = transf.par.cor_i)
cat("Computing variability matrix numerically ... \n")
Vi <- matrix(0, ncol = length(par), nrow = rho$n)
for (i in 1:rho$n) {
if (i %% 100 == 0) cat('Computed gradient for', i, 'out of', rho$n,'observations\n')
Vi[i, ] <- grad(function(par) rho$neglog_comp_i(par, rho, i), par, method = "Richardson")
}
cat("\n")
# rho$Vi1 <- Vi %*% J.inv
# all.equal(rho$Vi1, rho$Vi)
rho$Vi <- Vi %*% J.inv
}
rho$V <- crossprod(rho$Vi) # original variability matrix
cat("\nComputing Hessian ... \n")
Ht <- hessian(function(par) rho$ObjFun(par, rho), par,
method = "Richardson", method.args=list(eps=1e-6)) # Fisher matrix H(Gamma transf)
rho$V <- rho$n/(rho$n - length(par)) * rho$V ## see maop code, correct for degrees of freedom
rho$H.inv <- J %*% solve(Ht) %*% t(J)
rho$varGamma <- rho$H.inv %*% rho$V %*% rho$H.inv ## inverse godambe
rho$seGamma <- sqrt(diag(rho$varGamma))
rho$claic <- 2 * rho$objective + 2 * sum(diag(rho$V %*% rho$H.inv))
rho$clbic <- 2 * rho$objective + log(rho$n) * sum(diag(rho$V %*% rho$H.inv))
rho
}
###########################################################################
###### neg loglikelihood component for each subject i (gradient) ##############
###########################################################################
transf.par.cor_i <- function(par, rho, i) {
sigmas <- rho$transf.sigmas(par[rho$npar.thetas + rho$npar.betas + seq_len(rho$npar.sigmas)],rho)
#transform thresholds due to monotonicity
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]][i, ] %*% beta[j, ])
theta.lower <- sapply(1:rho$ndim, function(j) c(-rho$inf.value, theta[[j]])[rho$y[i, j]])
theta.upper <- sapply(1:rho$ndim, function(j) c(theta[[j]], rho$inf.value)[rho$y[i, 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_i <- function(par, rho, i) {
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) ])
lev <- match(rho$error.structure$x[i], rho$error.structure$levels)
#transform thresholds due to monotonicity
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]][i, ] %*% beta[j, ])
theta.lower <- sapply(1:rho$ndim, function(j) c(-rho$inf.value, theta[[j]])[rho$y[i, j]])
theta.upper <- sapply(1:rho$ndim, function(j) c(theta[[j]], rho$inf.value)[rho$y[i, j]])
pred.lower <- ((theta.lower - pred.fixed)/rho$sd.y)/exp.par.sd[[lev]]
pred.upper <- ((theta.upper - pred.fixed)/rho$sd.y)/exp.par.sd[[lev]]
sigmas <- lapply(sigmas, cov2cor)
list(U = pred.upper, L = pred.lower,
sigmas = sigmas)
}
neglogPL_comp_i.probit <- function(par, rho, i) {
# transform parameters and get upper and lower bounds
tmp <- rho$transf.par_i(par, rho, i)
U <- tmp$U; L <- tmp$L; sigmas <- tmp$sigmas
q <- which(!is.na(rho$y[i, ]))
if (length(q) == 1){
pr <- pnorm(U[q]) - pnorm(L[q])
logPLi <- rho$weights[i] * log(max(pr, .Machine$double.eps))
} else {
combis <- combn(q, 2)
combis <- combis[,which((combis[2,] - combis[1,]) <= 2 * rho$PL.lag), drop = FALSE]
logPLi <- 0
for (h in 1:ncol(combis)){
if(rho$error.structure$type %in% c("corGeneral", "covGeneral")){ #, "corAR1"
lev <- match(rho$error.structure$x[i], rho$error.structure$levels)
r <- sigmas[[lev]][combis[1, h], combis[2, h]]
} else if(rho$error.structure$type %in% c("corAR1")){
r <- sigmas[[i]][combis[1, h], combis[2, h]]
} else r <- sigmas[i]
pr <- rectbiv.norm.prob(U[combis[, h]], L[combis[, h]], r)
pr[pr < .Machine$double.eps] <- .Machine$double.eps
logPLi <- logPLi + rho$weights[i] * log(max(pr, .Machine$double.eps))
}
}
-logPLi
}
neglogPL_comp_i.logit <- function(par, rho, i) {
# transform parameters and get upper and lower bounds
tmp <- rho$transf.par_i(par, rho, i)
U <- tmp$U; L <- tmp$L; sigmas <- tmp$sigmas
q <- which(!is.na(rho$y[i, ]))
if (length(q) == 1){
pr <- pt(U[q], df = rho$df.t) - pt(L[q], df = rho$df.t)
logPLi <- rho$weights[i] * log(pr)
} else {
combis <- combn(q, 2)
combis <- combis[,which((combis[2,] - combis[1,]) <= 2 * rho$PL.lag), drop = FALSE]
logPLi <- 0
for (h in 1:ncol(combis)){
if(rho$error.structure$type %in% c("corGeneral", "covGeneral")){ #, "corAR1"
lev <- match(rho$error.structure$x[i], rho$error.structure$levels)
r <- sigmas[[lev]][combis[1, h], combis[2, h]]
} else if(rho$error.structure$type %in% c("corAR1")){
r <- sigmas[[i]][combis[1, h], combis[2, h]]
} else r <- sigmas[i]
pr <- biv.nt.prob2(rho$df.t,
lower = L[combis[,h]],
upper = U[combis[,h]],
r = r)
pr[pr < .Machine$double.eps] <- .Machine$double.eps
logPLi <- logPLi + rho$weights[i] * log(max(pr, .Machine$double.eps))
}
}
-logPLi
}
biv.prob <- function(link, U, L, r, df){
## U and L must have dimension 2 and are upper and lower bounds
## r is the correlation parameter
## df are the degrees of freedom for the t-distribution
if (is.null(dim(U)[2]) & is.null(dim(L)[2])) {
dim(U) <- dim(L) <- c(1, length(U))
}
if (length(r) == 1) r <- rep(r, nrow(U))
pr <- switch(link,
probit = rectbiv.norm.prob(U, L, r),
logit = sapply(seq_len(nrow(U)), function(i)
biv.nt.prob2(df = df,
lower = L[i, ],
upper = U[i, ],
r = r[i])))
pr[pr < .Machine$double.eps] <- .Machine$double.eps
pr
}
deriv_biv_norm <- function(A, B, r){
dnorm(A) * pnorm((B - r * A)/sqrt(1 - r^2))
}
deriv_corr_norm <- function(A, B, r){
1/(2 * pi * sqrt(1 - r^2)) *
exp(-(A^2 - 2 * r * A * B + B^2)/(2 * (1 - r^2)))
}
deriv_biv_t <- function(A, B, r, df){
mu_c <- r * A
sigma_c <- sqrt((df + A^2)/(df + 1) * (1 - r^2))
df_c <- df + 1
dt(A, df = df) * pt((B - mu_c)/sigma_c, df = df_c)
}
deriv_corr_t <- function(A, B, r, df){
1/(2 * pi * sqrt(1 - r^2)) *
(1 + (A^2 - 2 * r * A * B + B^2)/(df * (1 - r^2)))^(- df/2)
}
Vi_ana <- function(rho){
pfun <- switch(rho$link,
probit = pnorm,
logit = function(q) pt(q, df = rho$df.t))
dfun <- switch(rho$link,
probit = dnorm,
logit = function(q) dt(q, df = rho$df.t))
deriv_biv_fun <- switch(rho$link,
probit = deriv_biv_norm,
logit = function(A, B, r) deriv_biv_t(A, B, r, df = rho$df.t))
deriv_corr_fun <- switch(rho$link,
probit = deriv_corr_norm,
logit = function(A, B, r) deriv_corr_t(A, B, r, df = rho$df.t))
par <- rho$optpar
tmp <- rho$transf.par(par, rho)
U <- tmp$U; L <- tmp$L; sigmas <- tmp$sigmas;
rho$B2 <- lapply(1:rho$ndim, function(j)
1 * (col(matrix(0, rho$n, rho$ntheta[j] + 1)) ==
c(unclass(rho$y[, j]))))
rho$B1 <- lapply(1:rho$ndim, function(i) as.matrix(rho$B2[[i]][,-(rho$ntheta[i] + 1), drop = FALSE]))
rho$B2 <- lapply(1:rho$ndim, function(i) as.matrix(rho$B2[[i]][,-1, drop = FALSE]))
dbeta <- matrix(0, nrow = rho$n, ncol = rho$npar.betas)
dtheta <- matrix(0, nrow = rho$n, ncol = sum(rho$npar.theta.opt)) #rho$npar.theta
dcorr <- matrix(0, nrow = rho$n, ncol = rho$npar.cor * rho$ncor.levels)
if (rho$error.structure$type == "covGeneral") {
dstddev <- matrix(0, nrow = rho$n, ncol = rho$npar.cor.sd * rho$ncor.levels)
std.dev <- tmp$std.dev;
} else {
dstddev <- NULL
std.dev <- rep(list(rep(1, rho$ndim)), rho$ncor.levels)
}
npar.beta.opt <- apply(rho$ind.coef, 1, function(x) sum(!is.na(x)))
pick.col.beta <- lapply(seq_len(nrow(rho$ind.coef)), function(i)
which(!is.na(rho$ind.coef[i, ])))
for (k in seq_len(length(rho$y.NA.ind))) {
q <- as.numeric(strsplit(names(rho$y.NA.ind[k]), "")[[1]])# turn "101"to 1 0 1)
ind <- rho$y.NA.ind[[k]]
if(! rho$error.structure$type %in% c("corEqui", "corAR1")){
lev <- match(rho$error.structure$x[ind], rho$error.structure$levels)
} else {
lev <- rep(1, length(ind))
}
Uq <- U[ind, , drop = F]
Lq <- L[ind, , drop = F]
std.dev.mat <- do.call("rbind", std.dev[lev])
if (sum(q) == 1){
idq <- which(q == 1)
pr <- pfun(Uq[, idq]) - pfun(Lq[, idq])
pick.col.theta <- switch(rho$threshold,
flexible = seq_len(rho$ntheta[idq]),
fix1first = seq_len(rho$ntheta[idq])[-1],
fix2first = seq_len(rho$ntheta[idq])[-c(1,2)],
fix2firstlast = seq_len(rho$ntheta[idq] - 1)[-1])
if (length(pick.col.theta) != 0) {
dtheta[ind, rho$first.ind.theta[idq] + seq_len(rho$npar.theta[idq]) - 1] <- 1/pr *
(- (dfun(Uq[, idq]) * rho$B1[[idq]][ind, pick.col.theta] -
dfun(Lq[, idq]) * rho$B2[[idq]][ind, pick.col.theta]) *
1/rho$sd.y * 1/std.dev.mat[, idq])
}
dbeta[ind, (rho$first.ind.beta[idq] - rho$npar.thetas) + seq_len(npar.beta.opt[idq]) - 1] <-
1/pr * 1/rho$sd.y * 1/std.dev.mat[, idq] *
(dfun(Uq[, idq]) - dfun(Lq[, idq])) * rho$x[[idq]][ind, pick.col.beta[[idq]]]
if (rho$error.structure$type == "covGeneral"){
for (l in 1:rho$ncor.levels) {
idl <- which(lev == l)
dsd <- rep(0, length(ind))
dsd[idl] <- 1/rho$sd.y * 1/std.dev.mat[idl, idq] *
(dfun(Uq[idl, idq]) * Uq[idl, idq] - dfun(Lq[idl, idq]) * Lq[idl, idq])
dstddev[ind, (l - 1) * rho$npar.cor.sd + idq] <- dsd/pr
}
}
} else {
combis <- combn(which(q == 1), 2, simplify = F)
# combis <- combis[sapply(combis, function(x) (x[2] - x[1]) <= 2 * rho$PL.lag)] # TODO if length(combis) == 0!!!
deriv_theta_rect <- function(k, l, U, L, r){
## derivatives of the rectangle probabilities wrt to thresholds
UU <- deriv_biv_fun(U[, k], U[, l], r)
UL <- deriv_biv_fun(U[, k], L[, l], r)
LU <- deriv_biv_fun(L[, k], U[, l], r)
LL <- deriv_biv_fun(L[, k], L[, l], r)
pick.col.theta <- switch(rho$threshold,
flexible = seq_len(rho$ntheta[k]),
fix1first = seq_len(rho$ntheta[k])[-1],
fix2first = seq_len(rho$ntheta[k])[-c(1,2)],
fix2firstlast = seq_len(rho$ntheta[k] - 1)[-1])
if (length(pick.col.theta) != 0) {
dtheta <- - 1/rho$sd.y * 1/std.dev.mat[, k] *
(UU * rho$B1[[k]][ind, pick.col.theta, drop = F] - UL * rho$B1[[k]][ind, pick.col.theta, drop = F] -
LU * rho$B2[[k]][ind, pick.col.theta, drop = F] + LL * rho$B2[[k]][ind, pick.col.theta, drop = F])
} else {
dtheta <- NULL
}
return(dtheta)
}
deriv_beta_rect <- function(k, l, U, L, r){
## derivatives of the rectangle probabilities wrt to regr.coeff
UU <- deriv_biv_fun(U[, k], U[, l], r)
UL <- deriv_biv_fun(U[, k], L[, l], r)
LU <- deriv_biv_fun(L[, k], U[, l], r)
LL <- deriv_biv_fun(L[, k], L[, l], r)
1/rho$sd.y * 1/std.dev.mat[, k] * (UU - UL - LU + LL) * rho$x[[k]][ind, pick.col.beta[[k]]]
}
deriv_corr_rect <- function(k, l, U, L, r) {
-(deriv_corr_fun(A = U[, k], B = U[, l], r) -
deriv_corr_fun(A = U[, k], B = L[, l], r) -
deriv_corr_fun(A = L[, k], B = U[, l], r) +
deriv_corr_fun(A = L[, k], B = L[, l], r))
}
## gradient wrt threshold parameters
if (rho$npar.thetas != 0) {
dtheta[ind, ] <- do.call("cbind", lapply(unique(rho$threshold.constraints), function(j){
indj <- (1:rho$ndim)[rho$threshold.constraints == j]
combisj <- combis[sapply(combis, function(x) any(indj %in% x))]
## if (length(combisj) != 0) combisj <- combisj[sapply(combisj, function(x) (x[2] - x[1]) <= 2 * rho$PL.lag)]
if (length(combisj) != 0) {
Reduce("+", lapply(combisj, function(comb){
r <- switch(rho$error.structure$type,
corGeneral = sapply(sigmas[lev],'[', comb[1], comb[2]),
covGeneral = sapply(sigmas[lev],'[', comb[1], comb[2]),
corEqui = sigmas[ind, 1],
corAR1 = sapply(sigmas[ind],'[', comb[1], comb[2]))
if (sum(comb %in% indj) == 2){
ind1 <- 1; ind2 <- 2;
} else {
ind1 <- which(comb %in% indj)
ind2 <- which(! comb %in% indj)
}
pr <- biv.prob(link = rho$link,
U = Uq[, comb[c(ind1, ind2)], drop = F],
L = Lq[, comb[c(ind1, ind2)], drop = F],
r = r,
df = rho$df.t)
term1 <- 1/pr * deriv_theta_rect(k = comb[ind1], l = comb[ind2],
U = Uq, L = Lq, r = r)
if (rho$threshold.constraints[comb[1]] == rho$threshold.constraints[comb[2]]) {
term1 <- term1 + 1/pr * deriv_theta_rect(k = comb[ind2], l = comb[ind1],
U = Uq, L = Lq, r = r)
}
if (length(term1) != 0 & is.null(dim(term1)[2])) dim(term1) <- c(1, length(term1))
term1
}))
} else {
matrix(0, nrow = length(ind), ncol = rho$npar.theta.opt[indj[1]])
}}))
}
## gradient wrt regression coefficients
beta.constraints <- c(unique(rho$coef.constraints, MARGIN = 2))
dbeta[ind, ]<- do.call("cbind", lapply(unique(beta.constraints), function(j){
indj <- (1:rho$ndim)[beta.constraints == j]
combisj <- combis[sapply(combis, function(x) any(indj %in% x))]
## if (length(combisj) != 0) combisj <- combisj[sapply(combisj, function(x) (x[2] - x[1]) <= 2 * rho$PL.lag)]
if (length(combisj) != 0) {
Reduce("+", lapply(combisj, function(comb){
r <- switch(rho$error.structure$type,
corGeneral = sapply(sigmas[lev],'[', comb[1], comb[2]),
covGeneral = sapply(sigmas[lev],'[', comb[1], comb[2]),
corEqui = sigmas[ind, 1],
corAR1 = sapply(sigmas[ind],'[', comb[1], comb[2]))
if (sum(comb %in% indj) == 2){
ind1 <- 1; ind2 <- 2;
} else {
ind1 <- which(comb %in% indj)
ind2 <- which(! comb %in% indj)
}
pr <- biv.prob(link = rho$link,
U = Uq[, comb[c(ind1, ind2)], drop = F],
L = Lq[, comb[c(ind1, ind2)], drop = F],
r = r,
df = rho$df.t)
term1 <- 1/pr * deriv_beta_rect(k = comb[ind1], l = comb[ind2],
U = Uq, L = Lq, r = r)
if (beta.constraints[comb[1]] == beta.constraints[comb[2]]) {
term1 <- term1 + 1/pr * deriv_beta_rect(k = comb[ind2], l = comb[ind1],
U = Uq, L = Lq, r = r)
}
if (is.null(dim(term1)[2])) dim(term1) <- c(1, length(term1))
term1
}))
} else {
matrix(0, nrow = length(ind), ncol = npar.beta.opt[indj[1]])
}}))
# dbeta[ind, ] <- do.call("cbind", lapply(1:rho$ndim, function(j){
# combisj <- combis[sapply(combis, function(x) j %in% x)]
# ## if (length(combisj) != 0) combisj <- combisj[sapply(combisj, function(x) (x[2] - x[1]) <= 2 * rho$PL.lag)]
# if (length(combisj) != 0) {
# Reduce("+", lapply(combisj, function(comb){
# r <- switch(rho$error.structure$type,
# corGeneral = sapply(sigmas[lev],'[', comb[1], comb[2]),
# covGeneral = sapply(sigmas[lev],'[', comb[1], comb[2]),
# corEqui = sigmas[ind, 1],
# corAR1 = sapply(sigmas[ind],'[', comb[1], comb[2]))
# pr <- biv.prob(link = rho$link,
# U = Uq[, comb, drop = F],
# L = Lq[, comb, drop = F],
# r = r,
# df = rho$df.t)
# term1 <- 1/pr * deriv_beta_rect(j, comb[comb!=j], Uq, Lq, r)
# if (is.null(dim(term1)[2])) dim(term1) <- c(1, length(term1))
# term1
# }))
# } else {
# matrix(0, nrow = length(ind), ncol = rho$p)
# }
# }))
## gradient wrt correlation parameters
if (rho$error.structure$type == "corAR1") {
dcorr[ind, ] <-
Reduce("+",
lapply(combis, function(comb) {
rpowlag <- sapply(sigmas[ind],'[', comb[1], comb[2])
r <- rpowlag^(1/abs(comb[1] - comb[2]))
pr <- biv.prob(link = rho$link,
U = Uq[, comb, drop = F],
L = Lq[, comb, drop = F],
r = rpowlag,
df = rho$df.t)
dLdr <- abs(comb[1] - comb[2]) * r^(abs(comb[1] - comb[2]) - 1) * 1/pr *
deriv_corr_rect(comb[1], comb[2], Uq, Lq, r = rpowlag)
xbeta <- 0.5 * (log(1 + r) - log(1 - r))
dLdr * exp(2 * xbeta)/(exp(2 * xbeta) + 1)^2 * 4 * rho$error.structure$x[ind, ]
}))
} else {
if (rho$error.structure$type == "corEqui") {
dcorr[ind, ] <- Reduce("+",
lapply(combis, function(comb) {
r <- sigmas[ind, 1]
pr <- biv.prob(link = rho$link,
U = Uq[, comb, drop = F],
L = Lq[, comb, drop = F],
r = r, df = rho$df.t)
dLdr <- 1/pr * deriv_corr_rect(comb[1], comb[2], Uq, Lq, r = r)
xbeta <- 0.5 * (log(1 + r) - log(1 - r))
dLdr * exp(2 * xbeta)/(exp(2 * xbeta) + 1)^2 * 4 * rho$error.structure$x[ind, ]
}))
} else {
poslev <- which(combn(rho$ndim, 2, simplify=F) %in% combis)
for (l in 1:rho$ncor.levels){
dcorr[ind, (l-1) * rho$npar.cor + poslev] <- sapply(combis, function(comb) {
r <- sigmas[[l]][comb[1], comb[2]]
pr <- biv.prob(link = rho$link,
U = Uq[lev == l, comb, drop = F],
L = Lq[lev == l, comb, drop = F],
r = r, df = rho$df.t)
dc <- rep(0, length(ind))
dc[lev == l] <- 1/pr * deriv_corr_rect(k = comb[1], l = comb[2],
U = Uq[lev == l, , drop = F],
L = Lq[lev == l, , drop = F],
r = r)
dc
})
}
}
}
## gradient wrt standard deviation parameters
if (rho$error.structure$type == "covGeneral") {
deriv_stddev_rect <- function(k, l, U, L, r){
UU <- deriv_biv_fun(U[, k], U[, l], r)
UL <- deriv_biv_fun(U[, k], L[, l], r)
LU <- deriv_biv_fun(L[, k], U[, l], r)
LL <- deriv_biv_fun(L[, k], L[, l], r)
(UU * U[, k] - UL * U[, k] - LU * L[, k] + LL * L[, k])
}
for (l in 1:rho$ncor.levels) {
dstddev[ind, (l - 1) * rho$npar.cor.sd + seq_len(rho$npar.cor.sd)] <-
do.call("cbind", lapply(1:rho$ndim, function(j){
combisj <- combis[sapply(combis, function(x) j %in% x)]
## if (length(combisj) != 0) combisj <- combisj[sapply(combisj, function(x) (x[2] - x[1]) <= 2 * rho$PL.lag)]
if (length(combisj) != 0) {
Reduce("+", lapply(combisj, function(comb){
ind1 <- which(comb == j)
ind2 <- which(comb != j)
r <- sigmas[[l]][comb[1], comb[2]]
pr <- biv.prob(link = rho$link,
U = Uq[lev == l, comb, drop = F],
L = Lq[lev == l, comb, drop = F],
r = r, df = rho$df.t)
dsd <- rep(0, length(ind))
dsd[lev == l] <- 1/rho$sd.y * 1/std.dev.mat[lev==l,j] * 1/pr *
deriv_stddev_rect(j, comb[comb!=j],
U = Uq[lev==l,, drop = F],
L = Lq[lev == l, , drop = F], r)
dsd
}))
} else {
matrix(0, nrow = length(ind), ncol = 1)
}
}))
}
}
}
}
Vii <- cbind(dtheta, dbeta, dcorr, dstddev)
return(Vii)
}
# all.equal(Vii, rho$Vi)
#all.equal(c(dcorr), c(rho$Vi[, c(21:23)]))
#all.equal(c(dbeta), c(rho$Vi[, c((rho$npar.theta.opts+1):17)]))
#summary(abs(dcorr - rho$Vi[,21:23]))
#head(rho$Vi[, c(1:rho$npar.theta.opts)])
#all.equal(c(Vi[, c(1:rho$npar.theta.opts)]), c(rho$Vi[, c(1:rho$npar.theta.opts)]))
#summary(abs(Vi[ind, c(1:8)] - rho$Vi[ind, c(1:8)]))
# Hesse_ana_probit <- function(rho){
# par <- rho$optpar
# tmp <- rho$transf.par(par, rho)
# U <- tmp$U; L <- tmp$L; sigmas <- tmp$sigmas
# rho$B2 <- lapply(1:rho$ndim, function(j)
# 1 * (col(matrix(0, rho$n, rho$ntheta[j] + 1)) ==
# c(unclass(rho$y[, j]))))
# rho$B1 <- lapply(1:rho$ndim, function(i) as.matrix(rho$B2[[i]][,-(rho$ntheta[i] + 1), drop = FALSE]))
# rho$B2 <- lapply(1:rho$ndim, function(i) as.matrix(rho$B2[[i]][,-1, drop = FALSE]))
#
# H_comb<-V_comb <- list()
# for (i in 1:length(combis)) {
# comb <- combis[[i]]
# dbeta <- matrix(0, nrow = rho$n, ncol = rho$npar.betas)
# dtheta <- matrix(0, nrow = rho$n, ncol = sum(rho$npar.theta.opt))
# dcorr <- matrix(0, nrow = rho$n, ncol = rho$npar.cor * rho$ncor.levels)
# for (k in 1:length(rho$y.NA.list)){
# q <- as.numeric(strsplit(names(rho$y.NA.list[k]), "")[[1]])# turn "101"to 1 0 1)
# ind <- rho$y.NA.list[[k]]$ind
# #combis <- combn(which(q == 1), 2, simplify = F)
# if(rho$error.structure$type != "corEqui"){
# lev <- match(rho$error.structure$x[ind], rho$error.structure$levels)
# }
# Uq <- U[ind, , drop = F]
# Lq <- L[ind, , drop = F]
# if (sum(q) == 1){
# pr <- pnorm(Uq[, q == 1]) - pnorm(Lq[, q == 1])
# if (rho$model == "CMORcor") dtheta[ind, ] <- - (dnorm(Uq[, q == 1]) * rho$B1[[q == 1]][ind, ] - dnorm(Lq[, q == 1]) * rho$B2[[q == 1]][ind, ])/pr
# if (rho$model == "CMORcor2") dtheta[ind, ] <- - (dnorm(Uq[, q == 1]) * rho$B1[[q == 1]][ind, -1] - dnorm(Lq[, q == 1]) * rho$B2[[q == 1]][ind, -1])/pr
# dbeta[ind, ] <- (dnorm(Uq[, q == 1]) - dnorm(Lq[, q == 1])) * rho$x[[q == 1]]/pr
# } else {
# derivtheta2norm <- function(k,l, U, L,r){
# UU <- derivs2norm(U[, k], U[, l], r)
# UL <- derivs2norm(U[, k], L[, l], r)
# LU <- derivs2norm(L[, k], U[, l], r)
# LL <- derivs2norm(L[, k], L[, l], r)
# if (rho$model == "CMORcor") dth <- -(UU * rho$B1[[k]][ind, ] - UL * rho$B1[[k]][ind, ] - LU * rho$B2[[k]][ind, ] + LL * rho$B2[[k]][ind, ])
# if (rho$model == "CMORcor2") dth <- - (UU * rho$B1[[k]][ind, -1] - UL * rho$B1[[k]][ind, -1] - LU * rho$B2[[k]][ind, -1] + LL * rho$B2[[k]][ind, -1])
# dth
# }
#
# derivbeta2norm <- function(k,l, U, L,r){
# UU <- derivs2norm(U[, k], U[, l], r)
# UL <- derivs2norm(U[, k], L[, l], r)
# LU <- derivs2norm(L[, k], U[, l], r)
# LL <- derivs2norm(L[, k], L[, l], r)
# (UU - UL - LU + LL) * rho$x[[k]][ind, ]
# }
# derivcornorm <- function(k, l, r, U, L) {
# -(corrderivnorm(r, x = U[, k], y = U[, l]) -
# corrderivnorm(r, x = U[, k], y = L[, l]) -
# corrderivnorm(r, x = L[, k], y = U[, l]) +
# corrderivnorm(r, x = L[, k], y = L[, l]))
# }
#
# dtheta[ind,] <- do.call("cbind", lapply(1:rho$ndim, function(j){
# if (j %in% comb){
# r <- sapply(sigmas[lev],'[', comb[1], comb[2])
# pr <- rectbiv.norm.prob(Uq[, c(j, comb[comb!=j])], Lq[, c(j, comb[comb!=j])], r)
# pr[pr < .Machine$double.eps] <- .Machine$double.eps
# derivtheta2norm(j, comb[comb!=j], U=Uq, L=Lq, r=r)/pr
# } else {
# matrix(0, nrow = length(ind), ncol = rho$npar.theta.opt[j])
# } }))
# dbeta[ind, ] <- do.call("cbind", lapply(1:rho$ndim, function(j){
# if (j %in% comb){
# r <- sapply(sigmas[lev],'[', comb[1], comb[2])
# pr <- rectbiv.norm.prob(Uq[, c(j, comb[comb!=j])], Lq[, c(j, comb[comb!=j])], r)
# pr[pr < .Machine$double.eps] <- .Machine$double.eps
# derivbeta2norm(j, comb[comb!=j], Uq, Lq, r)/pr
# } else {
# matrix(0, nrow = length(ind), ncol = rho$p)
# }
# }))
# poslev <- which(combn(rho$npar.cor, 2, simplify=F) %in% list(comb))
# for (l in 1:rho$ncor.levels){
# r <- sigmas[[l]][comb[1], comb[2]]
# pr <- rectbiv.norm.prob(Uq[lev == l, comb], Lq[lev == l, comb], r)
# pr[pr < .Machine$double.eps] <- .Machine$double.eps
# dc <- rep(0, length(ind))
# dc[lev == l] <- derivcornorm(comb[1], comb[2], r = r, Uq[lev == l,], Lq[lev == l, ])/pr
# dcorr[ind, (l-1) * rho$npar.cor + poslev] <- dc
# }
# }
# }
# H_comb[[i]] <- crossprod(cbind(dtheta, dbeta, dcorr))
# }
# i<-1
# Hi_comb <- list()
# aaa <- Reduce("+", H_comb)
# HH_comb <- Reduce("+", Hi_comb)
# head(HH_comb)
# str(H_comb)
# H_comb[[1]][1:3,1:3]
# H_comb[[2]][1:3,1:3]
# H_comb[[3]][1:3,1:3]
# Hnum[1:3,1:3]
# aaa[1:3,1:3]
# H_comb[[1]][4:6,4:6]
# H_comb[[2]][4:6,4:6]
# H_comb[[3]][4:6,4:6]
# Hnum[4:6,4:6]
# aaa[4:6,4:6]
# H_comb[[1]][7:11,7:11]
# H_comb[[2]][7:11,7:11]
# H_comb[[3]][7:11,7:11]
# round(Hnum[7:11,7:11],3)
# aaa[7:11,7:11]
#
# aaa[1:11, 1:11]
# H_comb[[1]][12:15, 12:15]
# H_comb[[2]][12:15, 12:15]
# H_comb[[3]][12:15, 12:15]
# round(Hnum[12:15, 12:15],3)
# aaa[12:15, 12:15]
#
#
# Hnum <- t(J.inv) %*% Ht %*% J.inv
#
# head(aaa)
# all.equal(aaa,HH_comb)
# sum(solve(aaa) == rho$H.inv)
# Hnum <- solve(rho$H.inv)
# head(round(Hnum,3))
# head(round(aaa,3))
# Hnum[1:11, 1:11]
# aaa[1:11, 1:11]
# Hnum[12:20, 12:20]
# aaa[12:20, 12:20]
# Hnum[21:23, 21:23]
# aaa[21:23, 21:23]
# ## no constraints, nada !!!
#
# r <- sapply(sigmas[lev],'[', comb[1], comb[2])
# pr <- rectbiv.norm.prob(Uq[, c(j, comb[comb!=j])], Lq[, c(j, comb[comb!=j])], r)
# dtheta[ind, ] <- derivtheta2norm(j, comb[comb!=j], U=Uq, L=Lq, r=r)/pr
#
#
#
# dbeta[ind, ] <- do.call("cbind", lapply(1:rho$ndim, function(j){
# combisj <- combis[sapply(combis, function(x) j %in% x)]
# if (length(combisj) != 0) {
# Reduce("+", lapply(combisj, function(comb){
# r <- sapply(sigmas[lev],'[', comb[1], comb[2])
# pr <- rectbiv.norm.prob(Uq[, c(j, comb[comb!=j])], Lq[, c(j, comb[comb!=j])], r)
# pr[pr < .Machine$double.eps] <- .Machine$double.eps
# derivbeta2norm(j, comb[comb!=j], Uq, Lq, r)/pr}))
# } else {
# matrix(0, nrow = length(ind), ncol = rho$p)
# }
# }))
# poslev <- which(combn(rho$npar.cor, 2, simplify=F) %in% combis)
# for (l in 1:rho$ncor.levels){
# dcorr[ind, (l-1) * rho$npar.cor + poslev] <- sapply(combis, function(comb) {
# r <- sigmas[[l]][comb[1], comb[2]]
# pr <- rectbiv.norm.prob(Uq[lev == l, comb], Lq[lev == l, comb], r)
# pr[pr < .Machine$double.eps] <- .Machine$double.eps
# dc <- rep(0, length(ind))
# dc[lev == l] <- derivcornorm(comb[1], comb[2], r = r, Uq[lev == l,], Lq[lev == l, ])/pr
# dc
# })
# }
# }
# }
# Vii <- cbind(dtheta, dbeta, dcorr)
# return(Vii)
# }
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.