R/Qcoxph_ee.R

Defines functions A_beta_gamma_fun A_gamma_gamma_fun A_gamma_gamma_mix_fun A_beta_beta_fun A_gamma_fun A_beta_fun coxBB QCox.ee.ct QCox.ee.c

Documented in A_beta_beta_fun A_beta_fun A_beta_gamma_fun A_gamma_fun A_gamma_gamma_fun A_gamma_gamma_mix_fun coxBB QCox.ee.c QCox.ee.ct

# Qcoxph estimating equations (building blocks below)
# The argument "eps" is not used, but is required for compatibility with Qest.newton.

# Estimating equation + Jacobian of a Cox model with censored data.
QCox.ee.c <- function(theta, eps, z, y, d, X, w, knots, tau, J = FALSE, EE){

  n <- nrow(X)
  q <- ncol(X)
  k <- length(knots) - 2 # n. of internal knots
  Ltau <- function(opt, tau){opt$tau <- tau; opt}

  #############################################################################
  #############################################################################
  # Estimating equation #######################################################
  #############################################################################
  #############################################################################

  if(missing(EE)){

    # Building blocks
    BBy <- coxBB(theta, y, X, knots, tau)

    # A_beta(p) - to be A*X
    Ab <- A_beta_fun(BBy)
    A_beta <- Ab$A
    AA_beta <- Ab$AA
    AA1_beta <- Ab$AA1

    # A_gamma(p)
    Ag <- A_gamma_fun(BBy)
    A_gamma <- Ag$A
    AA_gamma <- Ag$AA
    AA1_gamma <- Ag$AA1

    #############################################################################
    # Estimating equation #######################################################
    #############################################################################

    Ay <- cbind(A_beta*X, A_gamma) # A_theta
    AAy <- cbind(AA_beta*X, AA_gamma) # AA_theta
    AA1 <- cbind(AA1_beta*X, AA1_gamma) # AA1_theta
    g.i <- -d*Ay + (1 - d)/BBy$Sy*(AAy - AA1) + AA1

    # Regularization. If F(y) = 1, the estimating equation has a 0/0 = NaN.

    outLy <- BBy$outL; outRy <- BBy$outR
    g.i[outLy,] <- (d*AA1)[outLy,]
    g.i[outRy,] <- (AA1 - d*Ay)[outRy,]
    g <- c(w %*% g.i)
  }
  else{
    g <- EE$g; g.i <- EE$g.i
    BBy <- EE$BBy; Ay <- EE$Ay; AAy <- EE$AAy; AA1 <- EE$AA1
  }

  #############################################################################
  #############################################################################
  # Jacobian ##################################################################
  #############################################################################
  #############################################################################

  if(J){

    # A_beta_beta(p) - to be t(A*X)%*%X
    Abb <- A_beta_beta_fun(BBy)

    # A_gamma_jh(p)
    if(k > 0){Agmix <- A_gamma_gamma_mix_fun(BBy)}

    # A_gamma_jj(p)
    Agg <- A_gamma_gamma_fun(BBy)

    # A_beta_gamma(p) - to be t(x)%*%A to get a matrix q*(k + 1)
    Abg <- A_beta_gamma_fun(BBy)


    #############################################################################
    #############################################################################
    # Assemble the Jacobian #####################################################
    #############################################################################
    #############################################################################

    genH <- function(A, d, U, outL, outR){
      H <- -d*A$A + U*(A$AA - A$AA1) + A$AA1
      H <- cbind(H)
      H[outLy,] <- d[outLy]*cbind(A$AA1)[outLy,]
      H[outRy,] <- cbind(A$AA1)[outRy,] - d[outRy]*cbind(A$A)[outRy,]
      H
    }

    Xw <- X*w
    Sy <- BBy$Sy; Fy <- BBy$Fy; fy <- BBy$fy
    U <- (1 - d)/Sy
    outLy <- BBy$outL; outRy <- BBy$outR

    H_beta <- t(c(genH(Abb, d, U, outLy, outRy))*X)%*%Xw
    H_gamma_mix <- (if(k > 0) .colSums(genH(Agmix, d, U, outLy, outRy)*w, n, k*(k + 1)/2) else NULL)
    H_gamma_jj <- .colSums(genH(Agg, d, U, outLy, outRy)*w, n, k + 1)
    H_beta_gamma <- t(Xw)%*%genH(Abg, d, U, outLy, outRy)


    # Assemble

    H0 <- matrix(NA, q + k + 1, q + k + 1)
    i_beta <- (1:q); i_gamma <- (q + 1):(q + k + 1)
    H0[i_beta, i_beta] <- H_beta
    H0[i_beta, i_gamma] <- H_beta_gamma; H0[i_gamma, i_beta] <- t(H_beta_gamma)
    H0gamma <- matrix(0, k + 1, k + 1)
    if(k > 0){H0gamma[lower.tri(H0gamma)] <- H_gamma_mix; H0gamma <- H0gamma + t(H0gamma)}
    diag(H0gamma) <- H_gamma_jj
    H0[i_gamma, i_gamma] <- H0gamma

    fy[c(outLy, outRy)] <- 0
    wy <- do.call(tau$wfun, Ltau(tau$opt, Fy))
    H1 <- -fy*BBy$Qtheta
    H2 <- -d*wy*BBy$Qtheta + U*(Ay + (AAy - AA1)/Sy)

    # Finish

    H <- crossprod(w*H1, H2) + H0
    J <- H # So if J = FALSE, return J = FALSE; and if J = TRUE, return the jacobian
  }

  list(g = g, g.i = g.i, J = J,
       Fy = BBy$Fy, fy = BBy$fy,
       Fz = NULL, fz = NULL,
       BBy = BBy, Ay = Ay, AAy = AAy,
       BBz = NULL, Az = NULL, AAz = NULL,
       AA1 = AA1
  )
}

