README

The ordgam package enables to fit an additive proportional odds model to ordinal data.

The combination of Laplace approximations and of Bayesian P-splines (named LPS) enable fast and flexible inference in a Bayesian framework. The Gaussian Markov random field prior assumed for the penalized spline parameters and the Bernstein-von Mises theorem typically provide reliable Laplace approximation to the posterior distribution of these quantities. However, this accuracy can be seriously compromised for some unpenalized parameters, especially when the information synthesized by the prior and the likelihood is sparse.

We propose a refined version of the LPS methodology by splitting the parameter space in two subsets, \(\pmb{\zeta}=(\pmb{\gamma}^\top,\pmb{\theta}^\top)^\top\) . The first set involves the non-penalized parameters \(\pmb{\gamma}\) for which the joint posterior distribution \(p(\pmb{\gamma}|\pmb{\lambda},{\cal D})\) is approached from a non-Gaussian perspective with an approximation scheme tailored to capture asymmetry and kurtosis, while the conditional posterior distribution \(p(\pmb{\theta}|\pmb{\gamma},\pmb{\lambda},{\cal D})\) for the penalized parameters \(\pmb{\theta}\) in the complementary set undergoes the LPS treatment with Laplace approximations, see Fig.1 (with \(\pmb{\eta}=\pmb{\lambda}\)).

The marginal posterior density \(p(\pmb{\nu}|{\cal D})\) for the log of the penalty parameters, \(\pmb{\nu}=\log(\pmb{\lambda})\), also could be approximated using a product of skew-t densities. More details can be found in Lambert & Gressani (2023) [1].

The ordgam package in action

Let us illustrate the use of the ordgam package on a data subset (\(n=552\)) from the European Social Survey (ESS 2018) specific to the French speaking respondents from Wallonia, one of the three regions in Belgium. Each of the participants (aged at least 15) was asked to react to the following statement, Gay men and lesbians should be free to live their own life as they wish, with a positioning on a Likert scale going from 1 (Agree strongly) to 5 (Disagree strongly).

## Package installation and loading
## install.packages("devtools")
## devtools::install_github("plambertULiege/ordgam")
library(ordgam) 
    
## Data reading
data(freehmsDataBE)
donnees = subset(freehmsDataBE,region=="WAL") ## Focus on Wallonia
head(donnees)

Model fitting

Let us fit a proportional odds model to these data with the number of completed years of education (\(14.1\pm 4.4\) years) and age (\(47.3\pm 18.5\) years) entering as additive terms, \(L=10\) recentered B-splines spanning each covariate range, and Gamma priors for the penalty parameters, \(\lambda_j\sim{\cal G}(1,10^{-4})\) (\(j=1,2\)).

## Model fit with a Gamma prior for the penalty parameters
mod = ordgam(freehms ~ s(eduyrs) + s(age), data=donnees, descending = FALSE) 
print(mod)

plot(mod, mfrow=c(1,2))

Asymmetric posterior for the non-penalized parameters

The asymmetry of the posterior for the non-penalized parameters \(\pmb{\gamma}\) can be visualized. The first step in the approach requires the projection of \(\pmb{\gamma}\) on the eigenvectors of the variance-covariance matrix \(\Sigma_\lambda^{\gamma\gamma}\), yielding \(\tilde{\pmb{\gamma}}\) and its approximately independent components. The posterior density \(p(\tilde\gamma_k|\lambda,{\cal D})\) is further approximated using a skew-t density, revealing a non-negligible asymmetry for the posterior of \(\tilde\gamma_1\) (corresponding to the direction of the eigenvector with the largest eigenvalue):

model = mod  ## Model for which the approximation is required
    
## Skew-t approximation to the marginal of gamma.tilde[k]
ngamma = with(model, nalpha+nfixed)  ## Number of non-penalized parms
gamt.ST = list()  ## Skew-t coefs --> dst(x,dp=coef)
for (k in 1:ngamma){ ## Loop over the gamma.tilde components
     x.grid = seq(-4,4,length=10)  ## Grid of values for gamma.tilde[k]
     lfy.grid = ordgam::lmarg.gammaTilde(x.grid,k=k,model)  ## log p(gamma.tilde[k] | D)
     gamt.ST[[k]] = ordgam:::STapprox(x.grid,lfy.grid)$dp  ## Approximate using a skew-t
}
    
## Visualization of the posterior for <gamma.tilde[k]>
par(mfrow=c(2,2),mar=c(4,5,1,1))
for (k in 1:ngamma){ ## Loop over the gamma.tilde components
     xlab = bquote(tilde(gamma)[.(k)])
     ylab = bquote(p(tilde(gamma)[.(k)]~ "|"~lambda~","~D))
     xlim = sn::qst(c(.0001,.9999),dp=gamt.ST[[k]])
     curve(sn::dst(x,dp=gamt.ST[[k]]), xlim=xlim,
           xlab=xlab,ylab=ylab,col="blue",lwd=2,lty=1) 
}

