diff --git a/R/utility.R b/R/utility.R
index 481d8a5..bd7ea99 100644
--- a/R/utility.R
+++ b/R/utility.R
@@ -45,7 +45,7 @@ function( object,
   # Loop through observations
   for(r in 1:nrow(df) ){
     ts_r = tsdata
-    ts_r[df[r,1],df[r,2]] = NA
+    ts_r[match(df[r,1],time(tsdata)),match(df[r,2],colnames(tsdata))] = NA
 
     #
     control = object$internal$control
@@ -84,8 +84,8 @@ function( object,
                   objective = obj$fn,
                   gradient = obj$gr )
     sdrep = TMB::sdreport( obj )
-    df[r,'est'] = as.list(sdrep, what="Estimate")$x_tj[df[r,1],df[r,2]]
-    df[r,'se'] = as.list(sdrep, what="Std. Error")$x_tj[df[r,1],df[r,2]]
+    df[r,'est'] = as.list(sdrep, what="Estimate")$x_tj[match(df[r,1],time(tsdata)),match(df[r,2],colnames(tsdata))]
+    df[r,'se'] = as.list(sdrep, what="Std. Error")$x_tj[match(df[r,1],time(tsdata)),match(df[r,2],colnames(tsdata))]
   }
 
   # Compute quantile residuals
@@ -93,10 +93,10 @@ function( object,
   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 )
+    resid_tj[cbind(match(df[,1],time(tsdata)),match(df[,2],colnames(tsdata)))] = 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'] )
+      resid_tjr[match(df[z,1],time(tsdata)),match(df[z,2],colnames(tsdata)),] = rnorm( n=nsim, mean=df[z,'est'], sd=df[z,'se'] )
     }
   }else{
     # Sample from leave-one-out predictive distribution for states
diff --git a/vignettes/vignette.Rmd b/vignettes/vignette.Rmd
index a1a3842..2abc1f8 100644
--- a/vignettes/vignette.Rmd
+++ b/vignettes/vignette.Rmd
@@ -59,6 +59,21 @@ knitr::kable( m1, digits=3)
 
 This shows that linear and dynamic structural equation models give identical estimates of the single path coefficient.
 
+We can also calculate leave-one-out residuals and display them using DHARMa:
+
+```{r, echo=TRUE, message=FALSE, fig.width=7, fig.height=7}
+# sample-based quantile residuals displayed using DHARMa
+samples = loo_residuals(fit, what="samples")
+which_use = which(!is.na(data))
+fittedPredictedResponse = loo_residuals(fit, what="loo")
+res = DHARMa::createDHARMa( simulatedResponse = apply(samples,MARGIN=3,FUN=as.vector)[which_use,],
+                            observedResponse = unlist(data)[which_use],
+                            fittedPredictedResponse = fittedPredictedResponse$est )
+plot(res)
+```
+
+This shows no pattern in the quantile-quantile plot, as expected given that we have a correctly specified model
+
 ## Comparison with dynamic linear models
 
 We first demonstrate that `dsem` gives identical results to `dynlm` for a well-known econometric model, the Klein-1 model.