# Estimating equation + Jacobian of a Cox model with censored and truncated data.
QCox.ee.ct <- function(theta, eps, z, y, d, X, w, knots, tau, J = FALSE, EE){

  n <- nrow(X)
  q <- ncol(X)
  k <- length(knots) - 2 # n. of internal knots
  Ltau <- function(opt, tau){opt$tau <- tau; opt}

  #############################################################################
  #############################################################################
  # Estimating equation #######################################################
  #############################################################################
  #############################################################################

  if(missing(EE)){

    # Building blocks
    BBy <- coxBB(theta, y, X, knots, tau)
    BBz <- coxBB(theta, z, X, knots, tau)

    # A_beta(p) - to be A*X
    Aby <- A_beta_fun(BBy); Abz <- A_beta_fun(BBz)
    A_betay <- Aby$A; A_betaz <- Abz$A
    AA_betay <- Aby$AA; AA_betaz <- Abz$AA
    AA1_beta <- Aby$AA1

    # A_gamma(p)
    Agy <- A_gamma_fun(BBy); Agz <- A_gamma_fun(BBz)
    A_gammay <- Agy$A; A_gammaz <- Agz$A
    AA_gammay <- Agy$AA; AA_gammaz <- Agz$AA
    AA1_gamma <- Agy$AA1

    #############################################################################
    # Estimating equation #######################################################
    #############################################################################

    Sy <- BBy$Sy; Sz <- BBz$Sy; Fy <- BBy$Fy; Fz <- BBz$Fy

    Ay <- cbind(A_betay*X, A_gammay) # A_theta(Fy)
    AAy <- cbind(AA_betay*X, AA_gammay) # AA_theta(Fy)
    Az <- cbind(A_betaz*X, A_gammaz) # A_theta(Fz)
    AAz <- cbind(AA_betaz*X, AA_gammaz) # AA_theta(Fz)
    AA1 <- cbind(AA1_beta*X, AA1_gamma) # AA1_theta
    gy.i <- -d*Ay + (1 - d)/Sy*(AAy - AA1)
    gz.i <- -1/Sz*(AAz - AA1)

    # Regularization. If F(y) = 1 or F(z) = 1, the estimating equation has a 0/0 = NaN.

    outLy <- BBy$outL; outRy <- BBy$outR
    outLz <- BBz$outL; outRz <- BBz$outR
    gy.i[outLy,] <- (-(1 - d)*AA1)[outLy,]
    gz.i[outLz,] <- AA1[outLz,]
    gy.i[outRy,] <- (-d*Ay)[outRy,]
    gz.i[outRz,] <- 0

    g.i <- gy.i + gz.i
    g <- c(w %*% g.i)
  }
  else{
    g <- EE$g; g.i <- EE$g.i
    BBy <- EE$BBy; BBz <- EE$BBz
    Ay <- EE$Ay; AAy <- EE$AAy; Az <- EE$Az; AAz <- EE$AAz; AA1 <- EE$AA1
  }

  #############################################################################
  #############################################################################
  # Jacobian ##################################################################
  #############################################################################
  #############################################################################

  if(J){

    # A_beta_beta(p) - to be t(A*X)%*%X
    Abby <- A_beta_beta_fun(BBy)
    Abbz <- A_beta_beta_fun(BBz)

    # A_gamma_jh(p)
    if(k > 0){
      Agmixy <- A_gamma_gamma_mix_fun(BBy)
      Agmixz <- A_gamma_gamma_mix_fun(BBz)
    }

    # A_gamma_jj(p)
    Aggy <- A_gamma_gamma_fun(BBy)
    Aggz <- A_gamma_gamma_fun(BBz)

    # A_beta_gamma(p) - to be t(x)%*%A to get a matrix q*(k + 1)
    Abgy <- A_beta_gamma_fun(BBy)
    Abgz <- A_beta_gamma_fun(BBz)

    #############################################################################
    #############################################################################
    # Assemble the Jacobian #####################################################
    #############################################################################
    #############################################################################

    genH <- function(Ay, Az, d, Uy, Uz, outLy, outLz, outRy, outRz){

      Hy <- -d*Ay$A + Uy*(Ay$AA - Ay$AA1)
      Hz <- -Uz*(Az$AA - Az$AA1)
      Hy <- cbind(Hy); Hz <- cbind(Hz)

      Hy[outLy,] <- -(1 - d[outLy])*cbind(Ay$AA1)[outLy,]
      Hz[outLz,] <- cbind(Az$AA1)[outLz,]
      Hy[outRy,] <- -d[outRy]*cbind(Ay$A)[outRy,]
      Hz[outRz,] <- 0

      Hy + Hz
    }

    Xw <- X*w
    Sy <- BBy$Sy; Sz <- BBz$Sy; Fy <- BBy$Fy; Fz <- BBz$Fy; fy <- BBy$fy; fz <- BBz$fy
    Uy <- (1 - d)/Sy; Uz <- 1/Sz
    outLy <- BBy$outL; outRy <- BBy$outR
    outLz <- BBz$outL; outRz <- BBz$outR

    H_beta <- t(c(genH(Abby, Abbz, d, Uy, Uz, outLy, outLz, outRy, outRz))*X)%*%Xw
    H_gamma_mix <- (if(k > 0) .colSums(genH(Agmixy, Agmixz, d, Uy, Uz, outLy, outLz, outRy, outRz)*w, n, k*(k + 1)/2) else NULL)
    H_gamma_jj <- .colSums(genH(Aggy, Aggz, d, Uy, Uz, outLy, outLz, outRy, outRz)*w, n, k + 1)
    H_beta_gamma <- t(Xw)%*%genH(Abgy, Abgz, d, Uy, Uz, outLy, outLz, outRy, outRz)


    # Assemble

    H0 <- matrix(NA, q + k + 1, q + k + 1)
    i_beta <- (1:q); i_gamma <- (q + 1):(q + k + 1)
    H0[i_beta, i_beta] <- H_beta
    H0[i_beta, i_gamma] <- H_beta_gamma; H0[i_gamma, i_beta] <- t(H_beta_gamma)
    H0gamma <- matrix(0, k + 1, k + 1)
    if(k > 0){H0gamma[lower.tri(H0gamma)] <- H_gamma_mix; H0gamma <- H0gamma + t(H0gamma)}
    diag(H0gamma) <- H_gamma_jj
    H0[i_gamma, i_gamma] <- H0gamma

    fy[c(outLy, outRy)] <- 0
    fz[c(outLz, outRz)] <- 0
    wy <- do.call(tau$wfun, Ltau(tau$opt, Fy))
    H1y <- -fy*BBy$Qtheta
    H2y <- -d*wy*BBy$Qtheta + Uy*(Ay + (AAy - AA1)/Sy)
    H1z <- -fz*BBz$Qtheta
    H2z <- Uz*(Az + (AAz - AA1)/Sz)

    # Finish

    H <- crossprod(w*H1y, H2y) - crossprod(w*H1z, H2z) + H0
    J <- H # So if J = FALSE, return J = FALSE; and if J = TRUE, return the jacobian
  }

  list(g = g, g.i = g.i, J = J,
       Fy = BBy$Fy, fy = BBy$fy,
       Fz = BBz$Fy, fz = BBz$fy,
       BBy = BBy, Ay = Ay, AAy = AAy,
       BBz = BBz, Az = Az, AAz = AAz,
       AA1 = AA1
  )
}