The results can be re-expressed in the original parametrization, yieding \(p(\gamma_k|\lambda,{\cal D})\) and a noticable asymmetry for the posterior density of \(\gamma_4\):

## Sampling <gamma.tilde> using the (approximate) independence of its components
M = 10000
sample.tilde = matrix(nrow=M,ncol=ngamma)
for (k in 1:ngamma){
    ## ... or using skew-t approximation
    sample.tilde[,k] = c(sn::rst(M, dp=gamt.ST[[k]]))
}
    
## Revert sample to the original <gamma> parametrization
gam.hat = model$theta[1:ngamma]  ## MAP estimate
Sig = model$Sigma.theta[1:ngamma,1:ngamma] ## Variance-covariance
sv = svd(Sig) ; V = sv$u ; omega = sv$d  ## SVD
sample.gam = t(gam.hat + V %*% (sqrt(omega) * t(sample.tilde)))
    
## Skew-t approximation to the marginal of gam[k]
gam.ST = gam.ST2 = list()
for (k in 1:ngamma){
    temp = hist(sample.gam[,k],plot=FALSE)
    gam.ST[[k]] = ordgam:::STapprox(temp$mids,log(temp$density+1e-6))$dp
    ## temp = selm(sample.gam[,k] ~ 1, family="ST")
    ## gam.ST2[[k]] = coef(temp,"DP")
}
    
## Visualize p(gamma[k] | lambda,data)
par(mfrow=c(2,2),mar=c(4,5,1,1))
for (k in 1:ngamma){
    xlab = bquote(gamma[.(k)])
    ylab = bquote(p(gamma[.(k)]~ "|"~lambda~","~D))
    xlim = sn::qst(c(.0001,.9999),dp=gam.ST[[k]])
    curve(sn::dst(x,dp=gam.ST[[k]]), lwd=2,col="blue",
          xlim=xlim, xlab=xlab, ylab=ylab)
}

Marginal posterior of the penalty parameters

Let us revisit the same model for the data collected in Flanders, the most populated region of Belgium.

## Model fit with a Gamma prior for the penalty parameters
mod.FL = ordgam(freehms ~ s(eduyrs) + s(age), data=subset(freehmsDataBE,region=="FL"), 
                descending = FALSE) 
print(mod.FL)

fhat.FL = plot(mod.FL, mfrow=c(1,2))

The estimated additive terms for education and age look a bit different from the estimates for Wallonia with, here, the suggestion of a statistically significant and linear effect of ‘eduyrs’.

The marginal posterior distribution of the log of the penalty parameters, \(\nu_j=\log\lambda_j\), can be visualized and compared to the their prior:

model = mod.FL
## log(lambda) marginal posteriors
par(mfrow=c(1,2),mar=c(4,5,3,1))
for (j in 1:2){ ## Loop over <nu> components
  ## Plot of the approximating skew-t to the marginal posterior of <nu[j]>
  xlims = sn::qst(c(.001,.999),dp=model$nu.dp[[j]])
  xlab = bquote(log(lambda[.(j)]))
  ylab = bquote(p(log(lambda[.(j)])~ "|"~D))
  curve(sn::dst(x,dp=model$nu.dp[[j]]),xlims[1],xlims[2],lwd=2,col="red",
        main=names(model$lambda)[j],xlab=xlab,ylab=ylab,ylim=c(0,.4))
  ## ... compared to the prior for <nu[j]>
  curve(exp(model$lprior.lambda(exp(x)))*exp(x),add=T,lty=3,lwd=2)
  legend("topright", legend=c("Posterior","Prior"),
         lwd=c(2,2), lty=c(1,3), col=c("red","black"),
         horiz=TRUE, bty="n",cex=1)
}

While the first posterior associated to penalty parameter for ‘eduyrs’ does not depart from its prior with a large prior mean (as \(\lambda_1 \sim {\cal G}(1,10^{-4})\) has prior mean \(10^4\)) and suggests linearity, the posterior associated to \(\nu_2=\log(\lambda_2)\) suggests a non-linear effect of ‘age’ by pointing to smaller values for the penalty parameter.

License

This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/.

References

[1] Lambert, P. and Gressani, O (2023). Penalty parameter selection and asymmetry corrections to Laplace approximations in Bayesian P-splines models. Statistical Modelling doi:10.1177/1471082X231181173. Preprint: arXiv:2210.01668

The ordgam R-package

The ordgam package in action

Model fitting

Asymmetric posterior for the non-penalized parameters

Marginal posterior of the penalty parameters

License

References