problem viewing summary for large number of params

[Marc Kéry]

Dear colleagues,

I am fitting a multi-species occupancy model with 104 species and 35 covariates in the occupancy part of the model. I can view on the screen the marginal posterior summaries for the community parameters, but not for the species level: on my console, the output of the latter (3780 params long) retains only the last couple 100 params  --- when I try to scroll up, I find the upper part of the summary is lost. 

Doing 'summary(out1) -> tmp' and then viewing 'tmp' selectively doesn't help (nothing gets assigned to 'tmp'). I can of course compute some or potentially all of the marginal posterior summaries for the species myself, but was wondering whether there was some ready solution or work-around.

Thank you and best regards  --- Marc

[Jeffrey Doser]

Hi Marc,

I unfortunately don't think I have a super great solution for you. You might be able to change some default printing behavior in R in order to prevent losing the upper part of the summary by changing some printing option or something like that, but I've never tried that and it may not be super easy to navigate with that many species/covariates. You could also modify your the summary.msPGOcc function to actually save the output instead of just printing it to the console. I pasted the code for the summary function below (if using one of the spatial models, you can find that code in the "R/generics.R" script in the spOccupancy github.

Cheers,

Jeff

summary.msPGOcc <- function(object,
   level = 'both',
   quantiles = c(0.025, 0.5, 0.975),
   digits = max(3L, getOption("digits") - 3L), ...) {

  print(object)

  n.post <- object$n.post
  n.samples <- object$n.samples
  n.burn <- object$n.burn
  n.thin <- object$n.thin
  n.chains <- object$n.chains
  run.time <- object$run.time[3] / 60 # minutes

  cat(paste("Samples per Chain: ", n.samples,"\n", sep=""))
  cat(paste("Burn-in: ", n.burn,"\n", sep=""))
  cat(paste("Thinning Rate: ",n.thin,"\n", sep=""))
  cat(paste("Number of Chains: ", n.chains, "\n", sep = ""))
  cat(paste("Total Posterior Samples: ",n.post * n.chains,"\n", sep=""))
  cat(paste("Run Time (min): ", round(run.time, digits), "\n\n", sep = ""))

  if (tolower(level) %in% c('community', 'both')) {

    cat("----------------------------------------\n");
    cat("\tCommunity Level\n");
    cat("----------------------------------------\n");

    # Occurrence
    cat("Occurrence Means (logit scale): \n")
    tmp.1 <- t(apply(object$beta.comm.samples, 2,
            function(x) c(mean(x), sd(x))))
    colnames(tmp.1) <- c("Mean", "SD")
    tmp <- t(apply(object$beta.comm.samples, 2,
          function(x) quantile(x, prob = quantiles)))
    diags <- matrix(c(object$rhat$beta.comm, round(object$ESS$beta.comm, 0)), ncol = 2)
    colnames(diags) <- c('Rhat', 'ESS')

    print(noquote(round(cbind(tmp.1, tmp, diags), digits)))

    cat("\nOccurrence Variances (logit scale): \n")
    tmp.1 <- t(apply(object$tau.sq.beta.samples, 2,
            function(x) c(mean(x), sd(x))))
    colnames(tmp.1) <- c("Mean", "SD")
    tmp <- t(apply(object$tau.sq.beta.samples, 2,
          function(x) quantile(x, prob = quantiles)))
    diags <- matrix(c(object$rhat$tau.sq.beta, round(object$ESS$tau.sq.beta, 0)), ncol = 2)
    colnames(diags) <- c('Rhat', 'ESS')
    print(noquote(round(cbind(tmp.1, tmp, diags), digits)))
    if (object$psiRE) {
      cat("\n")
      cat("Occurrence Random Effect Variances (logit scale): \n")
      tmp.1 <- t(apply(object$sigma.sq.psi.samples, 2,
              function(x) c(mean(x), sd(x))))
      colnames(tmp.1) <- c("Mean", "SD")
      tmp <- t(apply(object$sigma.sq.psi.samples, 2,
            function(x) quantile(x, prob = quantiles)))
      diags <- matrix(c(object$rhat$sigma.sq.psi, round(object$ESS$sigma.sq.psi, 0)), ncol = 2)
      colnames(diags) <- c('Rhat', 'ESS')

      print(noquote(round(cbind(tmp.1, tmp, diags), digits)))
    }
    cat("\n")
    # Detection
    cat("Detection Means (logit scale): \n")
    tmp.1 <- t(apply(object$alpha.comm.samples, 2,
            function(x) c(mean(x), sd(x))))
    colnames(tmp.1) <- c("Mean", "SD")
    tmp <- t(apply(object$alpha.comm.samples, 2,
          function(x) quantile(x, prob = quantiles)))
    diags <- matrix(c(object$rhat$alpha.comm, round(object$ESS$alpha.comm, 0)), ncol = 2)
    colnames(diags) <- c('Rhat', 'ESS')

    print(noquote(round(cbind(tmp.1, tmp, diags), digits)))
    cat("\nDetection Variances (logit scale): \n")
    tmp.1 <- t(apply(object$tau.sq.alpha.samples, 2,
            function(x) c(mean(x), sd(x))))
    colnames(tmp.1) <- c("Mean", "SD")
    tmp <- t(apply(object$tau.sq.alpha.samples, 2,
          function(x) quantile(x, prob = quantiles)))
    diags <- matrix(c(object$rhat$tau.sq.alpha, round(object$ESS$tau.sq.alpha, 0)), ncol = 2)
    colnames(diags) <- c('Rhat', 'ESS')

    print(noquote(round(cbind(tmp.1, tmp, diags), digits)))
    if (object$pRE) {
      cat("\n")
      cat("Detection Random Effect Variances (logit scale): \n")
      tmp.1 <- t(apply(object$sigma.sq.p.samples, 2,
              function(x) c(mean(x), sd(x))))
      colnames(tmp.1) <- c("Mean", "SD")
      tmp <- t(apply(object$sigma.sq.p.samples, 2,
            function(x) quantile(x, prob = quantiles)))
      diags <- matrix(c(object$rhat$sigma.sq.p, round(object$ESS$sigma.sq.p, 0)), ncol = 2)
      colnames(diags) <- c('Rhat', 'ESS')

      print(noquote(round(cbind(tmp.1, tmp, diags), digits)))
    }
  }

  if (tolower(level) %in% c('species', 'both')) {
    if (tolower(level) == 'both') cat("\n")
    cat("----------------------------------------\n");
    cat("\tSpecies Level\n");
    cat("----------------------------------------\n");
    cat("Occurrence (logit scale): \n")
    tmp.1 <- t(apply(object$beta.samples, 2,
            function(x) c(mean(x), sd(x))))
    colnames(tmp.1) <- c("Mean", "SD")
    tmp <- t(apply(object$beta.samples, 2,
          function(x) quantile(x, prob = quantiles)))
    diags <- matrix(c(object$rhat$beta, round(object$ESS$beta, 0)), ncol = 2)
    colnames(diags) <- c('Rhat', 'ESS')

    print(noquote(round(cbind(tmp.1, tmp, diags), digits)))
    cat("\n")
    # Detection
    cat("Detection (logit scale): \n")
    tmp.1 <- t(apply(object$alpha.samples, 2,
            function(x) c(mean(x), sd(x))))
    colnames(tmp.1) <- c("Mean", "SD")
    tmp <- t(apply(object$alpha.samples, 2,
          function(x) quantile(x, prob = quantiles)))
    diags <- matrix(c(object$rhat$alpha, round(object$ESS$alpha, 0)), ncol = 2)
    colnames(diags) <- c('Rhat', 'ESS')
    print(noquote(round(cbind(tmp.1, tmp, diags), digits)))

  }

  if (is(object, 'tMsPGOcc')) {
    if (object$ar1) {
      cat("\n")
      cat("Occurrence AR(1) Temporal Covariance: \n")
      tmp.1 <- t(apply(object$theta.samples, 2,
              function(x) c(mean(x), sd(x))))
      colnames(tmp.1) <- c("Mean", "SD")
      tmp <- t(apply(object$theta.samples, 2,
            function(x) quantile(x, prob = quantiles)))
      diags <- matrix(c(object$rhat$theta, round(object$ESS$theta, 0)), ncol = 2)
      colnames(diags) <- c('Rhat', 'ESS')
      print(noquote(round(cbind(tmp.1, tmp, diags), digits)))
    }
  }
}

--
You received this message because you are subscribed to the Google Groups "spOccupancy and spAbundance users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to spocc-spabund-u...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/spocc-spabund-users/CAP5oAEpmGHrZZNtsJufuN8g1e-yDJbrr-4%3DM%2BM7C-nLpOPZMCw%40mail.gmail.com.

--

Jeffrey W. Doser, Ph.D.

Assistant Professor

Department of Forestry and Environmental Resources

North Carolina State University

Statistical Ecology and Forest Science Lab

Pronouns: he/him/his

[Marc Kéry]

Dear Jeff and Michael,

thanks for your helpful suggestions.

Best regards  --- Marc

--

______________________________________________________________
 
Marc Kéry
Tel. ++41 41 462 97 93
marc...@vogelwarte.ch
www.vogelwarte.ch
 
Swiss Ornithological Institute | Seerose 1 | CH-6204 Sempach | Switzerland
______________________________________________________________

*** Hierarchical modeling in ecology ***

https://www.mbr-pwrc.usgs.gov/pubanalysis/roylebook/

[Marc Kéry]

Dear Jeff,

that may be a silly question ... but at least if it is, then there's also an easy solution: so, how can I change the function to print into a file instead of outputting on the console ?

Thanks and best regards  ---- Marc 

El lun, 21 jul 2025 a la(s) 6:56 p.m., Jeffrey Doser (jwd...@ncsu.edu) escribió:

[Jeffrey Doser]

Hi Marc,

I think the following should work, where "out" is the name of the model output.

sink("file_name.txt")

summary(out)

sink()

Jeff

[Marc Kéry]

Dear Jeff,

thanks for the R help --- that does work !

Best regards  -- Marc