diff --git a/R/fitINLA.R b/R/fitINLA.R
index eb9db6b..eb93bf9 100644
--- a/R/fitINLA.R
+++ b/R/fitINLA.R
@@ -170,11 +170,18 @@ smoothDirect <- function(data, Amat, formula = NULL, time.model = c("rw1", "rw2"
             "\n  No temporal components")
 
   }else{
+    T <- year_range[2] - year_range[1] + 1
+    if(!is.yearly){
+      T <- length(year_label)
+    }
     message("----------------------------------",
           "\nSmoothed Direct Model",
-            "\n  Main temporal model:        ", time.model, appendLF = FALSE)
-    msg <- paste0(msg, "\nSmoothed Direct Model",
-            "\n  Main temporal model:        ")
+            "\n  Main temporal model:        ", time.model, 
+            "\n  Number of time periods:     ", T, 
+            appendLF = FALSE)
+    msg <- paste0(msg, "\nSmoothed Direct Model", 
+            "\n  Main temporal model:        ", time.model, 
+            "\n  Number of time periods:     ", T)
     if(m == 1){
       if(is.yearly){
         message("\n  Temporal resolution:        period model (m = 1)", appendLF=FALSE)
@@ -221,8 +228,11 @@ smoothDirect <- function(data, Amat, formula = NULL, time.model = c("rw1", "rw2"
   }
 
   if(is.spatial){ 
-    message("\n  Spatial effect:             bym2", appendLF=FALSE) 
-    msg <- paste0(msg, "\n  Spatial effect:             bym2")
+    message("\n  Spatial effect:             bym2",
+            "\n  Number of regions:          ", dim(Amat)[1], 
+             appendLF=FALSE) 
+    msg <- paste0(msg, "\n  Spatial effect:             bym2",
+            "\n  Number of regions:          ", dim(Amat)[1])
   }
   if(is.spatial && is.temporal){
     message("\n  Interaction temporal model: ", st.time.model, 
diff --git a/R/fitINLA2.R b/R/fitINLA2.R
index 736ce74..bba66e4 100644
--- a/R/fitINLA2.R
+++ b/R/fitINLA2.R
@@ -15,7 +15,7 @@
 #' @param X Covariate matrix. It must contain either a column with name "region", or a column with name "years", or both. The covariates must not have missing values for all regions (if varying in space) and all time periods (if varying in time). The rest of the columns are treated as covariates in the mean model.
 #' @param age.groups a character vector of age groups in increasing order.
 #' @param age.n number of months in each age groups in the same order.
-#' @param age.rw.group vector indicating grouping of the ages groups in the temporal model. For example, if each age group is assigned a different temporal component, then set age.rw.group to c(1:length(age.groups)); if all age groups share the same random walk component, then set age.rw.group to a rep(1, length(age.groups)). The default for 6 age groups is c(1,2,3,3,3,3), which assigns a separate temporal trend to the first two groups and a common random walk for the rest of the age groups. The vector should contain values starting from 1.
+#' @param age.time.group vector indicating grouping of the ages groups in the temporal model. For example, if each age group is assigned a different temporal component, then set age.rw.group to c(1:length(age.groups)); if all age groups share the same random walk component, then set age.rw.group to a rep(1, length(age.groups)). The default for 6 age groups is c(1,2,3,3,3,3), which assigns a separate temporal trend to the first two groups and a common random walk for the rest of the age groups. The vector should contain values starting from 1. This argument replaces the previous \code{age.rw.group} argument. 
 #' @param age.strata.fixed.group vector indicating grouping of the ages groups for different strata in the intercept. The default is c(1:length(age.groups)), which correspond to each age group within each stratum receives a separate intercept. If several age groups are specified to be the same value in this vector, the stratum specific deviation from the baseline is assumed to be the same for these age groups. For example, if \code{age.strata.fixed.group = c(1, 2, 3, 3, 3, 3)}, then the intercept part of the linear predictor consists of 6 overall age-specific intercepts and 3 set of strata effects (where a base stratum is chosen internally), for age groups 1, 2, and the rest respectively. Notice that this argument does not control the linear temporal trends (which is also parameterized as fixed effect, but determined by the \code{age.rw.group} argument). The vector should contain values starting from 1.
 #' 
 #' More specific examples: (1) if each age group is assigned a different intercept, then set age.strata.fixed.group to c(1:length(age.groups)) (2) if all age groups share the same intercept, then set age.strata.fixed.group to a rep(1, length(age.groups)). The default for 6 age groups is the former. (3) If each temporal trend is associated with its own intercept, set it to be the same as \code{age.rw.group}. 
@@ -53,6 +53,7 @@
 #' @param ar Deprecated. Order of the autoregressive component. If ar is specified to be positive integer, the random walk components will be replaced by AR(p) terms in the interaction part. The main temporal trend remains to be random walk of order rw unless rw = 0.
 #' @param st.rw Deprecated. Take values 1 or 2, indicating the order of random walk for the interaction term. If not specified, it will take the same order as the argument rw in the main effect. Notice that this argument is only used if ar is set to 0.
 #' @param geo Deprecated. Spatial polygon file, legacy parameter from previous versions of the package.
+#' @param age.rw.group Deprecated. Legacy parameter replaced by \code{age.time.group}.
 #' @param ... arguments to be passed to the inla() function call.
 #' @seealso \code{\link{getDirect}}
 #' @import Matrix
@@ -140,7 +141,7 @@
 #' 
 #' 
 
-smoothCluster <- function(data, X = NULL, family = c("betabinomial", "binomial")[1], age.groups = c("0", "1-11", "12-23", "24-35", "36-47", "48-59"), age.n = c(1,11,12,12,12,12), age.rw.group = c(1,2,3,3,3,3), age.strata.fixed.group = c(1,2,3,4,5,6), time.model = c("rw1", "rw2", "ar1")[2], st.time.model = NULL, Amat, bias.adj = NULL, bias.adj.by = NULL, formula = NULL, year_label, type.st = 4, survey.effect = FALSE, linear.trend = TRUE, common.trend = FALSE, strata.time.effect = FALSE, hyper = "pc", pc.u = 1, pc.alpha = 0.01, pc.u.phi = 0.5, pc.alpha.phi = 2/3, pc.u.cor = 0.7, pc.alpha.cor = 0.9,  pc.st.u = NA, pc.st.alpha = NA, pc.st.slope.u = NA, pc.st.slope.alpha = NA, overdisp.mean = 0, overdisp.prec = 0.4, options = list(config = TRUE), control.inla = list(strategy = "adaptive", int.strategy = "auto"), control.fixed = list(), verbose = FALSE, geo = NULL, rw = NULL, ar = NULL, st.rw = NULL, ...){
+smoothCluster <- function(data, X = NULL, family = c("betabinomial", "binomial")[1], age.groups = c("0", "1-11", "12-23", "24-35", "36-47", "48-59"), age.n = c(1,11,12,12,12,12), age.time.group = c(1,2,3,3,3,3), age.strata.fixed.group = c(1,2,3,4,5,6), time.model = c("rw1", "rw2", "ar1")[2], st.time.model = NULL, Amat, bias.adj = NULL, bias.adj.by = NULL, formula = NULL, year_label, type.st = 4, survey.effect = FALSE, linear.trend = TRUE, common.trend = FALSE, strata.time.effect = FALSE, hyper = "pc", pc.u = 1, pc.alpha = 0.01, pc.u.phi = 0.5, pc.alpha.phi = 2/3, pc.u.cor = 0.7, pc.alpha.cor = 0.9,  pc.st.u = NA, pc.st.alpha = NA, pc.st.slope.u = NA, pc.st.slope.alpha = NA, overdisp.mean = 0, overdisp.prec = 0.4, options = list(config = TRUE), control.inla = list(strategy = "adaptive", int.strategy = "auto"), control.fixed = list(), verbose = FALSE, geo = NULL, rw = NULL, ar = NULL, st.rw = NULL, age.rw.group = NULL, ...){
 
   # if(family == "betabinomialna") stop("family = betabinomialna is still experimental.")
   # check region names in Amat is consistent
@@ -158,6 +159,20 @@ smoothCluster <- function(data, X = NULL, family = c("betabinomial", "binomial")
     message("Argument geo is deprecated in the smoothCluster function. Only Amat is needed.")
   }
 
+  if(sum(!unique(data$age) %in% age.groups) > 0){
+    stop("The data consist of more age groups than specified by the 'age.groups' argument. Please check your data and subset your data properly.")
+  }
+  if(!is.null(age.rw.group)){
+    age.time.group <- age.rw.group
+    message("Argument 'age.rw.group' have been deprecated and replaced by 'age.time.group' in version 1.4.0. The value for 'age.time.group' has been set to be the input argument 'age.rw.group'")
+  }else{
+    age.rw.group <- age.time.group
+  }
+  check.age <- c(length(age.groups), length(age.n), length(age.time.group), length(age.strata.fixed.group))
+  if(length(unique(check.age)) > 1){
+    stop("The arguments 'age.group', 'age.n', 'age.time.group', and 'age.strata.fixed.group' are of different lengths. Please check the specification of these arguments.")
+  }
+
   # get around CRAN check of using un-exported INLA functions
   rate0 <- shape0 <- my.cache <- inla.as.sparse <- type <- strata <- rescale_U <- sim_alpha <- pc.st.slope.prec.u <- NULL
 
@@ -181,6 +196,7 @@ smoothCluster <- function(data, X = NULL, family = c("betabinomial", "binomial")
   # define whether there is the temporal component
   is.temporal <- TRUE
 
+
   if(!is.null(rw) || !is.null(ar) || !is.null(st.rw)){
 
     if(is.null(ar)) ar <- 0
@@ -272,12 +288,21 @@ smoothCluster <- function(data, X = NULL, family = c("betabinomial", "binomial")
   }else{
     message("----------------------------------",
             "\nCluster-level model",
-            "\n  Main temporal model:        ", time.model, appendLF = FALSE)
+            "\n  Main temporal model:        ", time.model, 
+            "\n  Number of time periods:     ", length(year_label), 
+            appendLF = FALSE)
     msg <- paste0(msg, "\nCluster-level model",
-                       "\n  Main temporal model:        ", time.model)
+                       "\n  Main temporal model:        ", time.model, 
+                       "\n  Number of time periods:     ", length(year_label))
     if(linear.trend && !common.trend){
-        message("\n  Additional linear trends:   stratum-age-specific", appendLF = FALSE)
-        msg <- paste0(msg, "\n  Additional linear trends:   stratum-age-specific")
+        if("strata" %in% colnames(data) == FALSE || sum(!is.na(data$strata)) == 0){
+          message("\n  Additional linear trends:   age-specific", appendLF = FALSE)
+          msg <- paste0(msg, "\n  Additional linear trends:   age-specific")
+        }else{
+          message("\n  Additional linear trends:   stratum-age-specific", appendLF = FALSE)
+          msg <- paste0(msg, "\n  Additional linear trends:   stratum-age-specific")
+        }
+
     }
     if(linear.trend && common.trend){
         message("\n  Additional linear trends:   shared", appendLF = FALSE)
@@ -311,8 +336,11 @@ smoothCluster <- function(data, X = NULL, family = c("betabinomial", "binomial")
   }
 
   if(is.spatial){ 
-    message("\n  Spatial effect:             bym2", appendLF=FALSE) 
-    msg <- paste0(msg, "\n  Spatial effect:             bym2")
+    message("\n  Spatial effect:             bym2", 
+            "\n  Number of regions:          ", dim(Amat)[1], 
+            appendLF=FALSE) 
+    msg <- paste0(msg, "\n  Spatial effect:             bym2", 
+                   "\n  Number of regions:          ", dim(Amat)[1])
   }
   if(is.spatial && is.temporal){
     message("\n  Interaction temporal model: ", st.time.model, 
diff --git a/man/Benchmark.Rd b/man/Benchmark.Rd
index 15af94e..3d1380d 100644
--- a/man/Benchmark.Rd
+++ b/man/Benchmark.Rd
@@ -102,7 +102,7 @@ legend("topleft", c("Before benchmarking", "After benchmarking"), pch = c(1, 10)
 #  construct a simple national estimates
 national <- data.frame(years = periods, 
 					   est = seq(0.22, 0.1, length = 7), 
-					   sd = seq(.01, .1, len = 7))
+					   sd = runif(7, 0.01, 0.03))
 # national does not need to have all years
 national_sub <- national[1:3,]
 
@@ -111,7 +111,7 @@ est.bb.bench <- Benchmark(est.bb, national_sub,
 						weight.region = weight.region, 
 						estVar = "est", sdVar = "sd", timeVar = "years")
 
-# Sanity check: only benchmarked for two years
+# Sanity check: only benchmarked for three periods
 compare <- national
 compare$before <- NA
 compare$after <- NA
@@ -175,7 +175,7 @@ est.bb.new <- getSmoothed(fit.bb.new, nsim = 10000, CI = 0.95, save.draws=TRUE)
 #  construct a simple national estimates
 national <- data.frame(years = periods, 
 					   est = seq(0.22, 0.1, length = 7), 
-					   sd = seq(.01, .1, len = 7))
+					   sd = runif(7, 0.01, 0.03))
 weight.region <- expand.grid(region = unique(counts.all$region), 
 						 years = periods)
 weight.region$proportion <- 0.25					   
@@ -249,7 +249,7 @@ plot(national$est, national.est$mean,
 		 ylim = range(c(national$est, national.est$mean, national.est.bench$mean)),
 		 xlab = "External national estimates", 
 		 ylab = "Weighted posterior median after benchmarking",
-    main = "Sanity check: weighted average of area medians")
+    main = "Sanity check: weighted average of area means")
 points(national$est, national.est.bench$mean, col = 2, pch = 10)
 abline(c(0, 1))
 legend("topleft", c("Before benchmarking", "After benchmarking"), pch = c(1, 10), col = c(1, 2))
@@ -283,7 +283,7 @@ plot(national$est, national.est$mean,
 		 ylim = range(c(national$est, national.est$mean, national.est.bench$mean)),
 		 xlab = "External national estimates", 
 		 ylab = "Weighted posterior median after benchmarking",
-    main = "Sanity check: weighted average of area medians")
+    main = "Sanity check: weighted average of area means")
 points(national$est, national.est.bench$mean, col = 2, pch = 10)
 abline(c(0, 1))
 legend("topleft", c("Before benchmarking", "After benchmarking"), pch = c(1, 10), col = c(1, 2))
diff --git a/man/smoothCluster.Rd b/man/smoothCluster.Rd
index 65ad062..f1a4eb9 100644
--- a/man/smoothCluster.Rd
+++ b/man/smoothCluster.Rd
@@ -11,7 +11,7 @@ smoothCluster(
   family = c("betabinomial", "binomial")[1],
   age.groups = c("0", "1-11", "12-23", "24-35", "36-47", "48-59"),
   age.n = c(1, 11, 12, 12, 12, 12),
-  age.rw.group = c(1, 2, 3, 3, 3, 3),
+  age.time.group = c(1, 2, 3, 3, 3, 3),
   age.strata.fixed.group = c(1, 2, 3, 4, 5, 6),
   time.model = c("rw1", "rw2", "ar1")[2],
   st.time.model = NULL,
@@ -46,6 +46,7 @@ smoothCluster(
   rw = NULL,
   ar = NULL,
   st.rw = NULL,
+  age.rw.group = NULL,
   ...
 )
 
@@ -55,7 +56,7 @@ fitINLA2(
   family = c("betabinomial", "binomial")[1],
   age.groups = c("0", "1-11", "12-23", "24-35", "36-47", "48-59"),
   age.n = c(1, 11, 12, 12, 12, 12),
-  age.rw.group = c(1, 2, 3, 3, 3, 3),
+  age.time.group = c(1, 2, 3, 3, 3, 3),
   age.strata.fixed.group = c(1, 2, 3, 4, 5, 6),
   time.model = c("rw1", "rw2", "ar1")[2],
   st.time.model = NULL,
@@ -90,6 +91,7 @@ fitINLA2(
   rw = NULL,
   ar = NULL,
   st.rw = NULL,
+  age.rw.group = NULL,
   ...
 )
 }
@@ -113,7 +115,7 @@ fitINLA2(
 
 \item{age.n}{number of months in each age groups in the same order.}
 
-\item{age.rw.group}{vector indicating grouping of the ages groups in the temporal model. For example, if each age group is assigned a different temporal component, then set age.rw.group to c(1:length(age.groups)); if all age groups share the same random walk component, then set age.rw.group to a rep(1, length(age.groups)). The default for 6 age groups is c(1,2,3,3,3,3), which assigns a separate temporal trend to the first two groups and a common random walk for the rest of the age groups. The vector should contain values starting from 1.}
+\item{age.time.group}{vector indicating grouping of the ages groups in the temporal model. For example, if each age group is assigned a different temporal component, then set age.rw.group to c(1:length(age.groups)); if all age groups share the same random walk component, then set age.rw.group to a rep(1, length(age.groups)). The default for 6 age groups is c(1,2,3,3,3,3), which assigns a separate temporal trend to the first two groups and a common random walk for the rest of the age groups. The vector should contain values starting from 1. This argument replaces the previous \code{age.rw.group} argument.}
 
 \item{age.strata.fixed.group}{vector indicating grouping of the ages groups for different strata in the intercept. The default is c(1:length(age.groups)), which correspond to each age group within each stratum receives a separate intercept. If several age groups are specified to be the same value in this vector, the stratum specific deviation from the baseline is assumed to be the same for these age groups. For example, if \code{age.strata.fixed.group = c(1, 2, 3, 3, 3, 3)}, then the intercept part of the linear predictor consists of 6 overall age-specific intercepts and 3 set of strata effects (where a base stratum is chosen internally), for age groups 1, 2, and the rest respectively. Notice that this argument does not control the linear temporal trends (which is also parameterized as fixed effect, but determined by the \code{age.rw.group} argument). The vector should contain values starting from 1.
 
@@ -185,6 +187,8 @@ More specific examples: (1) if each age group is assigned a different intercept,
 
 \item{st.rw}{Deprecated. Take values 1 or 2, indicating the order of random walk for the interaction term. If not specified, it will take the same order as the argument rw in the main effect. Notice that this argument is only used if ar is set to 0.}
 
+\item{age.rw.group}{Deprecated. Legacy parameter replaced by \code{age.time.group}.}
+
 \item{...}{arguments to be passed to the inla() function call.}
 }
 \value{