diff --git a/NAMESPACE b/NAMESPACE
index 376a089..b91e8ad 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,10 +1,8 @@
 # Generated by roxygen2: do not edit by hand
 
 S3method(check,survey)
-S3method(cite,survey)
 S3method(clean,survey)
 export(check)
-export(cite)
 export(clean)
 export(contact_matrix)
 export(download_survey)
@@ -44,7 +42,6 @@ importFrom(httr,status_code)
 importFrom(httr,user_agent)
 importFrom(jsonlite,fromJSON)
 importFrom(jsonlite,toJSON)
-importFrom(lifecycle,deprecate_warn)
 importFrom(lubridate,period)
 importFrom(lubridate,period_to_seconds)
 importFrom(lubridate,years)
diff --git a/R/check.r b/R/check.r
index 98fef67..5fc563f 100644
--- a/R/check.r
+++ b/R/check.r
@@ -7,7 +7,6 @@ check <- function(x, ...) UseMethod("check")
 #' @description Checks that a survey fulfills all the requirements to work with the 'contact_matrix' function
 #'
 #' @param x A [survey()] object
-#' @param columns deprecated argument, ignored
 #' @param id.column the column in both the `participants` and `contacts` data frames that links contacts to participants
 #' @param participant.age.column the column in the `participants` data frame containing participants' age; if this does not exist, at least columns "..._exact", "..._est_min" and "..._est_max" must (see the `estimated.participant.age` option in [contact_matrix()])
 #' @param country.column the column in the `participants` data frame containing the country in which the participant was queried
@@ -19,20 +18,12 @@ check <- function(x, ...) UseMethod("check")
 #' data(polymod)
 #' check(polymod)
 #' @export