# Building blocks for Cox model

coxBB <- function(theta, y, X, knots, tau){

  n <- nrow(X)
  q <- ncol(X)
  k <- length(knots) - 2 # n. of internal knots
  beta <- (theta[1:q]) # Cox regression coefficients
  egamma <- exp(theta[(q + 1):length(theta)]) # coefficients of H0

  delta <- function(obj, e, k){t(t(obj[,2:(k + 2)] - obj[,1:(k + 1)])/e)}
  deltaj <- function(obj, e, b, k){o <- list(); for(j in 1:(k + 1)){o[[j]] <- delta(t(t(obj)*(e[j]*b[,j])), e, k)}; o}

  #############################################################################
  #############################################################################
  #############################################################################

  eta <- c(X%*%beta); HR <- exp(eta); HR1 <- 1/HR

  # Building blocks - part 1
  by <- plfcox(y, knots, deriv = 1) # b(y)
  H0y <- c(by$b%*%egamma) # H0(y)
  Fystar <- eta + log(H0y) # log(-log(1 - F(y | x)))
  outL <- which(Fystar <= tau$taustarL); outR <- which(Fystar >= tau$taustarR)
  Fystar <- pmax(pmin(Fystar, tau$taustarR), tau$taustarL)
  Hy <- exp(Fystar) # H(y | x)
  Sy <- exp(-Hy) # S(y | x)
  Fy <- 1 - Sy # F(y | x)


  # Building blocks - part 2
  bknots <- plfcox(knots, knots, deriv = 0) # b(knots)
  H0knots <- c(bknots%*%egamma) # H0(knots)
  Fknotstar <- eta + rep(1,n)%*%t(log(H0knots)) # log(-log(1 - F(knots | x)))
  Fknotstar[,k + 2] <- tau$taustarR # "as if" the last knot was Inf ("use all quantiles"). This is fundamental.
  Fknotstar <- pmax(pmin(Fknotstar, tau$taustarR), tau$taustarL)
  Hknots <- exp(Fknotstar) # H(knots | x)
  Sknots <- exp(-Hknots) # S(knots | x)
  Fknots <- 1 - Sknots # F(knots | x)


  # Building blocks - part 3 - Derivatives of the quantile function
  h0 <- c(by$b1 %*% egamma) # h0(y)
  hy <- h0*HR # h(y | x)
  fy <- hy*Sy # f(y | x)
  Qbeta <- -X*(Hy/hy)
  Qgamma <- -t(t(by$b/h0)*egamma)
  Qtheta <- cbind(Qbeta, Qgamma)


  # Building blocks - part 4

  ii <- (Fknotstar < Fystar)
  mFknotsp <- pmax0(Sknots - Sy)
  mFknots1 <- Sknots

  mins <- pmin(Fknots, Fy)
  minstar <- pmin(Fknotstar, Fystar)

  mins <- c(mins)
  minstar <- c(minstar)
  Fknots <- c(Fknots)
  Fknotstar <- c(Fknotstar)

  # w blocks

  w1mins <- matrix(tau$w1fun(mins), n, k + 2)
  wmins <- matrix(do.call(tau$wfun, Ltau(tau$opt, mins)), n, k + 2)
  Wmins <- matrix(tau$Wfun(mins), n, k + 2)
  iWmins <- matrix(tau$iWfun(mins), n, k + 2)

  w1knots <- matrix(tau$w1fun(Fknots), n, k + 2)
  w1knotsp <- w1knots*mFknotsp
  w1knots1 <- w1knots*mFknots1

  wknots <- matrix(do.call(tau$wfun, Ltau(tau$opt, Fknots)), n, k + 2)
  wknotsp <- wknots*mFknotsp
  wknots1 <- wknots*mFknots1

  Wknots <- matrix(tau$Wfun(Fknots), n, k + 2)
  Wknotsp <- Wknots*mFknotsp
  Wknots1 <- Wknots*mFknots1

  iWknots <- matrix(tau$iWfun(Fknots), n, k + 2)
  iWknotsp <- iWknots*mFknotsp

  wmins_2wknots <- wmins - 2*wknots
  Wmins_Wknots <- Wmins - Wknots
  w1knots1_wknots <- w1knots1 - wknots

  dWmins <- delta(Wmins, egamma, k)
  diWmins_dWknotsp <- delta(iWmins + Wknotsp, egamma, k)
  diWknots_dWknots1 <- delta(iWknots + Wknots1, egamma, k)


  # r blocks

  r1mins <- matrix(do.call(tau$r1fun, Ltau(tau$opt, mins)), n, k + 2)*HR1
  rmins <- matrix(do.call(tau$rfun, Ltau(tau$opt, minstar)), n, k + 2)*HR1
  Rmins <- matrix(tau$Rfun(minstar), n, k + 2)*HR1
  iRmins <- matrix(tau$iRfun(minstar), n, k + 2)*HR1

  r1knots <- matrix(do.call(tau$r1fun, Ltau(tau$opt, Fknots)), n, k + 2)*HR1
  r1knotsp <- r1knots*mFknotsp
  r1knots1 <- r1knots*mFknots1

  rknots <- matrix(do.call(tau$rfun, Ltau(tau$opt, Fknotstar)), n, k + 2)*HR1
  rknotsp <- rknots*mFknotsp
  rknots1 <- rknots*mFknots1

  Rknots <- matrix(tau$Rfun(Fknotstar), n, k + 2)*HR1
  Rknotsp <- Rknots*mFknotsp
  Rknots1 <- Rknots*mFknots1

  iRknots <- matrix(tau$iRfun(Fknotstar), n, k + 2)*HR1
  iRknotsp <- iRknots*mFknotsp

  rmins_2rknots <- rmins - 2*rknots
  Rmins_Rknots <- Rmins - Rknots
  r1knots1_rknots <- r1knots1 - rknots

  dRmins <- delta(Rmins, egamma, k)
  diRmins_dRknotsp <- delta(iRmins + Rknotsp, egamma, k)
  diRknots_dRknots1 <- delta(iRknots + Rknots1, egamma, k)


  # Building blocks with Vbeta

  vbeta <- t(t(Sknots)*c(H0knots))*HR
  Vbeta <- vbeta*ii

  dwminsVbeta <- delta(wmins*Vbeta, egamma, k)
  dwknots1vbeta <- delta(wknots1*vbeta, egamma, k)
  Owbeta <- delta(Wmins_Wknots*Vbeta + wknotsp*vbeta, egamma, k)

  drminsVbeta <- delta(rmins*Vbeta, egamma, k)
  drknots1vbeta <- delta(rknots1*vbeta, egamma, k)
  Orbeta <- delta(Rmins_Rknots*Vbeta + rknotsp*vbeta, egamma, k)


  # Building blocks with Vgamma
  # Note: in practice, I need Vgamma_hj = Vgamma*exp(gamma[j])*bj(knots)

  vgamma <- Sknots*HR
  Vgamma <- vgamma*ii

  dwminsVgamma <- deltaj(wmins*Vgamma, egamma, bknots, k)
  dwknots1vgamma <- deltaj(wknots1*vgamma, egamma, bknots, k)
  Owgamma <- deltaj(Wmins_Wknots*Vgamma + wknotsp*vgamma, egamma, bknots, k)

  drminsVgamma <- deltaj(rmins*Vgamma, egamma, bknots, k)
  drknots1vgamma <- deltaj(rknots1*vgamma, egamma, bknots, k)
  Orgamma <- deltaj(Rmins_Rknots*Vgamma + rknotsp*vgamma, egamma, bknots, k)


  #############################################################################
  # Output - only the strict necessary ########################################
  #############################################################################


  list(
   n = n, q = q, k = k,
   beta = beta, egamma = egamma, knots = knots,
   HR = HR, HR1 = HR1, fy = fy, Qtheta = Qtheta,
   by = by, Hy = Hy, Sy = Sy, Fy = Fy,
   bknots = bknots, H0knots = H0knots, Sknots = Sknots,
   outL = outL, outR = outR,
   ii = ii, c1 = knots[1:(k + 1)]*egamma - H0knots[1:(k + 1)],

   ###

   w1mins = w1mins, wmins = wmins, Wmins = Wmins,
   w1knotsp = w1knotsp, w1knots1 = w1knots1,
   wknots = wknots, wknotsp = wknotsp, wknots1 = wknots1,
   Wknots = Wknots,

   wmins_2wknots = wmins_2wknots,
   Wmins_Wknots = Wmins_Wknots,
   w1knots1_wknots = w1knots1_wknots,

   dWmins = dWmins,
   diWmins_dWknotsp = diWmins_dWknotsp,
   diWknots_dWknots1 = diWknots_dWknots1,

   ###

   r1mins = r1mins, rmins = rmins, Rmins = Rmins,
   r1knotsp = r1knotsp, r1knots1 = r1knots1,
   rknots = rknots, rknotsp = rknotsp, rknots1 = rknots1,
   Rknots = Rknots,

   rmins_2rknots = rmins_2rknots,
   Rmins_Rknots = Rmins_Rknots,
   r1knots1_rknots = r1knots1_rknots,

   dRmins = dRmins,
   diRmins_dRknotsp = diRmins_dRknotsp,
   diRknots_dRknots1 = diRknots_dRknots1,

   ###

   vbeta = vbeta,
   Vbeta = Vbeta,

   dwminsVbeta = dwminsVbeta,
   dwknots1vbeta = dwknots1vbeta,
   Owbeta = Owbeta,

   drminsVbeta = drminsVbeta,
   drknots1vbeta = drknots1vbeta,
   Orbeta = Orbeta,

   ###

   vgamma = vgamma,
   Vgamma = Vgamma,

   dwminsVgamma = dwminsVgamma,
   dwknots1vgamma = dwknots1vgamma,
   Owgamma = Owgamma,

   drminsVgamma = drminsVgamma,
   drknots1vgamma = drknots1vgamma,
   Orgamma = Orgamma
  )
}






