diff --git a/DESCRIPTION b/DESCRIPTION
index dc2f7d2..6da7dbc 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: ergm.ego
-Version: 1.1-729
-Date: 2025-01-23
+Version: 1.1-738
+Date: 2025-01-25
 Title: Fit, Simulate and Diagnose Exponential-Family Random Graph Models to Egocentrically Sampled Network Data
 Authors@R: c(
             person(c("Pavel", "N."), "Krivitsky", role=c("aut","cre"), email="pavel@statnet.org", comment=c(ORCID="0000-0002-9101-3362")),
diff --git a/R/degreedist.R b/R/degreedist.R
index 5434a26..2506f89 100644
--- a/R/degreedist.R
+++ b/R/degreedist.R
@@ -158,6 +158,10 @@ degreedist.egor <- function(object, freq = FALSE, prob = !freq,
 #' @param object A [`egor`] object.
 #' @param attrname A character vector containing the name of the network
 #' attribute whose mixing matrix is wanted.
+#' @param useNA One of `"ifany"`, `"no"` or `"always"`, interpreted as
+#'   in [table()]. By default (\code{useNA = "ifany"}) if there are
+#'   any \code{NA}s on the attribute corresponding row \emph{and}
+#'   column will be contained in the result. See Details.
 #' @param rowprob Whether the counts should be normalized by row sums. That is,
 #' whether they should be proportions conditional on the ego's group.
 #' @param weight Whether sampling weights should be incorporated into
@@ -180,18 +184,23 @@ degreedist.egor <- function(object, freq = FALSE, prob = !freq,
 #' (mm.ego <- mixingmatrix(fmh.ego,"Grade"))
 #' 
 #' @export
-mixingmatrix.egor <- function(object, attrname, rowprob = FALSE, weight = TRUE, ...){
+mixingmatrix.egor <- function(object, attrname, useNA = c("ifany", "no", "always"), rowprob = FALSE, weight = TRUE, ...){
+  useNA <- match.arg(useNA)
   ds <- .degreeseq(object)
+  if(! attrname %in% colnames(object$ego) || ! attrname %in% colnames(object$alter))
+    stop("vertex attribute ", sQuote(attrname), " not found in egocentric dataset ", sQuote(deparse1(substitute(object))))
   egos <- rep(.unfactor(as_tibble(object$ego)[[attrname]]), ds)
   alters <- .unfactor(object$alter[[attrname]])
-  levs <- sort(unique(c(egos,alters)))
+  levs <- sort(unique(c(egos, alters, if(useNA == "always") NA)), na.last = if(useNA == "no") NA else TRUE)
+
+  egos <- match(egos, levs, nomatch = 0)
+  alters <- match(alters, levs, nomatch = 0)
 
   w <- if(weight) rep(weights(object),ds) else rep(1, nrow(object$alter))
 
-  mxmat <- outer(levs, levs, Vectorize(function(l1, l2) sum(w[egos==l1&alters==l2])))
+  mxmat <- outer(seq_along(levs), seq_along(levs), Vectorize(function(l1, l2) sum(w[egos==l1&alters==l2])))
 
-  dimnames(mxmat) <- list(ego = levs,  
-                          alter = levs)
+  dimnames(mxmat) <- list(Ego = levs, Alter = levs)
   if(rowprob){
     mxmat <- mxmat/rowSums(mxmat)
   }
diff --git a/man/mixingmatrix.egor.Rd b/man/mixingmatrix.egor.Rd
index 725da98..ed31357 100644
--- a/man/mixingmatrix.egor.Rd
+++ b/man/mixingmatrix.egor.Rd
@@ -5,7 +5,14 @@
 \alias{mixingmatrix}
 \title{Summarizing the mixing among groups in an egocentric dataset}
 \usage{
-\method{mixingmatrix}{egor}(object, attrname, rowprob = FALSE, weight = TRUE, ...)
+\method{mixingmatrix}{egor}(
+  object,
+  attrname,
+  useNA = c("ifany", "no", "always"),
+  rowprob = FALSE,
+  weight = TRUE,
+  ...
+)
 }
 \arguments{
 \item{object}{A \code{\link[egor:egor]{egor}} object.}
@@ -13,6 +20,11 @@
 \item{attrname}{A character vector containing the name of the network
 attribute whose mixing matrix is wanted.}
 
+\item{useNA}{One of \code{"ifany"}, \code{"no"} or \code{"always"}, interpreted as
+in \code{\link[=table]{table()}}. By default (\code{useNA = "ifany"}) if there are
+any \code{NA}s on the attribute corresponding row \emph{and}
+column will be contained in the result. See Details.}
+
 \item{rowprob}{Whether the counts should be normalized by row sums. That is,
 whether they should be proportions conditional on the ego's group.}
 
diff --git a/tests/testthat/test-mixingmatrix.R b/tests/testthat/test-mixingmatrix.R
index 07b421f..d30b85d 100644
--- a/tests/testthat/test-mixingmatrix.R
+++ b/tests/testthat/test-mixingmatrix.R
@@ -20,6 +20,10 @@ for( v in varnames ) {
     })
 }
 
+test_that("Vertex attribute not found", {
+  expect_error(mixingmatrix(egor32, "abc"),
+               "vertex attribute 'abc' not found in egocentric dataset 'egor32'")
+})
 
 for( v in varnames ) {
   test_that(
@@ -94,3 +98,21 @@ test_that("mixing matrices for FMH and egoFMH are equivalent", {
     }
   )
 })
+
+test_that("mixing matrices for FMH and egoFMH are equivalent with missing data", {
+  data("faux.mesa.high")
+  set.vertex.attribute(faux.mesa.high, "Grade", NA, 1:40)
+  fmh.ego <- as.egor(faux.mesa.high)
+  expect_equal(
+    {
+      mm.ego <- mixingmatrix(fmh.ego, "Grade")
+      names(dimnames(mm.ego)) <- c("From", "To")
+      mm.ego
+    },
+    {
+      mm <- mixingmatrix(faux.mesa.high, "Grade")
+      diag(mm) <- diag(mm) * 2
+      mm
+    }
+  )
+})