diff --git a/R/utility.R b/R/utility.R
index e08b0ca..481d8a5 100644
--- a/R/utility.R
+++ b/R/utility.R
@@ -7,6 +7,9 @@
 #' @param object Output from \code{\link{dsem}}
 #' @param nsim Number of simulations to use if \code{family!="fixed"} for some variable,
 #'        such that simulation residuals are required.
+#' @param what whether to return quantile residuals, or samples from the leave-one-out predictive
+#'        distribution of data, or a table of leave-one-out predictions and standard errors for the
+#'        latent state
 #' @param ... Not used
 #'
 #' @details
@@ -27,9 +30,11 @@
 loo_residuals <-
 function( object,
           nsim = 100,
+          what = c("quantiles","samples","loo"),
           ... ){
 
   # Extract and make object
+  what = match.arg(what)
   tsdata = object$internal$tsdata
   parlist = object$obj$env$parList()
   df = expand.grid( time(tsdata), colnames(tsdata) )
@@ -84,13 +89,18 @@ function( object,
   }
 
   # Compute quantile residuals
+  resid_tjr = array( NA, dim=c(dim(tsdata),nsim) )
   if( all(object$internal$family == "fixed") ){
+    # analytical quantile residuals
     resid_tj = array( NA, dim=dim(tsdata) )
     resid_tj[cbind(df[,1],df[,2])] = pnorm( q=df$obs, mean=df$est, sd=df$se )
+    # Simulations from predictive distribution for use in DHARMa etc
+    for(z in 1:nrow(df) ){
+      resid_tjr[df[z,1],df[z,2],] = rnorm( n=nsim, mean=df[z,'est'], sd=df[z,'se'] )
+    }
   }else{
     # Sample from leave-one-out predictive distribution for states
     resid_rz = apply( df, MARGIN=1, FUN=\(vec){rnorm(n=nsim, mean=as.numeric(vec['est']), sd=as.numeric(vec['se']))} )
-    resid_tjr = array( NA, dim=c(dim(tsdata),nsim) )
     # Sample from predictive distribution of data given states
     for(r in 1:nrow(resid_rz) ){
       parameters = object$obj$env$parList()
@@ -107,5 +117,7 @@ function( object,
     # Calculate quantile residuals
     resid_tj = apply( resid_tjr > outer(tsdata,rep(1,nsim)), MARGIN=1:2, FUN=mean )
   }
-  return( resid_tj )
+  if(what=="quantiles") return( resid_tj )
+  if(what=="samples") return( resid_tjr )
+  if(what=="loo") return( df )
 }
diff --git a/man/loo_residuals.Rd b/man/loo_residuals.Rd
index ac55a80..260891e 100644
--- a/man/loo_residuals.Rd
+++ b/man/loo_residuals.Rd
@@ -4,7 +4,7 @@
 \alias{loo_residuals}
 \title{Calculate leave-one-out residuals}
 \usage{
-loo_residuals(object, nsim = 100, ...)
+loo_residuals(object, nsim = 100, what = c("quantiles", "samples", "loo"), ...)
 }
 \arguments{
 \item{object}{Output from \code{\link{dsem}}}
@@ -12,6 +12,10 @@ loo_residuals(object, nsim = 100, ...)
 \item{nsim}{Number of simulations to use if \code{family!="fixed"} for some variable,
 such that simulation residuals are required.}
 
+\item{what}{whether to return quantile residuals, or samples from the leave-one-out predictive
+distribution of data, or a table of leave-one-out predictions and standard errors for the
+latent state}
+
 \item{...}{Not used}
 }
 \value{