# A_beta
A_beta_fun <- function(BB){

  n <- BB$n; k <- BB$k

  c1 <- BB$c1
  c2A <- BB$dRmins
  c2AA <- BB$diRmins_dRknotsp
  c2AA1 <- BB$diRknots_dRknots1

  #########################################

  A <- .rowSums(t(t(BB$dwminsVbeta)*c1) - c2A + BB$drminsVbeta, n, k + 1)
  AA <- .rowSums(t(t(BB$Owbeta)*c1) - c2AA + BB$Orbeta, n, k + 1)
  AA1 <- .rowSums(t(t(BB$dwknots1vbeta)*c1) - c2AA1 + BB$drknots1vbeta, n, k + 1)

  list(A = A, AA = AA, AA1 = AA1)
}




# A_gamma
A_gamma_fun <- function(BB){

  egamma <- BB$egamma
  n <- BB$n; k <- BB$k
  H0knots <- BB$H0knots[1:(k + 1)]

  # A ################################################################

  cAw <- BB$dWmins
  cAr <- BB$dRmins
  A <- t(t(cAw)*H0knots) - cAr

  # AA ################################################################

  cAAw <- BB$diWmins_dWknotsp
  cAAr <- BB$diRmins_dRknotsp
  AA <- t(t(cAAw)*H0knots) - cAAr

  # AA1 ################################################################

  cAA1w <- BB$diWknots_dWknots1
  cAA1r <- BB$diRknots_dRknots1
  AA1 <- t(t(cAA1w)*H0knots) - cAA1r

  if(k == 0){return(list(A = A, AA = AA, AA1 = AA1))}
  for(j in 1:k){ # When j = k + 1, the result is zero.

    ej <- egamma[j]*BB$bknots[,j]

    indh <- j:(k + 1) # Which "h" I will use for this "l"
    ejh <- ej[indh]
    egammah <- egamma[indh]
    c1h <- BB$c1[indh]

    # A ################################################################

    t1 <- t(t(BB$dwminsVgamma[[j]][,indh])*c1h)
    t2 <- (BB$drminsVgamma[[j]][,indh])
    t3 <- -t(t(BB$dWmins[,indh])*ejh)
    A[,j] <- A[,j] + .rowSums(t1 + t2 + t3, n, k - j + 2)

    # AA ################################################################

    t1 <- t(t(BB$Owgamma[[j]][,indh])*c1h)
    t2 <- (BB$Orgamma[[j]][,indh])
    t3 <- -t(t(BB$diWmins_dWknotsp[,indh])*ejh)
    AA[,j] <- AA[,j] + .rowSums(t1 + t2 + t3, n, k - j + 2)

    # AA1 ################################################################

    t1 <- t(t(BB$dwknots1vgamma[[j]][,indh])*c1h)
    t2 <- (BB$drknots1vgamma[[j]][,indh])
    t3 <- -t(t(BB$diWknots_dWknots1[,indh])*ejh)
    AA1[,j] <- AA1[,j] + .rowSums(t1 + t2 + t3, n, k - j + 2)
  }
  list(A = A, AA = AA, AA1 = AA1)
}






