Proportional hazard model for ordered response

ordered_ph.R

library(dplyr)
library(ggplot2)
library(tidyr)
library(tibble)
library(MASS)

make_prob <- function(alpha, eta){
  xi <- plogis(outer(eta, alpha, "+"))
  for(j in 2:length(alpha)){
    xi[,j] <- xi[,j]*(1-rowSums(xi[,1:(j-1), drop=FALSE]))  
  }
  cbind(xi, 1-rowSums(xi))  
}

PHreg <-function(Y,X,prec){
  K <- ncol(X)
  J <- ncol(Y)
  M <- t(apply(Y,1, function(x)rev(cumsum(rev(x)))))
  LL <- function(par,Y,M,J,X,K,prec){
    alpha <- par[1:(J-1)]
    beta <- par[J:(J+K-1)]
    out <- 0
    Xbeta <- -X%*%beta 
    for(j in 1:(J-1)){
      out <- out +
        sum(Y[,j]*(alpha[j]+Xbeta) -
              M[,j]*log1p(exp(alpha[j]+Xbeta)))
    }
    out <- out - 0.5*sum(beta^2*prec)
    return(out)
  }
  
  dLL <- function(par,Y,M,J,X,K,prec){
    alpha <- par[1:(J-1)]
    beta <- par[J:(J+K-1)]
    dLLalpha <- function(alpha,beta,Y,M,J,X){
      Xbeta <- -X%*%beta
      out <- numeric(J-1)
      for(j in 1:(J-1)){
        out[j] <- sum(Y[,j] - M[,j]*plogis(alpha[j]+Xbeta))
      }
      return(out)
    }
    
    dLLbeta <- function(alpha,beta,Y,M,J,X,K,prec){
      Xbeta <- -drop(X%*%beta)
      out <- numeric(K)
      for(j in 1:(J-1)){
        out <- out + colSums(-Y[,j]*X + M[,j]*X*plogis(alpha[j]+Xbeta))
      }
      out <- out - beta*prec
      return(out)
    }
    
    return(
      c(dLLalpha(alpha, beta, Y, M, J, X),
        dLLbeta(alpha, beta, Y, M, J, X, K, prec))
    )
  }
  ini <-numeric(J+K-1)
  xnames <- colnames(X)
  if(is.null(xnames)){
    xnames <- paste0("X", seq_len(ncol(X)))
  }
  names(ini)<-c(paste0("alpha",1:(J-1)), xnames)
  opt <- optim(ini, LL, dLL, method="BFGS",
               control = list(fnscale=-1),
               J=J,Y=Y,M=M,X=X,K=K,prec=prec,hessian = TRUE)
  return(opt)
}



fit_l <- polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing, 
              method = "logistic", Hess = TRUE)

CI_l <- confint(fit_l)
colnames(CI_l) <- c("lower", "upper")
CI_l <- rownames_to_column(
    data.frame(CI_l, est=fit_l$coefficients, method="logistic")
  )

fit_p <- polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing, 
              method = "probit", Hess = TRUE)
CI_p <- confint(fit_p)
colnames(CI_p) <- c("lower", "upper")
CI_p <- rownames_to_column(
    data.frame(CI_p, est=fit_p$coefficients, method="probit")
  )

housing_wide <- pivot_wider(housing, names_from = Sat, values_from = Freq)
Y <- as.matrix(housing_wide[,4:6])
X <- model.matrix(~ Infl+Type+Cont, housing_wide)
X <- X[,-1]
res <- PHreg(Y,X,0)
se <- sqrt(-diag(solve(res$hessian)))
q <- qnorm(0.975)
CI_ph <- data.frame(lower=res$par - q*se,
                    upper=res$par + q*se,
                    est=res$par, method="PH")
CI_ph <- rownames_to_column(CI_ph[3:8,])

CI <- bind_rows(CI_p, CI_l, CI_ph)
ggplot(CI, aes(x= rowname, y=est, colour=method))+
  geom_pointrange(aes(ymin=lower,ymax=upper), position = position_dodge(width = 0.5))+
  geom_hline(yintercept = 0, linetype=2)+
  coord_flip()+
  theme_bw(12)