diff --git a/R/tongfen_estimate.R b/R/tongfen_estimate.R
index fa5842b..eb6f749 100644
--- a/R/tongfen_estimate.R
+++ b/R/tongfen_estimate.R
@@ -56,38 +56,26 @@ tongfen_estimate <- function(target,source,meta,na.rm=FALSE) {
 
   gc = which(st_is(i, "GEOMETRYCOLLECTION"))
   i[gc] = st_collection_extract(i[gc], "POLYGON")
-  two_d = which(st_dimension(i) == 2)
-  i[two_d] = st_cast(i[two_d], "MULTIPOLYGON")
 
   source <- source %>% rename(!!!safe_rename_vars)
+  source_area <- unclass(st_area(source))
 
   x_st <- source[idx[,1],, drop=FALSE] %>%
     select(names(safe_rename_vars)) %>%
     pre_scale(meta,meta_var = "var_name") %>%
     mutate(...area_st = st_area(i) %>% unclass,
-           ...area_s = unclass(st_area(.))) %>%
+           ...area_s = source_area[idx[,1]]) %>%
     mutate(...factor = .data$...area_st/.data$...area_s) %>%
     mutate(...partial = .data$...factor < 0.99) %>%
     st_drop_geometry()
 
-
-  for (var in meta$var_name) {
-    # for (ci in naive_CI) {
-    #   c <- ci/100
-    #   x_st <- x_st %>%
-    #     mutate(!!paste0(var,"_lower_",ci) := !!as.name(var) * ifelse(.data$...partial,(1-c) * .data$...factor, 1),
-    #            !!paste0(var,"_upper_",ci) := !!as.name(var) * ifelse(.data$...partial,(1-c) * .data$...factor + c, 1))
-    #
-    # }
-    x_st <- x_st %>%
-      mutate(!!var:=!!as.name(var) * .data$...factor)
-  }
+  x_st[meta$var_name] <- lapply(x_st[meta$var_name], `*`, x_st$...factor)
 
   x_st <- stats::aggregate(x_st, list(idx[,2]), sum, na.rm=na.rm)
 
   result <- target %>%
     left_join(x_st %>%
-                select(-.data$...factor,-.data$...partial,-.data$...area_s,-.data$...area_st) %>%
+                select(-all_of(c("...factor", "...partial", "...area_s", "...area_st"))) %>%
                 rename(!!unique_key:="Group.1"),
               by=unique_key) %>%
     select(-all_of(unique_key)) %>%
diff --git a/tests/testthat/test-estimate.R b/tests/testthat/test-estimate.R
new file mode 100644
index 0000000..c50efe7
--- /dev/null
+++ b/tests/testthat/test-estimate.R
@@ -0,0 +1,47 @@
+square <- function(xmin, xmax, ymin = 0, ymax = 1) {
+  sf::st_polygon(list(rbind(
+    c(xmin, ymin), c(xmax, ymin), c(xmax, ymax),
+    c(xmin, ymax), c(xmin, ymin)
+  )))
+}
+
+test_that("tongfen_estimate area-weights several variables", {
+  source <- sf::st_sf(
+    total = c(10, 20),
+    second = c(100, 200),
+    geometry = sf::st_sfc(square(0, 1), square(1, 2), crs = 3347)
+  )
+  target <- sf::st_sf(
+    region = c("left", "middle"),
+    geometry = sf::st_sfc(square(0, 0.5), square(0.5, 1.5), crs = 3347)
+  )
+  meta <- meta_for_additive_variables(
+    "synthetic",
+    c(total = "total", second = "second")
+  )
+
+  result <- tongfen_estimate(target, source, meta)
+
+  expect_s3_class(result, "sf")
+  expect_equal(result$region, c("left", "middle"))
+  expect_equal(result$total, c(5, 15))
+  expect_equal(result$second, c(50, 150))
+  expect_equal(sf::st_geometry(result), sf::st_geometry(target))
+})
+
+test_that("tongfen_estimate preserves missing-value aggregation semantics", {
+  source <- sf::st_sf(
+    total = c(10, NA_real_),
+    geometry = sf::st_sfc(square(0, 1), square(1, 2), crs = 3347)
+  )
+  target <- sf::st_sf(
+    geometry = sf::st_sfc(square(0.5, 1.5), crs = 3347)
+  )
+  meta <- meta_for_additive_variables("synthetic", c(total = "total"))
+
+  keep_na <- tongfen_estimate(target, source, meta, na.rm = FALSE)
+  drop_na <- tongfen_estimate(target, source, meta, na.rm = TRUE)
+
+  expect_true(is.na(keep_na$total))
+  expect_equal(drop_na$total, 5)
+})