# A_beta_beta
A_beta_beta_fun <- function(BB){

  delta <- function(obj, egamma){k <- ncol(obj); t(t(obj[,2:k] - obj[,1:(k - 1)])/egamma)}
  n <- BB$n; k <- BB$k
  egamma <- BB$egamma
  c1 <- BB$c1

  # Preliminary quantities

  HH <- 1 - BB$HR%*%t(BB$H0knots)
  Vbeta <- BB$Vbeta; vbeta <- BB$vbeta
  Vbeta2 <- Vbeta^2
  vbeta2 <- vbeta^2
  Vbetastar <- Vbeta*HH
  vbetastar <- vbeta*HH

  # A ################################################################

  uw <- BB$w1mins*Vbeta2
  Uw <- BB$wmins*Vbetastar
  dw <- delta(uw + Uw, egamma)

  ur <- BB$r1mins*Vbeta2
  Ur <- BB$rmins*Vbetastar
  dr <- delta(ur + Ur, egamma)

  cr <- BB$dRmins - 2*BB$drminsVbeta

  A <- .rowSums(t(t(dw)*c1) + (dr + cr), n, k + 1)

  # AA ################################################################

  uw <- (BB$wmins_2wknots)*Vbeta2 + BB$w1knotsp*vbeta2
  Uw <- (BB$Wmins_Wknots)*Vbetastar + BB$wknotsp*vbetastar
  dw <- delta(uw + Uw, egamma)

  ur <- (BB$rmins_2rknots)*Vbeta2 + BB$r1knotsp*vbeta2
  Ur <- (BB$Rmins_Rknots)*Vbetastar + BB$rknotsp*vbetastar
  dr <- delta(ur + Ur, egamma)

  cr <- BB$diRmins_dRknotsp - 2*BB$Orbeta

  AA <- .rowSums(t(t(dw)*c1) + (dr + cr), n, k + 1)

  # AA1 ################################################################

  uw <- (BB$w1knots1_wknots)*vbeta2
  Uw <- BB$wknots1*vbetastar
  dw <- delta(uw + Uw, egamma)

  ur <- (BB$r1knots1_rknots)*vbeta2
  Ur <- BB$rknots1*vbetastar
  dr <- delta(ur + Ur, egamma)

  cr <- BB$diRknots_dRknots1 - 2*BB$drknots1vbeta

  AA1 <- .rowSums(t(t(dw)*c1) + (dr + cr), n, k + 1)

  #######################

  list(A = A, AA = AA, AA1 = AA1)
}