-check.survey <- function(x, columns, id.column = "part_id", participant.age.column = "part_age", country.column = "country", year.column = "year", contact.age.column = "cnt_age", ...) {
+check.survey <- function(x, id.column = "part_id", participant.age.column = "part_age", country.column = "country", year.column = "year", contact.age.column = "cnt_age", ...) {
   chkDots(...)
   if (!is.data.frame(x$participants) || !is.data.frame(x$contacts)) {
     stop("The 'participants' and 'contacts' elements of 'x' must be data.frames")
   }
 
-  if (!missing(columns)) {
-    warning(
-      "The 'columns' argument is deprecated and will cause an error from ",
-      "version 1.0.0. The behaviour of the function now always corresponds ",
-      "to the previous documented case for `columns = TRUE`"
-    )
-  }
-
   x <- clean(x)
 
   success <- TRUE
@@ -49,7 +40,7 @@ check.survey <- function(x, columns, id.column = "part_id", participant.age.colu
     exact.column <- paste(participant.age.column, "exact", sep = "_")
     min.column <- paste(participant.age.column, "est_min", sep = "_")
     max.column <- paste(participant.age.column, "est_max", sep = "_")
-    
+
     if (!((exact.column %in% colnames(x$participants)) ||
       (min.column %in% colnames(x$participants) && max.column %in% colnames(x$participants)))) {
       warning(
diff --git a/R/cite.r b/R/cite.r
index 1a8f398..b3d35bb 100644
--- a/R/cite.r
+++ b/R/cite.r
@@ -1,30 +1,3 @@
-#' @export
-cite <- function(x, ...) UseMethod("cite")
-#' @name cite
-#' @rdname cite
-#' @title Citation for a survey
-#'
-#' @description Gets a full citation for a [survey()].
-#'
-#' @param x a character vector of surveys to cite
-#' @param ... ignored
-#' @return citation as bibentry
-#' @importFrom utils bibentry
-#' @importFrom httr GET content
-#' @examples
-#' data(polymod)
-#' cite(polymod)
-#' @export
-cite.survey <- function(x, ...) {
-  warning(
-    "The cite function is deprecated and will stop working in version 1.0.0. ",
-    "Please use get_citation() instead."
-  )
-  chkDots(...)
-
-  get_citation(x)
-}
-
 #' @title Citation for a survey
 #'
 #' @description Gets a full citation for a [survey()].
diff --git a/R/contact_matrix.r b/R/contact_matrix.r
index de28c1f..61311d8 100644
--- a/R/contact_matrix.r
+++ b/R/contact_matrix.r
@@ -1,14 +1,12 @@
 #' Generate a contact matrix from diary survey data
 #'
-#' Samples a contact survey using a bootstrap
+#' Samples a contact survey
 #'
 #' @param survey a [survey()] object
 #' @param countries limit to one or more countries; if not given, will use all countries in the survey; these can be given as country names or 2-letter (ISO Alpha-2) country codes
 #' @param survey.pop survey population -- either a data frame with columns 'lower.age.limit' and 'population', or a character vector giving the name(s) of a country or countries from the list that can be obtained via `wpp_countries`; if not given, will use the country populations from the chosen countries, or all countries in the survey if `countries` is not given
 #' @param age.limits lower limits of the age groups over which to construct the matrix
 #' @param filter any filters to apply to the data, given as list of the form (column=filter_value) - only contacts that have 'filter_value' in 'column' will be considered. If multiple filters are given, they are all applied independently and in the sequence given.
-#' @param n deprecated; number of bootstrap samples to generate
-#' @param bootstrap  deprecated; whether to bootstrap contact matrices
 #' @param counts whether to return counts (instead of means)
 #' @param symmetric whether to make matrix symmetric, such that \eqn{c_{ij}N_i = c_{ji}N_j}.
 #' @param split whether to split the contact matrix into the mean number of contacts, in each age group (split further into the product of the mean number of contacts across the whole population (`mean.contacts`), a normalisation constant (`normalisation`) and age-specific variation in contacts (`contacts`)), multiplied with an assortativity matrix (`assortativity`) and a population multiplier (`demograpy`). For more detail on this, see the "Getting Started" vignette.
@@ -21,7 +19,7 @@
 #' @param weigh.dayofweek whether to weigh social contacts data by the day of the week (weight (5/7 / N_week / N) for weekdays and (2/7 / N_weekend / N) for weekends)
 #' @param weigh.age whether to weigh social contacts data by the age of the participants (vs. the populations' age distribution)
 #' @param weight.threshold threshold value for the standardized weights before running an additional standardisation (default 'NA' = no cutoff)
-#' @param sample.all.age.groups what to do if bootstrapping fails to sample participants from one or more age groups; if FALSE (default), corresponding rows will be set to NA, if TRUE the sample will be discarded and a new one taken instead
+#' @param sample.all.age.groups what to do if sampling participants (with `sample.participants = TRUE`) fails to sample participants from one or more age groups; if FALSE (default), corresponding rows will be set to NA, if TRUE the sample will be discarded and a new one taken instead
 #' @param return.part.weights boolean to return the participant weights
 #' @param return.demography boolean to explicitly return demography data that corresponds to the survey data (default 'NA' = if demography data is requested by other function parameters)
 #' @param per.capita whether to return a matrix with contact rates per capita (default is FALSE and not possible if 'counts=TRUE' or 'split=TRUE')
@@ -30,7 +28,6 @@
 #' @importFrom stats xtabs runif median
 #' @importFrom utils data globalVariables
 #' @importFrom countrycode countrycode
-#' @importFrom lifecycle deprecate_warn
 #' @import data.table
 #' @export
 #' @autoglobal
@@ -38,7 +35,7 @@
 #' data(polymod)
 #' contact_matrix(polymod, countries = "United Kingdom", age.limits = c(0, 1, 5, 15))
 #' @author Sebastian Funk
-contact_matrix <- function(survey, countries = NULL, survey.pop, age.limits, filter, n = 1, bootstrap, counts = FALSE, symmetric = FALSE, split = FALSE, sample.participants = FALSE, estimated.participant.age = c("mean", "sample", "missing"), estimated.contact.age = c("mean", "sample", "missing"), missing.participant.age = c("remove", "keep"), missing.contact.age = c("remove", "sample", "keep", "ignore"), weights = NULL, weigh.dayofweek = FALSE, weigh.age = FALSE, weight.threshold = NA, sample.all.age.groups = FALSE, return.part.weights = FALSE, return.demography = NA, per.capita = FALSE, ...) {
+contact_matrix <- function(survey, countries = NULL, survey.pop, age.limits, filter, counts = FALSE, symmetric = FALSE, split = FALSE, sample.participants = FALSE, estimated.participant.age = c("mean", "sample", "missing"), estimated.contact.age = c("mean", "sample", "missing"), missing.participant.age = c("remove", "keep"), missing.contact.age = c("remove", "sample", "keep", "ignore"), weights = NULL, weigh.dayofweek = FALSE, weigh.age = FALSE, weight.threshold = NA, sample.all.age.groups = FALSE, return.part.weights = FALSE, return.demography = NA, per.capita = FALSE, ...) {
   surveys <- c("participants", "contacts")
 
   dot.args <- list(...)
@@ -57,41 +54,11 @@ contact_matrix <- function(survey, countries = NULL, survey.pop, age.limits, fil
   missing.participant.age <- match.arg(missing.participant.age)
   missing.contact.age <- match.arg(missing.contact.age)
 
-  error_string <-
-    "must be a survey object (created using `survey()` or `get_survey()`)"
-
-  if (is_doi(survey)) {
-    deprecate_warn(
-      "0.3.2", paste0("contact_matrix(survey = '", error_string, "')"),
-      details = "Passing a DOI will be removed in version 0.4.0."
+  if (!inherits(survey, "survey")) {
+    stop(
+      "`survey` must be a survey object (created using `survey()` ",
+      "or `get_survey()`)"
     )
-    survey <- get_survey(survey)
-  } else if (!inherits(survey, "survey")) {
-    stop(error_string)
-  }
-
-  if (!missing(n)) {
-    warning(
-      "The 'n' option is being deprecated and will be removed ",
-      "in version 1.0.0. Please see the ",
-      "'Bootstrapping' section in the vignette for an ",
-      "alternative approach."
-    )
-    if (n > 1) bootstrap <- TRUE
-  }
-  if (!missing(bootstrap)) {
-    warning(
-      "The 'bootstrap' option is being deprecated and will be removed ",
-      "in version 1.0.0. Please use the 'sample.participants'",
-      " option instead."
-    )
-    if (missing(sample.participants)) sample.participants <- bootstrap
-    if (bootstrap != sample.participants) {
-      stop(
-        "'bootstrap' (if given) and 'sample.participants must have the ",
-        "same value."
-      )
-    }
   }
 
   if (!missing(age.limits)) {
@@ -256,7 +223,7 @@ contact_matrix <- function(survey, countries = NULL, survey.pop, age.limits, fil
   if (missing.contact.age == "remove" &&
     nrow(survey$contacts[is.na(get(columns[["contact.age"]])) |
       get(columns[["contact.age"]]) < min(age.limits)]) > 0) {
-    if (n == 1 && !missing.contact.age.set) {
+    if (!missing.contact.age.set) {
       message(
         "Removing participants that have contacts without age information. ",
         "To change this behaviour, set the 'missing.contact.age' option"
@@ -274,7 +241,7 @@ contact_matrix <- function(survey, countries = NULL, survey.pop, age.limits, fil
   if (missing.contact.age == "ignore" &&
     nrow(survey$contacts[is.na(get(columns[["contact.age"]])) |
       get(columns[["contact.age"]]) < min(age.limits)]) > 0) {
-    if (n == 1 && !missing.contact.age.set) {
+    if (!missing.contact.age.set) {
       message(
         "Ignore contacts without age information. ",
         "To change this behaviour, set the 'missing.contact.age' option"
@@ -616,192 +583,179 @@ contact_matrix <- function(survey, countries = NULL, survey.pop, age.limits, fil
       right = FALSE
     )]
 
-  ## Bootstrap
-  if (n > 1 && !bootstrap) {
-    warning(
-      "n > 1 does not make sense if not bootstrapping. ",
-      "Will return just one sample."
-    )
-    n <- 1
-  }
   ret <- list()
-  for (i in seq_len(n)) {
-    if (sample.participants) {
-      good.sample <- FALSE
-      while (!good.sample) {
-        ## take a sample from the participants
-        part.sample <- sample(participant_ids, replace = TRUE)
-        part.age.limits <-
-          unique(survey$participants[
-            get(columns[["id"]]) %in% part.sample,
-            lower.age.limit
-          ])
-        good.sample <- !sample.all.age.groups ||
-          (length(setdiff(age.limits, part.age.limits)) == 0)
-
-        sample.table <-
-          data.table(id = part.sample, weight = 1)
-        sample.table <-
-          sample.table[, list(bootstrap.weight = sum(weight)), by = id]
-        setnames(sample.table, "id", columns[["id"]])
-        setkeyv(sample.table, columns[["id"]])
-
-        sampled.contacts <- merge(survey$contacts, sample.table)
-        sampled.contacts[, sampled.weight := weight * bootstrap.weight]
-
-        sampled.participants <-
-          merge(survey$participants, sample.table)
-        sampled.participants[, sampled.weight := weight * bootstrap.weight]
-      }
-    } else {
-      ## just use all participants
-      sampled.contacts <- survey$contacts
-      sampled.contacts[, sampled.weight := weight]
-      sampled.participants <- survey$participants
-      sampled.participants[, sampled.weight := weight]
+  if (sample.participants) {
+    good.sample <- FALSE
+    while (!good.sample) {
+      ## take a sample from the participants
+      part.sample <- sample(participant_ids, replace = TRUE)
+      part.age.limits <-
+        unique(survey$participants[
+          get(columns[["id"]]) %in% part.sample,
+          lower.age.limit
+        ])
+      good.sample <- !sample.all.age.groups ||
+        (length(setdiff(age.limits, part.age.limits)) == 0)
+
+      sample.table <-
+        data.table(id = part.sample, weight = 1)
+      sample.table <-
+        sample.table[, list(bootstrap.weight = sum(weight)), by = id]
+      setnames(sample.table, "id", columns[["id"]])
+      setkeyv(sample.table, columns[["id"]])
+
+      sampled.contacts <- merge(survey$contacts, sample.table)
+      sampled.contacts[, sampled.weight := weight * bootstrap.weight]
+
+      sampled.participants <-
+        merge(survey$participants, sample.table)
+      sampled.participants[, sampled.weight := weight * bootstrap.weight]
     }
+  } else {
+    ## just use all participants
+    sampled.contacts <- survey$contacts
+    sampled.contacts[, sampled.weight := weight]
+    sampled.participants <- survey$participants
+    sampled.participants[, sampled.weight := weight]
+  }
 
-    ## calculate weighted contact matrix
-    weighted.matrix <-
-      xtabs(
-        data = sampled.contacts,
-        formula = sampled.weight ~ age.group + contact.age.group,
-        addNA = TRUE
-      )
-
-    dims <- dim(weighted.matrix)
-    dim.names <- dimnames(weighted.matrix)
+  ## calculate weighted contact matrix
+  weighted.matrix <-
+    xtabs(
+      data = sampled.contacts,
+      formula = sampled.weight ~ age.group + contact.age.group,
+      addNA = TRUE
+    )
 
-    if (!counts) { ## normalise to give mean number of contacts
-      ## calculate normalisation vector
-      norm.vector <-
-        xtabs(
-          data = sampled.participants,
-          formula = sampled.weight ~ age.group, addNA = TRUE
-        )
+  dims <- dim(weighted.matrix)
+  dim.names <- dimnames(weighted.matrix)
 
-      ## normalise contact matrix
-      weighted.matrix <-
-        array(apply(weighted.matrix, 2, function(x) x / norm.vector),
-          dim = dims,
-          dimnames = dim.names
-        )
-      ## set non-existent data to NA
-      weighted.matrix[is.nan(weighted.matrix)] <- NA_real_
-    }
+  if (!counts) { ## normalise to give mean number of contacts
+    ## calculate normalisation vector
+    norm.vector <-
+      xtabs(
+        data = sampled.participants,
+        formula = sampled.weight ~ age.group, addNA = TRUE
+      )
 
-    ## construct a warning in case there are NAs
-    na.headers <- anyNA(dimnames(weighted.matrix), recursive = TRUE)
-    na.content <- anyNA(weighted.matrix)
-    na.present <- na.headers || na.content
-
-    if (na.present) {
-      warning.suggestion <- "  Consider "
-      if (na.headers) {
-        warning.suggestion <- paste0(warning.suggestion, "setting ")
-        suggested.options <- NULL
-        if (anyNA(rownames(weighted.matrix))) {
-          suggested.options <- c(suggested.options, "'missing.participant.age'")
-        }
-        if (anyNA(colnames(weighted.matrix))) {
-          suggested.options <- c(suggested.options, "'missing.contact.age'")
-        }
+    ## normalise contact matrix
+    weighted.matrix <-
+      array(apply(weighted.matrix, 2, function(x) x / norm.vector),
+        dim = dims,
+        dimnames = dim.names
+      )
+    ## set non-existent data to NA
+    weighted.matrix[is.nan(weighted.matrix)] <- NA_real_
+  }
 
-        warning.suggestion <-
-          paste0(warning.suggestion, paste(suggested.options, collapse = " and "))
-        if (na.content) {
-          warning.suggestion <- paste0(warning.suggestion, ", and ")
-        } else {
-          warning.suggestion <- paste0(warning.suggestion, ".")
-        }
+  ## construct a warning in case there are NAs
+  na.headers <- anyNA(dimnames(weighted.matrix), recursive = TRUE)
+  na.content <- anyNA(weighted.matrix)
+  na.present <- na.headers || na.content
+
+  if (na.present) {
+    warning.suggestion <- "  Consider "
+    if (na.headers) {
+      warning.suggestion <- paste0(warning.suggestion, "setting ")
+      suggested.options <- NULL
+      if (anyNA(rownames(weighted.matrix))) {
+        suggested.options <- c(suggested.options, "'missing.participant.age'")
       }
-      if (na.content) {
-        warning.suggestion <- paste0(warning.suggestion, "adjusting the age limits.")
+      if (anyNA(colnames(weighted.matrix))) {
+        suggested.options <- c(suggested.options, "'missing.contact.age'")
       }
-    }
 
-    if (symmetric && prod(dim(as.matrix(weighted.matrix))) > 1) {
-      if (counts) {
-        warning(
-          "'symmetric=TRUE' does not make sense with 'counts=TRUE'; ",
-          "will not make matrix symmetric."
-        )
-      } else if (na.present) {
-        warning(
-          "'symmetric=TRUE' does not work with missing data; ",
-          "will not make matrix symmetric\n",
-          warning.suggestion
-        )
+      warning.suggestion <-
+        paste0(warning.suggestion, paste(suggested.options, collapse = " and "))
+      if (na.content) {
+        warning.suggestion <- paste0(warning.suggestion, ", and ")
       } else {
-        ## set c_{ij} N_i and c_{ji} N_j (which should both be equal) to
-        ## 0.5 * their sum; then c_{ij} is that sum / N_i
-        normalised.weighted.matrix <- diag(survey.pop$population) %*% weighted.matrix
-        weighted.matrix <- 0.5 * diag(1 / survey.pop$population) %*%
-          (normalised.weighted.matrix + t(normalised.weighted.matrix))
+        warning.suggestion <- paste0(warning.suggestion, ".")
       }
     }
+    if (na.content) {
+      warning.suggestion <- paste0(warning.suggestion, "adjusting the age limits.")
+    }
+  }
 
-    ret[[i]] <- list()
+  if (symmetric && prod(dim(as.matrix(weighted.matrix))) > 1) {
+    if (counts) {
+      warning(
+        "'symmetric=TRUE' does not make sense with 'counts=TRUE'; ",
+        "will not make matrix symmetric."
+      )
+    } else if (na.present) {
+      warning(
+        "'symmetric=TRUE' does not work with missing data; ",
+        "will not make matrix symmetric\n",
+        warning.suggestion
+      )
+    } else {
+      ## set c_{ij} N_i and c_{ji} N_j (which should both be equal) to
+      ## 0.5 * their sum; then c_{ij} is that sum / N_i
+      normalised.weighted.matrix <- diag(survey.pop$population) %*% weighted.matrix
+      weighted.matrix <- 0.5 * diag(1 / survey.pop$population) %*%
+        (normalised.weighted.matrix + t(normalised.weighted.matrix))
+    }
+  }
 
-    if (split) {
-      if (counts) {
-        warning(
-          "'split=TRUE' does not make sense with 'counts=TRUE'; ",
-          "will not split the contact matrix."
-        )
-      } else if (na.present) {
-        warning(
-          "'split=TRUE' does not work with missing data; ",
-          "will not split contact.matrix.\n",
-          warning.suggestion
-        )
-        ret[[i]][["mean.contacts"]] <- NA
-        ret[[i]][["normalisation"]] <- NA
-        ret[[i]][["contacts"]] <- rep(NA, nrow(weighted.matrix))
-      } else {
-        ## get rid of name but preserve row and column names
-        weighted.matrix <- unname(weighted.matrix)
+  if (split) {
+    if (counts) {
+      warning(
+        "'split=TRUE' does not make sense with 'counts=TRUE'; ",
+        "will not split the contact matrix."
+      )
+    } else if (na.present) {
+      warning(
+        "'split=TRUE' does not work with missing data; ",
+        "will not split contact.matrix.\n",
+        warning.suggestion
+      )
+      ret[["mean.contacts"]] <- NA
+      ret[["normalisation"]] <- NA
+      ret[["contacts"]] <- rep(NA, nrow(weighted.matrix))
+    } else {
+      ## get rid of name but preserve row and column names
+      weighted.matrix <- unname(weighted.matrix)
 
-        nb.contacts <- rowSums(weighted.matrix)
-        mean.contacts <- sum(survey.pop$population * nb.contacts) /
+      nb.contacts <- rowSums(weighted.matrix)
+      mean.contacts <- sum(survey.pop$population * nb.contacts) /
           sum(survey.pop$population)
-        spectrum.matrix <- weighted.matrix
-        spectrum.matrix[is.na(spectrum.matrix)] <- 0
-        spectrum <- as.numeric(eigen(spectrum.matrix, only.values = TRUE)$values[1])
-        ret[[i]][["mean.contacts"]] <- mean.contacts
-        ret[[i]][["normalisation"]] <- spectrum / mean.contacts
-
-        age.proportions <- survey.pop$population / sum(survey.pop$population)
-        weighted.matrix <-
-          diag(1 / nb.contacts) %*% weighted.matrix %*% diag(1 / age.proportions)
-        nb.contacts <- nb.contacts / spectrum
-        ret[[i]][["contacts"]] <- nb.contacts
-      }
+      spectrum.matrix <- weighted.matrix
+      spectrum.matrix[is.na(spectrum.matrix)] <- 0
+      spectrum <- as.numeric(eigen(spectrum.matrix, only.values = TRUE)$values[1])
+      ret[["mean.contacts"]] <- mean.contacts
+      ret[["normalisation"]] <- spectrum / mean.contacts
+
+      age.proportions <- survey.pop$population / sum(survey.pop$population)
+      weighted.matrix <-
+        diag(1 / nb.contacts) %*% weighted.matrix %*% diag(1 / age.proportions)
+      nb.contacts <- nb.contacts / spectrum
+      ret[["contacts"]] <- nb.contacts
     }
+  }
+  # make sure the dim.names are retained after symmetric or split procedure
+  dimnames(weighted.matrix) <- dim.names
 
-    # make sure the dim.names are retained after symmetric or split procedure
-    dimnames(weighted.matrix) <- dim.names
-    
-    ret[[i]][["matrix"]] <- weighted.matrix
-
-    # option to add matrix per capita, i.e. the contact rate of age i with one individual of age j in the population.
-    if (per.capita) {
-      if (counts) {
-        warning(
-          "'per.capita=TRUE' does not make sense with 'counts=TRUE'; ",
-          "will not return the contact matrix per capita."
-        )
-      } else if (split) {
-        warning(
-          "'per.capita=TRUE' does not make sense with 'split=TRUE'; ",
-          "will not return the contact matrix per capita."
-        )
-      } else {
-        survey.pop$population
-        weighted.matrix.per.capita <- weighted.matrix / matrix(rep(survey.pop$population, nrow(survey.pop)), ncol = nrow(survey.pop), byrow = TRUE)
-        weighted.matrix.per.capita
-        ret[[i]][["matrix.per.capita"]] <- weighted.matrix.per.capita
-      }
+  ret[["matrix"]] <- weighted.matrix
+
+  # option to add matrix per capita, i.e. the contact rate of age i with one individual of age j in the population.
+  if (per.capita) {
+    if (counts) {
+      warning(
+        "'per.capita=TRUE' does not make sense with 'counts=TRUE'; ",
+        "will not return the contact matrix per capita."
+      )
+    } else if (split) {
+      warning(
+        "'per.capita=TRUE' does not make sense with 'split=TRUE'; ",
+        "will not return the contact matrix per capita."
+      )
+    } else {
+      survey.pop$population
+      weighted.matrix.per.capita <- weighted.matrix / matrix(rep(survey.pop$population, nrow(survey.pop)), ncol = nrow(survey.pop), byrow = TRUE)
+      weighted.matrix.per.capita
+      ret[["matrix.per.capita"]] <- weighted.matrix.per.capita
     }
   }
 
@@ -828,20 +782,13 @@ contact_matrix <- function(survey, countries = NULL, survey.pop, age.limits, fil
   setnames(part.pop, c("age.group", "participants"))
   part.pop[, proportion := participants / sum(participants)]
 
-  # set function output
-  if (length(ret) > 1) {
-    return_value <- list(matrices = ret)
-  } else {
-    return_value <- ret[[1]]
-  }
-
-  if (!is.null(return_value)) {
+  if (!is.null(ret)) {
     if (need.survey.pop && (is.na(return.demography) || return.demography)) {
       # change survey.pop$age.group factors into characters (cfr. part.pop)
       survey.pop[, age.group := as.character(age.group)]
-      return_value[["demography"]] <- survey.pop[]
+      ret[["demography"]] <- survey.pop[]
     }
-    return_value[["participants"]] <- part.pop[]
+    ret[["participants"]] <- part.pop[]
   }
 
   # option to return participant weights
@@ -865,10 +812,10 @@ contact_matrix <- function(survey, countries = NULL, survey.pop, age.limits, fil
     # set name of last column
     names(part.weights)[ncol(part.weights)] <- "participants"
 
-    # add proportion and add to return_value
+    # add proportion and add to ret
     part.weights[, proportion := participants / sum(participants)]
-    return_value[["participants.weights"]] <- part.weights[]
+    ret[["participants.weights"]] <- part.weights[]
   }
 
-  return(return_value)
+  return(ret)
 }
diff --git a/R/load_survey.r b/R/load_survey.r
index 0c0e151..193df16 100644
--- a/R/load_survey.r
+++ b/R/load_survey.r
@@ -24,12 +24,12 @@ load_survey <- function(files, ...) {
   }
   survey_files <- grep("csv$", files, value = TRUE) # select csv files
   reference_file <- grep("json$", files, value = TRUE) # select json file
-  if(length(reference_file)>0){
+  if (length(reference_file) > 0) {
     reference <- fromJSON(reference_file)
   } else {
     reference <- NULL
   }
-  
+
   contact_data <- lapply(survey_files, fread)
   names(contact_data) <- survey_files
 
diff --git a/R/matrix_plot.R b/R/matrix_plot.R
index f9e195b..a742c89 100644
--- a/R/matrix_plot.R
+++ b/R/matrix_plot.R
@@ -29,7 +29,7 @@
 #' }
 #' @author Lander Willem
 matrix_plot <- function(mij, min.legend = 0, max.legend = NA, num.digits = 2, num.colors = 50, main, xlab, ylab, legend.width, legend.mar, legend.shrink, cex.lab, cex.axis, cex.text, color.palette = heat.colors) {
-  
+
   # check function arguments
   xlab <- ifelse(!missing(xlab), xlab, "Age group (year)")
   ylab <- ifelse(!missing(ylab), ylab, "Contact age group (year)")
@@ -40,7 +40,7 @@ matrix_plot <- function(mij, min.legend = 0, max.legend = NA, num.digits = 2, nu
   legend.width <- ifelse(!missing(legend.width), legend.width, 1)
   legend.mar <- ifelse(!missing(legend.mar), legend.mar, 5.1)
   legend.shrink <- ifelse(!missing(legend.shrink), legend.shrink, 0.9)
-  
+
   # set colors
   redc <- rev(color.palette(num.colors))
 
@@ -51,47 +51,51 @@ matrix_plot <- function(mij, min.legend = 0, max.legend = NA, num.digits = 2, nu
 
   # set breaks and midpoints
   breaks <- seq(zlim[1], zlim[2], length = num.colors + 1)
-  midpoints <- matrix(breaks[-length(breaks)] + diff(breaks)/2, nrow = 1, ncol = length(breaks)-1)
-  
+  midpoints <- matrix(
+    breaks[-length(breaks)] + diff(breaks) / 2,
+    nrow = 1, ncol = length(breaks) - 1
+  )
+
   # get plot region for matrix and legend based on current graphical parameters
   # note: based on layout from fields::imagePlot
-  char.size <- par()$cin[1]/par()$din[1] # get text character size
+  char.size <- par()$cin[1] / par()$din[1] # get text character size
   offset <- char.size * par()$mar[4] # space between legend and main plot
 
   # set legends' plot region
   legend_plot_region <- par()$plt
   legend_plot_region[2] <- 1 - (legend.mar * char.size)
   legend_plot_region[1] <- legend_plot_region[2] - (legend.width * char.size)
-  
+
   # account for legend.shrink
-  pr <- (legend_plot_region[4] - legend_plot_region[3]) * ((1 - legend.shrink)/2)
+  pr <- (legend_plot_region[4] - legend_plot_region[3]) *
+    ((1 - legend.shrink) / 2)
   legend_plot_region[4] <- legend_plot_region[4] - pr
   legend_plot_region[3] <- legend_plot_region[3] + pr
-  
-  # set main matrix' plot region 
+
+  # set main matrix' plot region
   main_plot_region    <- par()$plt
   main_plot_region[2] <- min(main_plot_region[2], legend_plot_region[1] - offset)
-  
-  # defensive check for main and legends' plot region 
+
+  # defensive check for main and legends' plot region
   dp <- legend_plot_region[2] - legend_plot_region[1]
   legend_plot_region[1] <- min(main_plot_region[2] + offset, legend_plot_region[1])
   legend_plot_region[2] <- legend_plot_region[1] + dp
-  
+
   # store old graphical parameters, and initiate the ones for the main plot
   old.par <- par(no.readonly = TRUE)
   par(plt = main_plot_region)
-  
+
   # add image plot
-  image(mij, 
+  image(mij,
         xlab = xlab,
         ylab = ylab,
         main = main,
         cex.lab = cex.lab,
         breaks = breaks,
-        col = redc, 
+        col = redc,
         xaxt = "n",
         yaxt = "n")
-  
+
   # add axis labels
   plt_ticks <- seq(0, 1, length = nrow(mij))
   axis(2, at = plt_ticks, labels = c(colnames(mij)), cex.axis = cex.axis, tick = FALSE, las = 1)
@@ -110,17 +114,17 @@ matrix_plot <- function(mij, min.legend = 0, max.legend = NA, num.digits = 2, nu
     e_grid <- expand.grid(plt_ticks, plt_ticks)
     text(e_grid, labels = mij, cex = cex.text)
   }
-  
-  # set graphical parameters for the legend 
+
+  # set graphical parameters for the legend
   par(new = TRUE, pty = "m", plt = legend_plot_region, err = -1)
-  
+
   # include legend bar with axis
-  image(x = 1:2, y=breaks, z=midpoints, 
-        xaxt = "n", yaxt = "n", xlab = "", 
-        ylab = "", col = redc, 
-        breaks=breaks)
+  image(x = 1:2, y = breaks, z = midpoints,
+        xaxt = "n", yaxt = "n", xlab = "",
+        ylab = "", col = redc,
+        breaks = breaks)
   axis(side = 4, mgp = c(3, 1, 0), las = 2)
-  
+
   # restore original graphical parameters
-  par(old.par) 
-}
\ No newline at end of file
+  par(old.par)
+}
diff --git a/man/check.Rd b/man/check.Rd
index cd38faf..15ce71d 100644
--- a/man/check.Rd
+++ b/man/check.Rd
@@ -7,7 +7,6 @@
 \usage{
 \method{check}{survey}(
   x,
-  columns,
   id.column = "part_id",
   participant.age.column = "part_age",
   country.column = "country",
@@ -19,8 +18,6 @@
 \arguments{
 \item{x}{A \code{\link[=survey]{survey()}} object}
 
-\item{columns}{deprecated argument, ignored}
-
 \item{id.column}{the column in both the \code{participants} and \code{contacts} data frames that links contacts to participants}
 
 \item{participant.age.column}{the column in the \code{participants} data frame containing participants' age; if this does not exist, at least columns "..._exact", "..._est_min" and "..._est_max" must (see the \code{estimated.participant.age} option in \code{\link[=contact_matrix]{contact_matrix()}})}
diff --git a/man/cite.Rd b/man/cite.Rd
deleted file mode 100644
index aace7b6..0000000
--- a/man/cite.Rd
+++ /dev/null
@@ -1,24 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/cite.r
-\name{cite}
-\alias{cite}
-\alias{cite.survey}
-\title{Citation for a survey}
-\usage{
-\method{cite}{survey}(x, ...)
-}
-\arguments{
-\item{x}{a character vector of surveys to cite}
-
-\item{...}{ignored}
-}
-\value{
-citation as bibentry
-}
-\description{
-Gets a full citation for a \code{\link[=survey]{survey()}}.
-}
-\examples{
-data(polymod)
-cite(polymod)
-}
diff --git a/man/contact_matrix.Rd b/man/contact_matrix.Rd
index 9139acd..635b371 100644
--- a/man/contact_matrix.Rd
+++ b/man/contact_matrix.Rd
@@ -10,8 +10,6 @@ contact_matrix(
   survey.pop,
   age.limits,
   filter,
-  n = 1,
-  bootstrap,
   counts = FALSE,
   symmetric = FALSE,
   split = FALSE,
@@ -42,10 +40,6 @@ contact_matrix(
 
 \item{filter}{any filters to apply to the data, given as list of the form (column=filter_value) - only contacts that have 'filter_value' in 'column' will be considered. If multiple filters are given, they are all applied independently and in the sequence given.}
 
-\item{n}{deprecated; number of bootstrap samples to generate}
-
-\item{bootstrap}{deprecated; whether to bootstrap contact matrices}
-
 \item{counts}{whether to return counts (instead of means)}
 
 \item{symmetric}{whether to make matrix symmetric, such that \eqn{c_{ij}N_i = c_{ji}N_j}.}
@@ -70,7 +64,7 @@ contact_matrix(
 
 \item{weight.threshold}{threshold value for the standardized weights before running an additional standardisation (default 'NA' = no cutoff)}
 
-\item{sample.all.age.groups}{what to do if bootstrapping fails to sample participants from one or more age groups; if FALSE (default), corresponding rows will be set to NA, if TRUE the sample will be discarded and a new one taken instead}
+\item{sample.all.age.groups}{what to do if sampling participants (with \code{sample.participants = TRUE}) fails to sample participants from one or more age groups; if FALSE (default), corresponding rows will be set to NA, if TRUE the sample will be discarded and a new one taken instead}
 
 \item{return.part.weights}{boolean to return the participant weights}
 
@@ -84,7 +78,7 @@ contact_matrix(
 a contact matrix, and the underlying demography of the surveyed population
 }
 \description{
-Samples a contact survey using a bootstrap
+Samples a contact survey
 }
 \examples{
 data(polymod)
diff --git a/tests/testthat/test-agegroups.r b/tests/testthat/test-agegroups.r
index efb3595..b7d257a 100644
--- a/tests/testthat/test-agegroups.r
+++ b/tests/testthat/test-agegroups.r
@@ -21,7 +21,7 @@ test_that("age groups are ordered factors", {
   ages <- seq_len(50)
   age_limits <- c(0, 5, 10)
   groups <- reduce_agegroups(ages, age_limits)
-  age_groups <- limits_to_agegroups(groups)
+  age_groups <- limits_to_agegroups(groups, notation = "dashes")
   expect_s3_class(age_groups, "ordered")
   expect_s3_class(age_groups, "factor")
 })
diff --git a/tests/testthat/test-checks.r b/tests/testthat/test-checks.r
index 7eedb15..85e110e 100644
--- a/tests/testthat/test-checks.r
+++ b/tests/testthat/test-checks.r
@@ -18,15 +18,11 @@ erroneous_structure1$participants$part_id <- NULL
 erroneous_structure2 <- copy(erroneous_survey)
 erroneous_structure2$participants$part_age <- NULL
 erroneous_structure3 <- copy(erroneous_survey)
-erroneous_structure3$contacts$cnt_age_exact <- NULL 
+erroneous_structure3$contacts$cnt_age_exact <- NULL
 erroneous_structure3$contacts$cnt_age_est_min <- NULL
 
 test_that("incorrect structure of data frames is correctly identified", {
   expect_warning(check(erroneous_structure1), "does not exist")
-  expect_warning(check(erroneous_structure2), "do not exist") 
-  expect_warning(check(erroneous_structure3), "do not exist") 
-})
-
-test_that("deprecated arguments are warned about", {
-  expect_warning(check(polymod, columns = TRUE), "deprecated")
+  expect_warning(check(erroneous_structure2), "do not exist")
+  expect_warning(check(erroneous_structure3), "do not exist")
 })
diff --git a/tests/testthat/test-matrix.r b/tests/testthat/test-matrix.r
index 52f8e12..3383e3c 100644
--- a/tests/testthat/test-matrix.r
+++ b/tests/testthat/test-matrix.r
@@ -540,12 +540,3 @@ test_that("Contact matrices per capita are also generated when bootstrapping", {
     ), 2)
   })
 })
-
-test_that("passing a DOI for survey is deprecated", {
-  skip_if_offline("zenodo.org")
-  skip_on_cran()
-  skip_on_ci()
-  lifecycle::expect_deprecated(
-    contact_matrix(survey = "10.5281/zenodo.1095664") # nolint
-  )
-})
diff --git a/tests/testthat/test-surveys.r b/tests/testthat/test-surveys.r
index 0c3020d..1a93b56 100644
--- a/tests/testthat/test-surveys.r
+++ b/tests/testthat/test-surveys.r
@@ -31,7 +31,3 @@ test_that("missing surveys can't be cited", {
 test_that("multiple DOI's cannot be loaded", {
   expect_error(suppressMessages(suppressWarnings(get_survey(c("10.5281/zenodo.1095664", "10.5281/zenodo.1127693"))))) # nolint
 })
-
-test_that("deprecated functions are warned about", {
-  expect_warning(cite(polymod), "deprecated")
-})
diff --git a/vignettes/socialmixr.Rmd b/vignettes/socialmixr.Rmd
index acf70f9..8f7c1d4 100644
--- a/vignettes/socialmixr.Rmd
+++ b/vignettes/socialmixr.Rmd
@@ -16,6 +16,8 @@ vignette: >
 
 ```{r setup, include = FALSE}
 library("knitr")
+library("socialmixr")
+library("ggplot2")
 knitr::opts_chunk$set(
   collapse = TRUE,
   comment = "#>"
@@ -27,28 +29,6 @@ data.table::setDTthreads(1)
 
 # Usage
 
-The latest stable version of the `socialmixr` package is installed via
-
-```{r, eval=FALSE}
-install.packages("socialmixr")
-```
-
-The latest development version of the `socialmixr` package can be installed via
-
-```r
-devtools::install_github("epiforecasts/socialmixr")
-```
-
-To load the package, use
-
-```{r, eval=FALSE}
-library("socialmixr")
-```
-```{r, echo=FALSE}
-suppressWarnings(library("socialmixr"))
-```
-
-
 At the heart of the `socialmixr` package is the `contact_matrix()` function. This extracts a contact matrix from survey data. You can use the `R` help to find out about usage of the `contact_matrix()` function, including a list of examples:
 
 ```{r eval=FALSE}
@@ -99,10 +79,10 @@ survey_countries(polymod)
 If one wishes to get a contact matrix for one or more specific countries, a `countries` argument can be passed to `contact_matrix()`. If this is not done, the different surveys contained in a dataset are combined as if they were one single sample (i.e., not applying any population-weighting by country or other correction).
 
 
-By default, socialmixr uses the POLYMOD survey. A reference for any given survey can be obtained using `cite()`, e.g.
+By default, socialmixr uses the POLYMOD survey. A reference for any given survey can be obtained using `get_citation()`, e.g.
 
 ```{r}
-cite(polymod)
+get_citation(polymod)
 ```
 
 # Bootstrapping
@@ -336,7 +316,7 @@ survey_data["m_i * w_tilde"] <- survey_data$m_i * survey_data$w_tilde
 kable(survey_data)
 
 # remove the weighted number of contacts
-survey_data$"m_i * w_tilde" <- NULL
+survey_data$`m_i * w_tilde` <- NULL
 
 print(paste("weighted average number of contacts:", round(mean(survey_data$m_i * survey_data$w_tilde), digits = 2)))
 ```
@@ -478,9 +458,7 @@ print(paste("weighted average number of contacts after truncation:", round(mean(
 The contact matrices can be plotted by using the `geom_tile()` function of the `ggplot2` package.
 
 ```{r fig.width=4, fig.height=4}
-library("reshape2")
-library("ggplot2")
-df <- melt(mr, varnames = c("age.group", "age.group.contact"), value.name = "contacts")
+df <- reshape2::melt(mr, varnames = c("age.group", "age.group.contact"), value.name = "contacts")
 ggplot(df, aes(x = age.group, y = age.group.contact, fill = contacts)) +
   theme(legend.position = "bottom") +
   geom_tile()