diff --git a/NEWS.md b/NEWS.md
index 2dd1570b127f4cf597633b8adbf810e47e82d485..6827e8d0d2a708d26a2255caa36a8535c2b960cd 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,6 +2,7 @@
 - Change Regression() parameter 'time_dim' to 'reg_dim', and enable the inputs to be vectors.
 - Change Trend() parameter 'time_dim' default value from 'sdate' to 'ftime'.
 - Change the default of Season() parameter 'time_dim' from 'sdate' to 'ftime'.
+- Bugfix for Regression() na.action. 'na.action = na.fail' was not functional before.
 - Add p-value by ANOVA in Trend().
 - Change MeanDims() na.rm default to FALSE to be in line with mean()
 - Remove unecessary parameter checks in Clim()
diff --git a/R/Regression.R b/R/Regression.R
index 8b2bd9c3c73d4b90122f0f6623957882ddd1c410..8e5d8af41cb26c09ef37bff974d7584d06067a8f 100644
--- a/R/Regression.R
+++ b/R/Regression.R
@@ -208,10 +208,7 @@ Regression <- function(datay, datax, reg_dim = 'sdate', formula = y ~ x,
     check_na <- TRUE
   }
 
-  # remove NAs for potential poly()
-  x2 <- x[!is.na(NApos)]
-  y2 <- y[!is.na(NApos)]
-  lm.out <- lm(formula, data = data.frame(x = x2, y = y2), na.action = na.action)
+  lm.out <- lm(formula, data = data.frame(x = x, y = y), na.action = na.action)
   coeff <- lm.out$coefficients
   if (conf) {
     conf.lower <- confint(lm.out, level = conf.lev)[, 1]
diff --git a/tests/testthat/test-Regression.R b/tests/testthat/test-Regression.R
index 85e57822d74f1c861dea5c0952268501ebe8e8f9..a29076f5ed9be19cb0cb7b776fc9d7c2344ac349 100644
--- a/tests/testthat/test-Regression.R
+++ b/tests/testthat/test-Regression.R
@@ -24,7 +24,7 @@ context("s2dv::Regression tests")
                  dim = c(date = 5, ftime = 2, lon = 2, lat = 4))
   set.seed(2)
   datax3 <- array(c(1:80) + rnorm(80),
-                 dim = c(date = 5, lon = 2, lat = 4,ftime = 2))
+                 dim = c(date = 5, lon = 2, lat = 4, ftime = 2))
 
 
 ##############################################
@@ -158,17 +158,12 @@ test_that("3. Output checks: dat2", {
     2
   )
   expect_equal(
-    length(which(is.na(Regression(datay2, datax2, na.action = na.pass)$p.val))),
+    length(which(is.na(Regression(datay2, datax2, na.action = na.omit)$p.val))),
     0
   )
   expect_equal(
     which(is.na(Regression(datay2, datax2, na.action = 1)$p.val)),
-    c(3,15)
-  )
-  expect_equal(
-    which(is.na(Regression(datay2, datax2, na.action = 1, 
-                           formula = y~poly(x, 2, raw = T))$p.val)),
-    c(3,15)
+    c(3, 15)
   )
 })