# A_(gamma_j, gamma_l) with l > j
A_gamma_gamma_mix_fun <- function(BB){

  delta <- function(obj, egamma){k <- ncol(obj); t(t(obj[,2:k] - obj[,1:(k - 1)])/egamma)}
  n <- BB$n; k <- BB$k
  egamma <- BB$egamma
  H0knots <- BB$H0knots[1:(k + 1)]

  A <- AA <- AA1 <- NULL
  for(j in 1:k){
    indl <- (j + 1):(k + 1) # Which "l" I will use for this "j"
    ej <- egamma[j]*BB$bknots[,j]

    # A ################################################################

    cjA.1 <- BB$dwminsVgamma[[j]]
    cjA.2 <- BB$dWmins
    cjA.3 <- BB$drminsVgamma[[j]]
    cjA <- t(t(cjA.1)*H0knots) + t(t(cjA.2)*ej[1:(k + 1)]) - cjA.3

    # AA ################################################################

    cjAA.1 <- BB$Owgamma[[j]]
    cjAA.2 <- BB$diWmins_dWknotsp
    cjAA.3 <- BB$Orgamma[[j]]
    cjAA <- t(t(cjAA.1)*H0knots) + t(t(cjAA.2)*ej[1:(k + 1)]) - cjAA.3

    # AA1 ################################################################

    cjAA1.1 <- BB$dwknots1vgamma[[j]]
    cjAA1.2 <- BB$diWknots_dWknots1
    cjAA1.3 <- BB$drknots1vgamma[[j]]
    cjAA1 <- t(t(cjAA1.1)*H0knots) + t(t(cjAA1.2)*ej[1:(k + 1)]) - cjAA1.3

    for(l in indl){

      A_jl <- cjA[,l]
      AA_jl <- cjAA[,l]
      AA1_jl <- cjAA1[,l]

      if(l < k + 1){ # When l = k + 1, the result is zero.

        indh <- l:(k + 1) # Which "h" I will use for this "l"
        indh2 <- l:(k + 2) # Again, including last
        el <- egamma[l]*BB$bknots[,l]

        v_hlj <- -BB$vgamma[,indh2]*(BB$HR%*%t(ej[indh2]*el[indh2]))
        v_hl_hj <- -v_hlj*BB$Sknots[,indh2]
        V_hlj <- v_hlj*BB$ii[,indh2]
        V_hl_hj <- v_hl_hj*BB$ii[,indh2]

        ejh <- ej[indh]; elh <- el[indh]
        egammah <- egamma[indh]
        c1h <- BB$c1[indh]

        # A ################################################################

        w_V_hl_hj <- BB$w1mins[,indh2]*V_hl_hj
        w_V_hlj <- BB$wmins[,indh2]*V_hlj

        t1 <- delta(w_V_hl_hj + w_V_hlj, egammah)
        t1 <- t(t(t1)*c1h)

        r_V_hl_hj <- BB$r1mins[,indh2]*V_hl_hj
        r_V_hlj <- BB$rmins[,indh2]*V_hlj
        t2 <- (r_V_hl_hj + r_V_hlj)
        t2 <- delta(t2, egammah)

        dw_mixj <- BB$dwminsVgamma[[j]][,indh]
        dw_mixl <- BB$dwminsVgamma[[l]][,indh]
        t3 <- -(t(t(dw_mixj)*elh) + t(t(dw_mixl)*ejh))

        A_jl <- A_jl + .rowSums(t1 + t2 + t3, n, k - l + 2)

        # AA ################################################################

        w_V_hl_hj <- (BB$wmins_2wknots)[,indh2]*V_hl_hj
        w_V_hlj <- (BB$Wmins_Wknots)[,indh2]*V_hlj
        w_v_hl_hj <- BB$w1knotsp[,indh2]*v_hl_hj
        w_v_hlj <- BB$wknotsp[,indh2]*v_hlj
        t1 <- delta(w_V_hl_hj + w_V_hlj + w_v_hl_hj + w_v_hlj, egammah)
        t1 <- t(t(t1)*c1h)

        r_V_hl_hj <- (BB$rmins_2rknots)[,indh2]*V_hl_hj
        r_V_hlj <- (BB$Rmins_Rknots)[,indh2]*V_hlj
        r_v_hl_hj <- BB$r1knotsp[,indh2]*v_hl_hj
        r_v_hlj <- BB$rknotsp[,indh2]*v_hlj
        t2 <- (r_V_hl_hj + r_V_hlj + r_v_hl_hj + r_v_hlj)
        t2 <- delta(t2, egammah)

        dw_mixj <- BB$Owgamma[[j]][,indh]
        dw_mixl <- BB$Owgamma[[l]][,indh]
        t3 <- -(t(t(dw_mixj)*elh) + t(t(dw_mixl)*ejh))

        AA_jl <- AA_jl + .rowSums(t1 + t2 + t3, n, k - l + 2)

        # AA1 ################################################################

        w_v_hl_hj <- (BB$w1knots1_wknots)[,indh2]*v_hl_hj
        w_v_hlj <- BB$wknots1[,indh2]*v_hlj
        t1 <- delta(w_v_hl_hj + w_v_hlj, egammah)
        t1 <- t(t(t1)*c1h)

        r_v_hl_hj <- (BB$r1knots1_rknots)[,indh2]*v_hl_hj
        r_v_hlj <- BB$rknots1[,indh2]*v_hlj
        t2 <- (r_v_hl_hj + r_v_hlj)
        t2 <- delta(t2, egammah)

        dw_mixj <- BB$dwknots1vgamma[[j]][,indh]
        dw_mixl <- BB$dwknots1vgamma[[l]][,indh]
        t3 <- -(t(t(dw_mixj)*elh) + t(t(dw_mixl)*ejh))

        AA1_jl <- AA1_jl + .rowSums(t1 + t2 + t3, n, k - l + 2)
      }

      A <- cbind(A, A_jl)
      AA <- cbind(AA, AA_jl)
      AA1 <- cbind(AA1, AA1_jl)
    }
  }

  list(A = A, AA = AA, AA1 = AA1)
}




