diff --git a/NAMESPACE b/NAMESPACE
index 08b5091..8da178b 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -30,6 +30,7 @@ importFrom(splines, interpSpline)
 importFrom(LearnBayes, rwmetrop)
 importFrom(nlme, fdHess)
 importFrom(multcomp, glht)
+importFrom("utils", "packageVersion")
 
 ##################
 # ML SE
diff --git a/R/predict.sarlm.R b/R/predict.sarlm.R
index 626e2b3..e3bd1d5 100644
--- a/R/predict.sarlm.R
+++ b/R/predict.sarlm.R
@@ -193,7 +193,11 @@ predict.Sarlm <- function(object, newdata=NULL, listw=NULL, pred.type="TS", all.
           W <- as(listw, "CsparseMatrix")
           W <- W[region.id, region.id]
           style <- listw$style
-          listw <- mat2listw(W, row.names = region.id, style = style, zero.policy=zero.policy) # re-normalize to keep the style
+          if (packageVersion("spdep") >= "1.3.1") {
+              listw <- mat2listw(W, row.names = region.id, style = style, zero.policy=zero.policy) # re-normalize to keep the style
+          } else {
+              listw <- mat2listw(W, row.names = region.id, style = style) # re-normalize to keep the style
+          }
           rm(W) # avoid the use of a wrong W
         }
       }
@@ -247,7 +251,11 @@ predict.Sarlm <- function(object, newdata=NULL, listw=NULL, pred.type="TS", all.
           W <- as(listw.mixed, "CsparseMatrix")
           W <- W[region.id.mixed, region.id.mixed]
           style <- listw.mixed$style
-          listw.mixed <- mat2listw(W, row.names = region.id.mixed, style = style, zero.policy=zero.policy) # re-normalize to keep the style
+          if (packageVersion("spdep") >= "1.3.1") {
+              listw.mixed <- mat2listw(W, row.names = region.id.mixed, style = style, zero.policy=zero.policy) # re-normalize to keep the style
+          } else {
+              listw.mixed <- mat2listw(W, row.names = region.id.mixed, style = style) # re-normalize to keep the style
+          }
           rm(W) # avoid the use of a wrong W
         }
       }
@@ -540,7 +548,11 @@ predict.Sarlm <- function(object, newdata=NULL, listw=NULL, pred.type="TS", all.
         for (i in 1:nrow(newdata)) {
           region.id.temp <- c(region.id.data, region.id.newdata[i])
           Wi <- W[region.id.temp, region.id.temp]
-          listwi <- mat2listw(Wi, row.names = region.id.temp, style = style, zero.policy=zero.policy) # re-normalize
+          if (packageVersion("spdep") >= "1.3.1") {
+              listwi <- mat2listw(Wi, row.names = region.id.temp, style = style, zero.policy=zero.policy) # re-normalize
+          } else {
+              listwi <- mat2listw(Wi, row.names = region.id.temp, style = style) # re-normalize to keep the style
+          }
           if (power)
             Wi <- as(listwi, "CsparseMatrix")
           Xi <- rbind(Xs, Xo[i,])
@@ -564,7 +576,11 @@ predict.Sarlm <- function(object, newdata=NULL, listw=NULL, pred.type="TS", all.
         for (i in 1:nrow(newdata)) {
           region.id.temp <- c(region.id.data, region.id.newdata[i])
           Wi <- W[region.id.temp, region.id.temp]
-          listwi <- mat2listw(Wi, row.names = region.id.temp, style = style, zero.policy=zero.policy) # re-normalize
+          if (packageVersion("spdep") >= "1.3.1") {
+              listwi <- mat2listw(Wi, row.names = region.id.temp, style = style, zero.policy=zero.policy) # re-normalize
+          } else {
+              listwi <- mat2listw(Wi, row.names = region.id.temp, style = style) # re-normalize to keep the style
+          }
           Wi <- as(listwi, "CsparseMatrix")
           Xi <- rbind(Xs, Xo[i,])
           # compute TC1 for S and o units
@@ -599,7 +615,11 @@ predict.Sarlm <- function(object, newdata=NULL, listw=NULL, pred.type="TS", all.
         for (i in 1:nrow(newdata)) {
           region.id.temp <- c(region.id.data, region.id.newdata[i])
           Wi <- W[region.id.temp, region.id.temp]
-          listwi <- mat2listw(Wi, row.names = region.id.temp, style = style, zero.policy=zero.policy) # re-normalize
+          if (packageVersion("spdep") >= "1.3.1") {
+              listwi <- mat2listw(Wi, row.names = region.id.temp, style = style, zero.policy=zero.policy) # re-normalize
+          } else {
+              listwi <- mat2listw(Wi, row.names = region.id.temp, style = style) # re-normalize to keep the style
+          }
           Wi <- as(listwi, "CsparseMatrix")
           Xi <- rbind(Xs, Xo[i,])
           is.data <- 1:length(ys)
@@ -633,7 +653,11 @@ predict.Sarlm <- function(object, newdata=NULL, listw=NULL, pred.type="TS", all.
         for (i in 1:nrow(newdata)) {
           region.id.temp <- c(region.id.data, region.id.newdata[i])
           Wi <- W[region.id.temp, region.id.temp]
-          listwi <- mat2listw(Wi, row.names = region.id.temp, style = style, zero.policy=zero.policy) # re-normalize
+          if (packageVersion("spdep") >= "1.3.1") {
+              listwi <- mat2listw(Wi, row.names = region.id.temp, style = style, zero.policy=zero.policy) # re-normalize
+          } else {
+              listwi <- mat2listw(Wi, row.names = region.id.temp, style = style) # re-normalize to keep the style
+          }
           Wi <- as(listwi, "CsparseMatrix")
           Xi <- rbind(Xs, Xo[i,])
           # compute TC1 for S and o units
diff --git a/vignettes/sids_models.Rmd b/vignettes/sids_models.Rmd
index 6267e3d..a431c80 100644
--- a/vignettes/sids_models.Rmd
+++ b/vignettes/sids_models.Rmd
@@ -97,7 +97,11 @@ n <- nc$BIR74
 el1 <- min(dij)/dij
 el2 <- sqrt(n[sids.nhbr$to]/n[sids.nhbr$from])
 sids.nhbr$weights <- el1*el2
-sids.nhbr.listw <- sn2listw(sids.nhbr, style="B", zero.policy=TRUE)
+if (packageVersion("spdep") >= "1.3.1") {
+  sids.nhbr.listw <- sn2listw(sids.nhbr, style="B", zero.policy=TRUE)
+} else {
+  sids.nhbr.listw <- sn2listw(sids.nhbr)
+}
 ```
 
 The first model (I) is a null model with just an intercept, the second