R/AlternativeEstimator.R

In Francesco16p/brrr: Accurate estimation in relative risk regression

meth_alternative <- function(x,z=NULL,y,t,maxit=200, start=NULL, method= "Mle",slowit =1)
{
  # controll z
  if(is.null(z)) z <- x

  # transorm x,z in matrix
  x <- as.matrix(x)
  z <- as.matrix(z)
  #print(z)
  # re-ordering the data
  if(dim(x)[2]>1)
  {
    x <- rbind(x[t==0,], x[t==1,])
  }
  else
  {
    x <- matrix(c(x[t==0,], x[t==1,]), ncol = 1)
  }
  if(dim(z)[2]>1)
  {
    z <- rbind(z[t==0,], z[t==1,])
  }
  else
  {
    z <- matrix(c(z[t==0,], z[t==1,]), ncol = 1)
  }
  y <- c(y[t==0], y[t==1])
  t <- c(t[t==0], t[t==1])

  # getting indicies
  nobs <- nrow(x)
  nex <- nobs-sum(t) # number of exposed units
  nvarx <- ncol(x)
  nvarz <- ncol(z)
  if(is.null(start))
  {
    start<- rep(0,nvarx + nvarz)
  }

  # define x_delta and z_delta
  #x_delta <- array(NA, dim = c(n,n,p))
  #for(i in 1:nvarx) x_delta[,,i] <- diag(x[,i])
  #z_delta <- array(NA, dim = c(n,n,q))
  #for(i in 1:nvarz) z_delta[,,i] <- diag(z[,i])

  # defining functions for probabilities
  # function for Richardson's nusiance parameter
  prob <- function(eta1,eta2)
  {
    p0 <- rep(0,length(eta1))
    num <- (1+exp(eta2)*(1+exp(eta1))-sqrt((1+exp(eta2)*(1+exp(eta1)))^2-4*exp(eta1+2*eta2)))
    den<- 2*exp(eta1+eta2)
    p0 <- num/den
    #p0[p0< ep]<- ep
    #p0[p0>(1-ep)]<-(1-ep)
    return(p0)
  }


  # derivatives of prob.standard
  # derivative of pi0 with respect to eta1
  p0_dp0_eta1 <- function(eta1, eta2)
  {
    dp0.eta1 <-  ((exp(-eta2-eta1)*(exp(eta2+eta1)-(2*((exp(eta1)+1)*exp(eta2)+1)*exp(eta2+eta1)-4*exp(2*eta2+eta1))/
                                      (2*sqrt(((exp(eta1)+1)*exp(eta2)+1)^2-4*exp(2*eta2+eta1)))))/2-
                    (exp(-eta2-eta1)*(-sqrt(((exp(eta1)+1)*exp(eta2)+1)^2
                                            -4*exp(2*eta2+eta1))+(exp(eta1)+1)*exp(eta2)+1))/2)
    return(dp0.eta1)
  }
  # derivative of pi0 with respect to eta2
  # written in a different way compared to thesis
  p0_dp0_eta2 <- function(eta1,eta2)
  {
    dp0.eta2 <- ((exp(-eta2-eta1)*((exp(eta1)+1)*exp(eta2)-(2*(exp(eta1)+1)*exp(eta2)*((exp(eta1)+1)*exp(eta2)+1)-8*exp(2*eta2+eta1))/
                                     (2*sqrt(((exp(eta1)+1)*exp(eta2)+1)^2-4*exp(2*eta2+eta1)))))/2-
                   (exp(-eta2-eta1)*(-sqrt(((exp(eta1)+1)*exp(eta2)+1)^2-4*exp(2*eta2+eta1))+(exp(eta1)+1)*exp(eta2)+1))/2)
    return(dp0.eta2)
  }
  # Second derivative of pi0 with respect to eta1
  p0_d2p0_eta1 <- function(eta1,eta2)
  {
    d2p0.eta1<-((exp(-eta2-eta1)*(-(2*exp(2*eta2+2*eta1)-4*exp(2*eta2+eta1)+2*((exp(eta1)+1)*exp(eta2)+1)*exp(eta2+eta1))/
                                    (2*sqrt(((exp(eta1)+1)*exp(eta2)+1)^2-4*exp(2*eta2+eta1)))+(2*((exp(eta1)+1)*exp(eta2)+1)*exp(eta2+eta1)-4*exp(2*eta2+eta1))^2/
                                    (4*(((exp(eta1)+1)*exp(eta2)+1)^2-4*exp(2*eta2+eta1))^(3/2))+exp(eta2+eta1)))/2-
                  exp(-eta2-eta1)*(exp(eta2+eta1)-(2*((exp(eta1)+1)*exp(eta2)+1)*exp(eta2+eta1)-4*exp(2*eta2+eta1))/
                                     (2*sqrt(((exp(eta1)+1)*exp(eta2)+1)^2-4*exp(2*eta2+eta1))))+(exp(-eta2-eta1)*(-sqrt(((exp(eta1)+1)*exp(eta2)+1)^2-
                                                                                                                           4*exp(2*eta2+eta1))+(exp(eta1)+1)*exp(eta2)+1))/2)
    return(d2p0.eta1)
  }
  # Second derivative of pi0 with respect to eta2
  p0_d2p0_eta2 <- function(eta1,eta2)
  {
    d2p0.eta2<- ((exp(-eta2-eta1)*(-sqrt(((exp(eta1)+1)*exp(eta2)+1)^2-4*exp(2*eta2+eta1))+(exp(eta1)+1)*exp(eta2)+1))/
                   2-exp(-eta2-eta1)*((exp(eta1)+1)*exp(eta2)-(2*(exp(eta1)+1)*exp(eta2)*((exp(eta1)+1)*exp(eta2)+1)-8*exp(2*eta2+eta1))/
                                        (2*sqrt(((exp(eta1)+1)*exp(eta2)+1)^2-4*exp(2*eta2+eta1))))+(exp(-eta2-eta1)*
                                                                                                       (-(-16*exp(2*eta2+eta1)+2*(exp(eta1)+1)^2*exp(2*eta2)+2*(exp(eta1)+1)*exp(eta2)*((exp(eta1)+1)*exp(eta2)+1))/
                                                                                                          (2*sqrt(((exp(eta1)+1)*exp(eta2)+1)^2-4*exp(2*eta2+eta1)))+(2*(exp(eta1)+1)*exp(eta2)*((exp(eta1)+1)*exp(eta2)+1)-
                                                                                                                                                                        8*exp(2*eta2+eta1))^2/(4*(((exp(eta1)+1)*exp(eta2)+1)^2-4*exp(2*eta2+eta1))^(3/2))+(exp(eta1)+1)*exp(eta2)))/(2))
    return(d2p0.eta2)
  }
  # Second derivative of p0 respect to eta1 and eta2
  p0_d2p0_eta1_eta2 <- function(eta1,eta2)
  {
    d2p0.eta12<- ((exp(-eta2-eta1)*(-(2*(exp(eta1)+1)*exp(2*eta2+eta1)-8*exp(2*eta2+eta1)+2*((exp(eta1)+1)*exp(eta2)+1)*
                                        exp(eta2+eta1))/(2*sqrt(((exp(eta1)+1)*exp(eta2)+1)^2-4*exp(2*eta2+eta1)))+
                                      ((2*(exp(eta1)+1)*exp(eta2)*((exp(eta1)+1)*exp(eta2)+1)-8*exp(2*eta2+eta1))*(2*((exp(eta1)+1)*exp(eta2)+1)*
                                                                                                                     exp(eta2+eta1)-4*exp(2*eta2+eta1)))/(4*(((exp(eta1)+1)*exp(eta2)+1)^2-4*exp(2*eta2+eta1))^(3/2))+
                                      exp(eta2+eta1)))/(2)-(exp(-eta2-eta1)*(exp(eta2+eta1)-(2*((exp(eta1)+1)*exp(eta2)+1)*exp(eta2+eta1)-
                                                                                               4*exp(2*eta2+eta1))/(2*sqrt(((exp(eta1)+1)*exp(eta2)+1)^2-4*exp(2*eta2+eta1)))))/2+
                    (exp(-eta2-eta1)*(-sqrt(((exp(eta1)+1)*exp(eta2)+1)^2-4*exp(2*eta2+eta1))+(exp(eta1)+1)*exp(eta2)+1))/2-
                    (exp(-eta2-eta1)*((exp(eta1)+1)*exp(eta2)-(2*(exp(eta1)+1)*exp(eta2)*((exp(eta1)+1)*exp(eta2)+1)-8*exp(2*eta2+eta1))/
                                        (2*sqrt(((exp(eta1)+1)*exp(eta2)+1)^2-4*exp(2*eta2+eta1)))))/2)
    return(d2p0.eta12)
  }

  # Key quantities
  key_quantities <- function(pars)
  {
    # Basic quantities
    gammas <- pars[1:nvarx]
    betas <- pars[-c(1:nvarx)]
    eta1 <- drop(x %*% gammas)
    eta2 <- drop(z %*% betas)
    pi <- prob(eta1,eta2) * exp(eta1 * t)
    #pi[pi< ep] <- ep
    #pi[pi>1-ep] <- 1-ep
    d.1.eta1 <- (p0_dp0_eta1(eta1,eta2) * exp(eta1 * t) + pi * t)
    d.1.eta2 <- (p0_dp0_eta2(eta1,eta2) * exp(eta1 * t))
    d.2.eta1 <- (p0_d2p0_eta1(eta1,eta2) * exp(eta1 * t) + 2 *
                   p0_dp0_eta1(eta1,eta2) * exp(eta1 * t) * t + pi * t)
    d.2.eta2 <- (p0_d2p0_eta2(eta1,eta2) * exp(eta1 * t))
    d.2.eta12 <- (p0_d2p0_eta1_eta2(eta1,eta2) * exp(eta1 * t) +
                    p0_dp0_eta2(eta1,eta2) * exp(eta1 * t) * t)

    # Diagonal matrices

    # d.1.eta1 <- diag(d.1.eta1)
    #d.1.eta2 <- diag(d.1.eta2)
    #d.2.eta1 <- diag(d.2.eta1)
    #d.2.eta2 <- diag(d.2.eta2)
    #d.2.eta12 <- diag(d.2.eta12)

    # Quanties related to the variance

    v <- pi*(1-pi) # variance
    d1v <- 1-2*pi  # variance's derivative w.r.t pi
    #V <- diag(v)
    #d1v <- diag(d1v)

    out <- list(betas = betas, gammas = gammas, eta1 = eta1, eta2 = eta2,
                pi = pi, d.1.eta1 = d.1.eta1, d.1.eta2 = d.1.eta2,
                d.2.eta1 = d.2.eta1, d.2.eta2 = d.2.eta2,  d.2.eta12 = d.2.eta12,
                v = v, d1v = d1v)
    return(out)
  }
  # Score-function
  score <- function(pars,fit=NULL)
  {
    if(is.null(fit))
    {
      fit <- key_quantities(pars)
    }

    out <-with( fit,
                # score for each observation
                {
                  score_gamma <- ( x * d.1.eta1 * (y-pi) / v)
                  score_beta <-  ( z * d.1.eta2 * (y-pi) / v)
                  score_tot <- colSums(cbind(score_gamma, score_beta))
                  score_tot
                }
    )
    out
  }

  # Fisher expected information
  information <- function(pars,inverse = FALSE,fit = NULL)
  {
    if(is.null(fit))
    {
      fit <- key_quantities(pars)
    }
    info<- with(fit,{
      # defining the blocks of the expected information
      igg <- crossprod( x , (d.1.eta1 ^ 2 * x / v))
      igb <- crossprod( x , (d.1.eta1 * d.1.eta2 * z / v))
      ibb <- crossprod( z , (d.1.eta2 ^ 2 * z / v))
      info <- as.matrix(rbind( cbind(igg, igb),
                               cbind(t(igb), ibb)))
      info
    })
    if(!inverse)
    {
      return(info)
    }
    else
    {
      return(chol2inv(chol(info)))
    }
  }

  # nvarx + nvarz gamma matrix
  P_Q_gamma <- function(pars, fit = NULL)
  {
    if(is.null(fit))
    {
      fit <- key_quantities(pars)
    }
    PQg <- with(fit,{
      # parts costant in the matrix
      PQg <- array(data = NA, dim= c(nvarx + nvarz, nvarx + nvarz, nvarx))
      for(t in 1:nvarx)
      {
        xt =x[, t]
        L.gg.t <- crossprod( x ,( d.2.eta1 * d.1.eta1 / v ) * xt * x)
        L.gb.t <- crossprod( x ,( d.2.eta12 * d.1.eta1 / v ) * xt * z)
        L.bb.t <- crossprod( z ,( d.2.eta2 * d.1.eta1 / v  ) * xt * z)

        PQg[,,t] <- as.matrix( rbind( cbind(L.gg.t,L.gb.t),
                                      cbind(t(L.gb.t),L.bb.t ) ))


      }
      PQg
    })
    return(PQg)
  }

  # P+ Q beta matrix
  P_Q_beta <- function(pars, fit = NULL)
  {
    if(is.null(fit))
    {
      fit <- key_quantities(pars)
    }
    PQb <- with(fit,{
      # parts costant in the matrix
      PQb <- array(data = NA, dim= c(nvarx + nvarz,nvarx + nvarz, nvarz))

      for(s in 1:nvarz)
      {
        zs <- z[,s]
        L.gg.s <- crossprod( x ,( d.2.eta1 * d.1.eta2 / v )*zs * x)
        L.gb.s <- crossprod( x ,( d.2.eta12 * d.1.eta2 / v ) * zs * z)
        L.bb.s <- crossprod( z ,( d.2.eta2 * d.1.eta2 / v) * zs * z)
        PQb[,,s] <- as.matrix( rbind( cbind(L.gg.s,L.gb.s),
                                      cbind(t(L.gb.s),L.bb.s)))
      }
      PQb
    })
    return(PQb)
  }
  # P/3+Q/2 gamma matrix
  P_Q_3_2_gamma <- function(pars, fit = NULL)
  {
    if(is.null(fit))
    {
      fit <- key_quantities(pars)
    }
    PQ32g <- with(fit,{
      # parts constant in the matrix
      PQ32g <- array(data = NA, dim= c(nvarx + nvarz,nvarx + nvarz, nvarx))

      for(t in 1:nvarx)
      {
        xt <- x[,t]
        L.gg.t <- crossprod( x ,( 1 / 3 * d1v * d.1.eta1 ^ 3 / v ^ 2 +
                                    1 / 2 * (d.2.eta1 * d.1.eta1 / v - d1v *
                                               d.1.eta1 ^ 3 / v ^ 2)) * xt * x)
        L.gb.t <- crossprod( x ,( 1 / 3 * d1v * d.1.eta1 ^ 2 * d.1.eta2 / v ^ 2 +
                                    1 / 2 * ( d.2.eta12 * d.1.eta1 / v - d1v *
                                                d.1.eta1 ^ 2 * d.1.eta2 / v ^ 2)) * xt * z)
        L.bb.t <- crossprod( z ,( 1 / 3 * d1v * d.1.eta1 * d.1.eta2 ^ 2 / v ^ 2 +
                                    1 / 2 * (d.2.eta2 * d.1.eta1 / v - d1v *
                                               d.1.eta1 * d.1.eta2 ^ 2 / v ^ 2)) * xt * z)

        PQ32g[,,t] <- as.matrix( rbind( cbind(L.gg.t,L.gb.t),

                                        cbind(t(L.gb.t),L.bb.t ) ))
      }
      PQ32g
    })
    return(PQ32g)
  }

  # P/3+ Q/2 beta matrix
  P_Q_3_2_beta <- function(pars, fit = NULL)
  {
    if(is.null(fit))
    {
      fit <- key_quantities(pars)
    }
    PQ32b <- with(fit,{
      # parts constant in the matrix
      PQ32b <- array(data = NA, dim= c(nvarx + nvarz,nvarx + nvarz,nvarz))

      for(s in 1:nvarz)
      {
        zs <- z[,s]
        L.gg.s <- crossprod( x ,(1 / 3 * d1v * d.1.eta1 ^ 2 * d.1.eta2 / v ^ 2 +
                                   1 / 2 * (d.2.eta1 * d.1.eta2 / v - d1v *
                                              d.1.eta1 ^ 2 * d.1.eta2 / v ^ 2)) * zs * x)
        L.gb.s <- crossprod( x ,(1 / 3 * d1v * d.1.eta1 * d.1.eta2 ^ 2 / v ^ 2 +
                                   1 / 2 * (d.2.eta12 * d.1.eta2 / v - d1v *
                                              d.1.eta1 * d.1.eta2 ^ 2 / v ^ 2)) * zs * z)
        L.bb.s <- crossprod( z ,(1 / 3 * d1v * d.1.eta2 ^ 3 / v ^ 2 +
                                   1 / 2 * (d.2.eta2 * d.1.eta2 / v - d1v *
                                              d.1.eta2 ^ 3 / v ^ 2)) * zs * z)

        PQ32b[,,s] <- as.matrix( rbind( cbind(L.gg.s,L.gb.s),
                                        cbind(t(L.gb.s),L.bb.s ) ))
      }
      PQ32b
    })
    return(PQ32b)

  }
  # Matrix F1
  F1 <- function(pars, fit = NULL)
  {
    if(is.null(fit))
    {
      fit <- key_quantities(pars)
    }
    f1 <- rep(NA,nvarx + nvarz)
    p.q.gamma <- P_Q_gamma(fit=fit)
    p.q.beta <- P_Q_beta(fit=fit)
    info.inv <- information(inverse = TRUE, fit = fit)
    for(k in 1:(nvarx+nvarz))
    {
      if( k<=nvarx)
      {
        f1[k] <- sum(diag(info.inv %*% p.q.gamma[,, k])) # maybe change %*%
      }
      else
      {
        f1[k] <- sum(diag(info.inv %*% p.q.beta[,, k-nvarx] ))
      }

    }
    return(f1)
  }

  # h function
  h_function <- function(pars, fit = NULL)
  {
    if(is.null(fit))
    {
      fit <- key_quantities(pars)
    }
    h <- with(fit,{
      h <- array(data = NA, dim= c(nvarx + nvarz, nvarx + nvarz, nvarx + nvarz))
      info.inv <- information(inverse= TRUE, fit=fit)
      for( k in 1:(nvarx + nvarz) )
      {
        h[, , k] <-  tcrossprod( info.inv[, k], info.inv[, k]) / info.inv[k, k]
      }
      h
    })
    return(h)
  }

  # Creating F2 a function that calculate all the 2*(p-1) F2 vectors


  # Creating F2 a function that calculate all the 2*(p-1) F2 vectors
  F2_tilde <- function(pars, fit = NULL)
  {
    if(is.null(fit))
    {
      fit <- key_quantities(pars)
    }
    f2 <- matrix(NA,ncol = nvarx + nvarz, nrow = nvarx + nvarz )
    f2.tilde <- rep(NA, nvarx + nvarz)
    p.q.3.2.gamma <- P_Q_3_2_gamma(fit = fit)
    p.q.3.2.beta <- P_Q_3_2_beta(fit = fit)
    h<- h_function(fit=fit)
    inv.info <- information(inverse = TRUE, fit = fit)
    for(k in 1:(nvarx + nvarz))
    {
      if(k <= nvarx){
        for(j in 1:(nvarx + nvarz)){
          f2[j, k] <- sum(diag(h[, , j] %*% p.q.3.2.gamma[, , k]))
        }
      }
      else
      {
        for(j in 1:(nvarx + nvarz))
        {
          f2[j,k] <- sum(diag(h[,,j]%*%p.q.3.2.beta[,, k - nvarx ]))
        }
      }

    }
    for(l in 1:(nvarx + nvarz))
    {
      f2.tilde[l] <- inv.info[l, ]%*%f2[l,]
    }
    return(f2.tilde)
  }
  ## median  adjustment term


  #Estimation
  max.step.factor<- 12
  epsilon <- 1e-5
  par <- start ## initial value
  # Control for mle, mean bias reduction and median
  # bias reduction
  controlA <- 0
  controlB <- 0
  # Mean bias reduction
  if(method=="MeanBr")
  {
    controlA <- 1
  }
  # Median bias reduction
  if(method=="MedianBr")
  {
    controlA <- 1
    controlB <- 1
  }
  fit <- key_quantities(par)
  adjusted.grad <- (score(fit=fit) + 1/2 * F1(fit=fit) * controlA -
                      information(fit=fit) %*%  F2_tilde(fit=fit) * controlB)
  step.par <- information(fit=fit, inverse = TRUE) %*% adjusted.grad
  for (iter in seq.int(maxit))
  {
    step.factor <- 0
    testhalf <- TRUE
    while (testhalf & step.factor < max.step.factor)
    {
      step.par.previous <- step.par
      par <- par + slowit * 2 ^ ( - step.factor) * step.par
      fit <- key_quantities(par)
      adjusted.grad <- (score(fit=fit) + 1 / 2 * F1(fit=fit) * controlA -
                          information(fit=fit) %*% F2_tilde(fit=fit) * controlB)
      step.par <- information(fit=fit, inverse = TRUE) %*% adjusted.grad
      if (step.factor == 0 & iter == 1)
      {
        testhalf <- TRUE
      }
      else
      {
        testhalf <- sum(abs(step.par),na.rm=TRUE) > sum(abs(step.par.previous),na.rm = TRUE)
        # condition not very precise
      }
      step.factor <- step.factor + 1
    }

    if ( (sum(abs(adjusted.grad),na.rm = TRUE) < epsilon) | ( any(abs(par)>15) ) ) {
      break # was (sum(abs(step.par),na.rm = TRUE) < epsilon)
    }
  }


  if(iter >= maxit| (any(abs(par) > 15)))
  {
    convergence <- 0
    warning("optimization failed to converge")
  }
  else
  {
    convergence <- 1
  }
  Point.est <- par
  var.est <- information(Point.est, inverse = TRUE)
  se <- sqrt(diag(var.est))
  return(list(Point.est=Point.est, se=se,convergence=convergence))

}