# A_(gamma_j, gamma_j)
A_gamma_gamma_fun <- function(BB){

  delta <- function(obj, egamma){k <- ncol(obj); t(t(obj[,2:k] - obj[,1:(k - 1)])/egamma)}
  listj <- function(x){out <- NULL; for(j in 1:length(x)){out <- cbind(out, x[[j]][,j])}; out}
  n <- BB$n; k <- BB$k

  egamma <- BB$egamma
  H0knots <- BB$H0knots[1:(k + 1)]
  c1 <- BB$c1

  # Preliminary quantities

  e <- egamma*BB$bknots[k + 2,] # Final value reached by each b(y)

  vgammaHR <- BB$vgamma*BB$HR
  vgamma2 <- BB$vgamma^2
  VgammaHR <- BB$Vgamma*BB$HR
  Vgamma2 <- BB$Vgamma^2

  v <- t(t(BB$vgamma[,-1])*(e/egamma)) # V(j + 1, j)/exp(gamma)
  V <- t(t(BB$Vgamma[,-1])*(e/egamma)) # V(j + 1, j)/exp(gamma)

  # A ################################################################

  cAw <- 2*(V*BB$wmins[,-1]) - BB$dWmins
  cAr <- 2*(V*BB$rmins[,-1]) - BB$dRmins
  A <- t(t(cAw)*H0knots) - cAr

  # AA ################################################################

  cAAw <- 2*(listj(BB$Owgamma)) - BB$diWmins_dWknotsp
  cAAr <- 2*(listj(BB$Orgamma)) - BB$diRmins_dRknotsp
  AA <- t(t(cAAw)*H0knots) - cAAr

  # AA1 ################################################################

  cAA1w <- 2*(listj(BB$dwknots1vgamma)) - BB$diWknots_dWknots1
  cAA1r <- 2*(listj(BB$drknots1vgamma)) - BB$diRknots_dRknots1
  AA1 <- t(t(cAA1w)*H0knots) - cAA1r


  if(k == 0){return(list(A = A, AA = AA, AA1 = AA1))}
  for(j in 1:k){ # When j = k + 1, the result is zero.

    indh <- j:(k + 1) # Which "h" I will use for this "l"
    indh2 <- j:(k + 2) # Again, including last

    ej <- egamma[j]*BB$bknots[,j]
    v_hjj <- -t(t(vgammaHR[,indh2])*ej[indh2])
    v2_hj <- t(t(vgamma2[,indh2])*(ej[indh2]^2))
    V_hjj <- -t(t(VgammaHR[,indh2])*ej[indh2])
    V2_hj <- t(t(Vgamma2[,indh2])*(ej[indh2]^2))

    ejh <- ej[indh]
    egammah <- egamma[indh]
    c1h <- c1[indh]

    # A ################################################################

    w_V_hjj <- BB$wmins[,indh2]*V_hjj
    w_V2_hj <- BB$w1mins[,indh2]*V2_hj
    t1 <- delta(w_V_hjj + w_V2_hj, egammah)
    t1 <- t(t(t1)*c1h)

    r_V_hjj <- BB$rmins[,indh2]*V_hjj
    r_V2_hj <- BB$r1mins[,indh2]*V2_hj
    t2 <- (r_V_hjj + r_V2_hj)
    t2 <- delta(t2, egammah)

    dw_mixj <- 2*BB$dwminsVgamma[[j]][,indh] + BB$dWmins[,indh]
    t3 <- -t(t(dw_mixj)*ejh)

    A[,j] <- A[,j] + .rowSums(t1 + t2 + t3, n, k - j + 2)

    # AA ################################################################

    w_V_hjj <- (BB$Wmins_Wknots[,indh2])*V_hjj
    w_V2_hj <- (BB$wmins_2wknots[,indh2])*V2_hj
    w_v_hjj <- BB$wknotsp[,indh2]*v_hjj
    w_v2_hj <- BB$w1knotsp[,indh2]*v2_hj
    t1 <- delta(w_V_hjj + w_V2_hj + w_v_hjj + w_v2_hj, egammah)
    t1 <- t(t(t1)*c1h)


    r_V_hjj <- (BB$Rmins_Rknots[,indh2])*V_hjj
    r_V2_hj <- (BB$rmins_2rknots[,indh2])*V2_hj
    r_v_hjj <- BB$rknotsp[,indh2]*v_hjj
    r_v2_hj <- BB$r1knotsp[,indh2]*v2_hj
    t2 <- (r_V_hjj + r_V2_hj + r_v_hjj + r_v2_hj)
    t2 <- delta(t2, egammah)

    dw_mixj <- 2*BB$Owgamma[[j]][,indh] + BB$diWmins_dWknotsp[,indh]
    t3 <- -t(t(dw_mixj)*ejh)

    AA[,j] <- AA[,j] + .rowSums(t1 + t2 + t3, n, k - j + 2)

    # AA1 ################################################################

    w_v_hjj <- BB$wknots1[,indh2]*v_hjj
    w_v2_hj <- (BB$w1knots1_wknots[,indh2])*v2_hj
    t1 <- delta(w_v_hjj + w_v2_hj, egammah)
    t1 <- t(t(t1)*c1h)

    r_v_hjj <- BB$rknots1[,indh2]*v_hjj
    r_v2_hj <- (BB$r1knots1_rknots[,indh2])*v2_hj
    t2 <- (r_v_hjj + r_v2_hj)
    t2 <- delta(t2, egammah)

    dw_mixj <- 2*BB$dwknots1vgamma[[j]][,indh] + BB$diWknots_dWknots1[,indh]
    t3 <- -t(t(dw_mixj)*ejh)

    AA1[,j] <- AA1[,j] + .rowSums(t1 + t2 + t3, n, k - j + 2)
  }
  list(A = A, AA = AA, AA1 = AA1)
}







