Nothing
#' @title Ordinal regression for continuous scales
#' @description Continuous ordinal regression with logit link using I-splines to model the g function.
#' @param formula a formula expression as for regression models, of the form
#' response ~ predictors. Only fixed effects are supported.
#' The model must have an intercept: attempts to remove one will lead to a warning and will be
#' ignored.
#' @param data an optional data frame in which to interpret the variables occurring in the
#' formulas
#' @param scale a vector of length 2 with the boundaries of the ordinal scale used. If not specified, the range of the data is used, and a warning is displayed.
#' @param weights optional case weights in fitting. Defaults to 1.
#' @param link link function, i.e. the type of location-scale distribution assumed for
#' the latent distribution. The default ``logit'' link gives the proportional odds model.
#' Other options are "logit", "probit", "cloglog", "loglog", "cauchit".
#' @param niters a vector of length 2 with the maximimum number of external and internal
#' iterations used in the fitting algorithm. The internal algorithm estimates the parameters
#' of the model conditional on the current values of \eqn{\lambda}s, the smoothing parameters.
#' The external algorithm estimates the values of \eqn{\lambda}s conditional on the current
#' estimates of the parameters of the model. Default is \code{c(500,500)}
#' @param conv_crit the smoothing parameters \eqn{\lambda}'s convergence criteria for the iterative process.
#' Default is \eqn{0.01}
#' @param n.int.knots the number of internal knots used to compute the spline bases. The default (NULL) is round((n-1-order)*0.8) if in the interval [8,15], and 8 or 15 otherwise.
#' @param order the order of the spline functions. The default is 4 (cubic splines).
#' @param lambdas NA (the default) or a vector of length equal to the number of smoothing terms, including the g function and, optionally, the random effect terms and the smooters. If ``lambdas'' is a vector, each element \eqn{\lambda_i} can be a number, in which case the corresponding term is penalized using \eqn{\lambda_i} as smoothing parameter, zero, in which case the corresponding term is unpenalized, or NA, in which case the value of \eqn{\lambda_i} is estimated maximmizing the marginal posterior function.
#' @concept likelihood log-likelihood ordinal regression
#' @details Fits a continuous ordinal regression model using penalized maximum likelihood.
#' The model can contain fixed effects and optionally mixed effects and smoothers.
#' The g function is estimated using monotone increasing I-splines, and the link function is the logit,
#' implying the standard logistic distribution for the latent variable. Penalized maximum likelihood
#' estimation is performed using the \code{MI} algorithm and the splines smoothing parameters are estimated
#' maximizing the marginal posterior (details of the iterative process are printed out during the fit).
#'
#' @return an object of type \code{ocm} with the components listed below. Parameter estimates are in \code{coefficients}.
#' \item{coefficients}{parameter estimates}
#' \item{pars_obj}{an object of class \code{ocmpars} carrying the parameter estimates and other properties of the regression terms}
#' \item{vcov}{variance-covariance matrix}
#' \item{H}{the Hessian matrix}
#' \item{logLik}{value of the log-likelihood at the estimated optimum}
#' \item{penlogLik}{value of the lenalized log-likelihood at the estimated optimum}
#' \item{v}{vector of continuous scores}
#' \item{sample.size}{sample size (can differ from the number of observations if the weights are different from 1)}
#' \item{edf}{estimated degrees of freedom}
#' \item{df.residual}{the residual degrees of freedom}
#' \item{nobs}{number of observations}
#' \item{terms}{model terms}
#' \item{call}{call to fit the model}
#' \item{data}{the data frame as in input, ordered by the outcome values}
#' \item{model.frame}{the model.frame used in the fit}
#' \item{model.matrix}{the model.matrix used in the fit}
#' \item{weights}{case weights in fitting}
#' \item{sorting}{the ordinal score v sorting vector}
#' \item{link}{link function used}
#' \item{formula}{formula used}
#' \item{scale}{the boundaries of the ordinal scale used}
#' @references Manuguerra M, Heller GZ (2010). Ordinal Regression Models for Continuous
#' Scales, \emph{The International Journal of Biostatistics}: 6(1), Article 14.
#' @references Manuguerra M, Heller GZ, Ma J (2017). Semi-parametric Ordinal Regression Models for Continuous
#' Scales, \emph{Proceedings of the 32nd International Workshop on Statistical Modelling}. July 3-7, 2017, Groningen, Netherlands.
#' @references Manuguerra M, Heller GZ, Ma J (2020). Continuous Ordinal Regression for Analysis of Visual Analogue Scales: The R Package ordinalCont, \emph{Journal of Statistical Software}. 96(8). doi:10.18637/jss.v096.i08
#' @author Maurizio Manuguerra, Gillian Heller
#' @import grDevices
#' @import graphics
#' @import stats
#' @import utils
#' @export
#' @examples
#' fit.overall <- ocm(overall ~ cycleno + age + bsa + treatment, data=ANZ0001.sub, scale=c(0,100))
#' summary(fit.overall)
#' \dontrun{
#' plot(fit.overall)
#' ## Smoothers and complete data set
#' fit.overall.smooth <- ocm(overall ~ age + treatment : s(cycleno), data=ANZ0001, scale=c(0,100))
#' summary(fit.overall.smooth)
#' plot(fit.overall.smooth)
#' }
ocm <- function(formula, data=NULL, scale=NULL, weights, link = c("logit", "probit", "cloglog", "loglog", "cauchit"), niters=c(500,500), conv_crit=1e-2, n.int.knots=NULL, order=4, lambdas=NA)
{
lambda=0
if (missing(formula) | length(formula)<3)
stop("Model needs a formula")
link <- match.arg(link)
if(is.null(data)) data <- model.frame(formula=formula, data=parent.frame(n=1))
if(missing(weights)) weights <- rep(1, nrow(data))
keep <- weights > 0
data <- data[keep,]
weights <- weights[keep]
data=cbind(Intercept=rep(1,nrow(data)), data)
#Model Response
v = eval(formula[[2]], envir = data)
v <- as.numeric(v)
if (is.null(scale)) {scale=range(v); warning(paste("The range of the data [", min(v),", ",max(v),"] is used to scale data to the interval (0,1).\n", sep=''))}
if (min(v)<scale[1] | max(v)>scale[2]) {stop(paste("Ordinal scores (v) outside the ",scale[1],"-",scale[2]," range have been found. Please either rescale or use the 'scale' parameter when calling ocm.",sep=''))}
v <- (v-scale[1])/(scale[2]-scale[1])
v <- ((length(v)-1)*v+0.5)/length(v)
n <- length(unique(v))
### Sort data set by v #FIXME: do we need this?
sorting_ind <- order(v)
v <- v[sorting_ind]
data <- data[sorting_ind,]
weights <- weights[sorting_ind]
### Create ocmPARS object
pars_obj=ocmPars(formula, data, v, n.int.knots, order)
### Terms
terms_mf_mm=ocmTerms(formula, data)
###
pen_index = which(sapply(pars_obj, function(x)x$use_lambda))
n_lambda = length(pen_index)
if (any(!is.na(lambdas))){
il = which(!sapply(lambdas,is.na))
ipen = pen_index[il]
lambdas = lambdas[il]
for (ii in 1:length(ipen)) {pars_obj[[ipen[ii]]]$estimate_lambda=F; pars_obj[[ipen[ii]]]$lambda=lambdas[ii]}
}
#Fit
regression_edf0 <- sum(sapply(pars_obj, function(x)ifelse(x$type=="fix.eff",x$len,0)))
pen_index = which(sapply(pars_obj, function(x)x$estimate_lambda))
##pen_index at least =1 as g function is now non-parametric
if (length(pen_index)>0) cat("Ext.iters\tInt.iters\tlambda",rep("\t",2*length(pen_index)),"Convergence (<",conv_crit,")\n", sep='')
for (iter in 1:niters[1]){
convergence <- NULL
conv_val <- NULL
regression_edf <- regression_edf0
##
if (link == "logit" | link=="probit" | link=="cloglog" | link=="loglog" | link=="cauchit"){
deriv_funs <- deriv_link(link=link)
} else {
stop("link function not implemented.")
}
est <- ocmEst4(v, weights, pars_obj, deriv_funs, niters[2], conv_crit=conv_crit)
pars_obj <- est[["pars_obj"]]
Ginv = est$vcov
if (length(pen_index)>0){
for (ipen in pen_index){
oo = pars_obj[[ipen]]
lambda_old = oo$lambda
sigma2_old = 1/(2*lambda_old)
Q <- oo$Rstar/sigma2_old
edf <- oo$len-(sum(diag(Ginv%*%Q)))
sigma2 = c(t(oo$pars)%*%oo$R%*%oo$pars/edf)
lambda_old = ifelse(iter>1,lambda_old,-1)
lambda <- pars_obj[[ipen]]$lambda <- 1/(2*sigma2)
conv_val <- c(conv_val, abs(lambda-lambda_old)/abs(lambda_old))
convergence = c(convergence, (abs(lambda-lambda_old)/abs(lambda_old))<conv_crit)
regression_edf <- regression_edf + edf
}
# check for convergence
#cat(iter,"\t\t",est$iter,"\t\t",conv_val,"\n")
cat(iter,"\t\t",est$iter,"\t\t",sapply(pars_obj[pen_index],function(x)x$lambda),"\t\t",round(conv_val,4),"\t\t","\n")
if(all(convergence)){break}
} else {
break
}
}
if (iter==niters[1]) {
warn_msg="The process did not converge. Try increasing the number of iterations (eg niters=c(500,1000))."
warning(warn_msg)
}
est$v <- v
est$sample.size <- nrow(data)
est$edf <- regression_edf
est$df.residual <- sum(weights)-regression_edf
est$nobs <- sum(weights)
est$terms <- terms_mf_mm$terms
est$call <- match.call()
est$data <- data
est$model.frame <- terms_mf_mm$mf
est$model.matrix <- terms_mf_mm$mm
est$weights <- weights
est$sorting <- sorting_ind
est$link <- link
est$formula <- formula
est$scale <- scale
est$iter <- NULL
class(est) <- "ocm"
est
}
ocmEst4 <- function(v, weights, pars_obj, deriv_funs, int_iters, conv_crit=1e-3){
n <- length(v)
start <- split_obj2pars(pars_obj)
#fit0 <- NewtonMi(pars_obj, maxiters=int_iters, wts=weights, convVal=conv_crit)
fit <- NewtonMi_gen(pars_obj, deriv_funs=deriv_funs, maxiters=int_iters, wts=weights, convVal=conv_crit)
coef <- fit$par
names(coef) <- names(start)
pars_obj <- split_pars2obj(pars_obj, coef)
# curval <- current.values(pars_obj, deriv_funs)
#H <- -hessian(pars_obj)
# H <- -hessian_gen(pars_obj, curval)
# G <- compute_G(H, pars_obj)
vcov <- solve(fit$hessian)
## degrees of freedom and standard deviation of residuals
fix.eff_index = which(sapply(pars_obj, function(x)x$type)=="fix.eff")
colnames(vcov) <- rownames(vcov) <- names(coef)
#logLik <- sum(weights * llik(pars_obj))
logLik <- sum(weights * llik_gen(pars_obj, fit$curval))
list(coefficients = coef,
pars_obj = pars_obj,
vcov = vcov,
H = fit$hessian,
iter=fit$iter,
logLik = logLik,
penlogLik = -fit$value)
}
penalty <- function(pars_obj){
pen_index = which(sapply(pars_obj, function(x)x$use_lambda))
if (length(pen_index)==0){
pen=0
} else {
pen=sapply(pen_index, function(i){with(pars_obj[[i]], lambda * t(pars) %*% (R %*% pars))})
}
if (any(pen<0)) {
print(pen)
print(sapply(pen_index, function(i){with(pars_obj[[i]], print(list(lambda,pars)))}))
warning("Something went wrong with the penalty term.")
pen[which(pen<0)] <- 0
}
return(pen)
}
lin_regressor <- function(pars_obj){
# if (!intercept){
# fix_index = which(sapply(pars_obj, function(x)x$type)=="fix.eff")
# pars_obj <- pars_obj[-fix_index]
# }
# if (length(pars_obj)==0) return(0)
h_comp = sapply(pars_obj, function(x)x$mat %*% x$pars)
apply(h_comp,1,sum)
}
CDF <- function(pars_obj){
h <- lin_regressor(pars_obj)
return(exp(h)/(1+exp(h)))
}
#Only loglik, no pen
hessian <- function(pars_obj){
gfun_index = which(sapply(pars_obj, function(x)x$type)=="gfun")
dg <- dgfun(pars_obj)
h <- lin_regressor(pars_obj)
F <- CDF(pars_obj)
nr <- sum(sapply(pars_obj,function(oi)oi$len))
H_mat <- matrix(NA, nrow=nr, ncol=nr)
H_mat2 <- matrix(NA, nrow=nr, ncol=nr)
rubric <- make_rubric(pars_obj)
D0 <- -2*exp(h)/(1+exp(h))^2
D <- diag(D0)
E0 <- as.numeric(-1/dg^2)
E <- diag(E0)
#FIXME optimise/compare H-mat vs H_mat2
for (i in 1:(nrow(rubric)-1)){
for (j in (i+1):nrow(rubric)){
t1=myTXAY(pars_obj[[i]]$mat, D0, pars_obj[[j]]$mat)
#t2=t(pars_obj[[i]]$mat) %*% D %*% pars_obj[[j]]$mat
##
H_mat[rubric[i,1]:rubric[i,2],rubric[j,1]:rubric[j,2]] <- t1
H_mat[rubric[j,1]:rubric[j,2],rubric[i,1]:rubric[i,2]] <- t(t1)
##
#H_mat2[rubric[i,1]:rubric[i,2],rubric[j,1]:rubric[j,2]] <- t2
#H_mat2[rubric[j,1]:rubric[j,2],rubric[i,1]:rubric[i,2]] <- t(t2)
}
}
for (i in 1:nrow(rubric)){
t1=myTXAY(pars_obj[[i]]$mat, D0, pars_obj[[i]]$mat)
#t2=t(pars_obj[[i]]$mat) %*% D %*% pars_obj[[i]]$mat
if (i==gfun_index) {
t1=t1+myTXAY(pars_obj[[i]]$mat1, E0, pars_obj[[i]]$mat1)
#t2=t2+t(pars_obj[[i]]$mat1) %*% E %*% pars_obj[[i]]$mat1
}
H_mat[rubric[i,1]:rubric[i,2],rubric[i,1]:rubric[i,2]] <- t1
#H_mat2[rubric[i,1]:rubric[i,2],rubric[i,1]:rubric[i,2]] <- t2
}
return(H_mat)
}
compute_G <- function(H, pars_obj){
pen_index = which(sapply(pars_obj, function(x)x$use_lambda))
if (length(pen_index)==0){
return(H)
} else {
Qs=lapply(pen_index, function(i){with(pars_obj[[i]], 2*lambda*Rstar)})
}
G <- H
for (Qi in Qs) {G <- G + Qi}
return(G)
}
compute_Rstar <- function(pars_obj){
#FIXME optimise
pen_index = which(sapply(pars_obj, function(x)x$use_lambda))
if (length(pen_index)>0){
Rs=lapply(pen_index, function(i){with(pars_obj[[i]], R)})
rubric <- make_rubric(pars_obj)
m = max(rubric)
Rstar0=matrix(0,nrow=m,ncol=m)
for (i in 1:length(pen_index)) {
inds=rubric[pen_index[i],1]:rubric[pen_index[i],2]
pars_obj[[pen_index[i]]]$Rstar = Rstar0
pars_obj[[pen_index[i]]]$Rstar[inds,inds] <- Rs[[i]]
}
}
return(pars_obj)
}
compute_gmat1_ext <- function(pars_obj){
rubric <- make_rubric(pars_obj)
gfun_index = which(sapply(pars_obj, function(x)x$type)=="gfun")
n = length(pars_obj[[gfun_index]]$v)
gfun_range = rubric[gfun_index,]
gfun_inds = gfun_range[1]:gfun_range[2]
nr <- sum(sapply(pars_obj,function(oi)oi$len))
gmat1_ext = matrix(0, nrow=n, ncol=nr)
# print(str(pars_obj[[gfun_index]]$mat2))
# print(gfun_inds)
gmat1_ext[, gfun_inds] = pars_obj[[gfun_index]]$mat1
pars_obj[[gfun_index]]$mat1_ext = gmat1_ext
return(pars_obj)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.