diff --git a/DESCRIPTION b/DESCRIPTION
index fefa6c9..cd82eea 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: bayesPop
 Type: Package
 Title: Probabilistic Population Projection
-Version: 10.0-1.9021
+Version: 10.0-1.9022
 Date: 2024-12-10
 Author: Hana Sevcikova, Adrian Raftery, Thomas Buettner
 Maintainer: Hana Sevcikova <hanas@uw.edu>
diff --git a/R/get_outputs.R b/R/get_outputs.R
index 3e9bb2f..4ee3c67 100644
--- a/R/get_outputs.R
+++ b/R/get_outputs.R
@@ -1784,19 +1784,22 @@ cohorts <- function(pop.pred, country=NULL, expression=NULL, pi=c(80, 95)) {
 	age.index <- as.integer(dimnames(alldata)[[1]])
 	nage <- dim(alldata)[1]
 	years <- dimnames(alldata)[[2]]
-	last.observed.cohort <- pop.pred$proj.years.pop[1]-age.index[1]*5
-	from.cohorts <- seq(last.observed.cohort, length=nage-1, by=-5)
-	observed.cohorts <- paste(from.cohorts, '-', from.cohorts+5, sep="")
+	step <- if(pop.pred$annual) 1 else 5
+	last.observed.cohort <- pop.pred$proj.years.pop[1]-age.index[1]*step
+	from.cohorts <- seq(last.observed.cohort, length=nage-1, by=-step)
+	observed.cohorts <- if(step == 5) paste(from.cohorts, '-', from.cohorts+step, sep="") else as.character(from.cohorts)
 	for(cohort in length(observed.cohorts):1) {
 		cohort.traj <- apply(alldata[cohort:nage,,,drop=FALSE], 3, 'diag')
+		if(is.null(dim(cohort.traj))) next
 		result[[observed.cohorts[cohort]]] <- .get.quantiles.from.cohort.data(cohort.traj)
 		colnames(result[[observed.cohorts[cohort]]]) <- years[1:ncol(result[[observed.cohorts[cohort]]])]
 	}
 	nyears <- length(pop.pred$proj.years.pop)
-	from.cohorts <- seq(last.observed.cohort + 5, length=nyears-2, by=5)
-	projected.cohorts <- paste(from.cohorts, '-', from.cohorts+5, sep="")	
-	for(cohort in 1:length(projected.cohorts)) {
+	from.cohorts <- seq(last.observed.cohort + step, length=if(step == 5) nyears-2 else nyears, by=step)
+	projected.cohorts <- if(step == 5) paste(from.cohorts, '-', from.cohorts+step, sep="")	else as.character(from.cohorts)
+	for(cohort in 1:(length(projected.cohorts) - 1)) {
 		cohort.traj <- apply(alldata[,(cohort+1):nyears,,drop=FALSE], 3, 'diag')
+		if(is.null(dim(cohort.traj))) next
 		result[[projected.cohorts[cohort]]] <- .get.quantiles.from.cohort.data(cohort.traj)
 		colnames(result[[projected.cohorts[cohort]]]) <- years[(cohort+1):(cohort+ncol(result[[projected.cohorts[cohort]]]))]
 	}
diff --git a/R/plot_functions.R b/R/plot_functions.R
index fa35efd..0914e2c 100644
--- a/R/plot_functions.R
+++ b/R/plot_functions.R
@@ -1233,13 +1233,16 @@ pop.cohorts.plot <- function(pop.pred, country=NULL, expression=NULL, cohorts=NU
 		cohort.data <- cohorts(pop.pred, country=country, expression=expression, pi=pi)
 	all.cohorts <- names(cohort.data)[-which(names(cohort.data) == 'last.observed')]
 	all.cohorts.num.start <- as.integer(substr(all.cohorts, 1, 4))
+	step <- if(pop.pred$annual) 1 else 5
 	if(is.null(cohorts)) 
-		cohorts <- seq(cohort.data[['last.observed']], by=5, 
+		cohorts <- seq(cohort.data[['last.observed']], by=step, 
 							length=min(10, sum(all.cohorts.num.start > cohort.data[['last.observed']])))
-	if(any(is.numeric(cohorts))) {
-		# convert to the from-to format, e.g. 2000-2005
-		from.cohorts <- .round.to.lower5(cohorts)
-		cohorts <- paste(from.cohorts, '-', from.cohorts+5, sep="")
+	if(any(is.numeric(cohorts))){
+	    if(!pop.pred$annual) {
+		    # convert to the from-to format, e.g. 2000-2005
+		    from.cohorts <- .round.to.lower5(cohorts)
+		    cohorts <- paste(from.cohorts, '-', from.cohorts+5, sep="")
+	    } else cohorts <- as.character(cohorts)
 	}
 	if(is.null(xlim))
 		xlim <- range(unlist(sapply(cohorts, function(x) range(as.integer(colnames(cohort.data[[x]]))))))