# A_beta_gamma
A_beta_gamma_fun <- function(BB){

  delta <- function(obj, egamma){k <- ncol(obj); t(t(obj[,2:k] - obj[,1:(k - 1)])/egamma)}
  n <- BB$n; q <- BB$q; k <- BB$k
  egamma <- BB$egamma
  H0knots <- BB$H0knots[1:(k + 1)]
  c1 <- BB$c1

  # A ################################################################

  cAw <- BB$dwminsVbeta
  cAr <- BB$dRmins - BB$drminsVbeta
  A <- t(t(cAw)*H0knots) + cAr

  # AA ################################################################

  cAAw <- BB$Owbeta
  cAAr <- BB$diRmins_dRknotsp - BB$Orbeta
  AA <- t(t(cAAw)*H0knots) + cAAr

  # AA1 ################################################################

  cAA1w <- BB$dwknots1vbeta
  cAA1r <- BB$diRknots_dRknots1 - BB$drknots1vbeta
  AA1 <- t(t(cAA1w)*H0knots) + cAA1r

  if(k == 0){return(list(A = A, AA = AA, AA1 = AA1))}
  for(j in 1:k){ # When j = k + 1, the result is zero.

    ej <- egamma[j]*BB$bknots[,j]

    indh <- j:(k + 1) # Which "h" I will use for this "l"
    indh2 <- j:(k + 2) # Again, including last

    v_mixj <- t(t(-BB$vbeta[,indh2]*BB$HR + BB$vgamma[,indh2])*(ej[indh2]))
    vv_j <- t(t(BB$vbeta[,indh2]*BB$vgamma[,indh2])*(ej[indh2]))
    V_mixj <- v_mixj*BB$ii[,indh2]
    VV_j <- vv_j*BB$ii[,indh2]

    ejh <- ej[indh]
    egammah <- egamma[indh]
    c1h <- c1[indh]

    # A ################################################################

    w_VV_j <- BB$w1mins[,indh2]*VV_j
    w_V_mixj <- BB$wmins[,indh2]*V_mixj
    t1 <- delta(w_VV_j + w_V_mixj, egammah)
    t1 <- t(t(t1)*c1h)

    r_VV_j <- BB$r1mins[,indh2]*VV_j
    r_V_mixj <- BB$rmins[,indh2]*V_mixj
    t2 <- (r_VV_j + r_V_mixj)
    t2 <- delta(t2, egammah)

    dw_j <- BB$dwminsVbeta[,indh]
    dr_j <- BB$drminsVgamma[[j]][,indh]
    t3 <- -t(t(dw_j)*ejh) - dr_j

    A[,j] <- A[,j] + .rowSums(t1 + t2 + t3, n, k - j + 2)

    # AA ################################################################

    w_VV_j <- (BB$wmins[,indh2] - 2*BB$wknots[,indh2])*VV_j
    w_V_mixj <- (BB$Wmins[,indh2] - BB$Wknots[,indh2])*V_mixj
    w_vv_j <- BB$w1knotsp[,indh2]*vv_j
    w_v_mixj <- BB$wknotsp[,indh2]*v_mixj
    t1 <- delta(w_VV_j + w_V_mixj + w_vv_j + w_v_mixj, egammah)
    t1 <- t(t(t1)*c1h)

    r_VV_j <- (BB$rmins[,indh2] - 2*BB$rknots[,indh2])*VV_j
    r_V_mixj <- (BB$Rmins[,indh2] - BB$Rknots[,indh2])*V_mixj
    r_vv_j <- BB$r1knotsp[,indh2]*vv_j
    r_v_mixj <- BB$rknotsp[,indh2]*v_mixj
    t2 <- (r_VV_j + r_V_mixj + r_vv_j + r_v_mixj)
    t2 <- delta(t2, egammah)

    dw_j <- BB$Owbeta[,indh]
    dr_j <- BB$Orgamma[[j]][,indh]
    t3 <- -t(t(dw_j)*ejh) - dr_j

    AA[,j] <- AA[,j] + .rowSums(t1 + t2 + t3, n, k - j + 2)

    # AA1 ################################################################

    w_vv_j <- (BB$w1knots1[,indh2] - BB$wknots[,indh2])*vv_j
    w_v_mixj <- BB$wknots1[,indh2]*v_mixj
    t1 <- delta(w_vv_j + w_v_mixj, egammah)
    t1 <- t(t(t1)*c1h)

    r_vv_j <- (BB$r1knots1[,indh2] - BB$rknots[,indh2])*vv_j
    r_v_mixj <- BB$rknots1[,indh2]*v_mixj
    t2 <- (r_vv_j + r_v_mixj)
    t2 <- delta(t2, egammah)

    dw_j <- BB$dwknots1vbeta[,indh]
    dr_j <- BB$drknots1vgamma[[j]][,indh]
    t3 <- -t(t(dw_j)*ejh) - dr_j

    AA1[,j] <- AA1[,j] + .rowSums(t1 + t2 + t3, n, k - j + 2)

  }
  list(A = A, AA = AA, AA1 = AA1)
}

Try the Qest package in your browser

Any scripts or data that you put into this service are public.

Qest documentation built on May 29, 2024, 8:57 a.m.