diff --git a/R/LCBD.comp.R b/R/LCBD.comp.R
index 6e2fc96..88a94b5 100644
--- a/R/LCBD.comp.R
+++ b/R/LCBD.comp.R
@@ -98,7 +98,8 @@
 #' col = c("white", "brown"), main = "Doubs fish abundance diff. LCBD")
 #' }
 #' 
-#' \dontrun{
+#' \donttest{
+#' if(require("ade4", quietly = TRUE) & require("betapart", quietly = TRUE)){
 #' ### Example 3
 #' ### This example requires packages \code{"betapart"} and \code{"ade4"} for data. 
 #' ### For the Baselga-family indices, the same partitioning results are obtained using
@@ -131,6 +132,7 @@
 #' out3.AbDiff$beta
 #' # Compare BDtotal here to RichDiff/Nes in out.comp$part (above)
 #' }
+#' }
 #'
 #' @importFrom stats as.dist
 #' @export LCBD.comp
diff --git a/R/Tiahura.R b/R/Tiahura.R
index c6d769d..3b0c0f0 100644
--- a/R/Tiahura.R
+++ b/R/Tiahura.R
@@ -108,7 +108,7 @@
 #'                               chron=TRUE)
 #' 
 #' ## Plotting the results
-#' par(mfrow=c(3,1))
+#' oldpar <- par(mfrow=c(3,1))
 #' 
 #' ## First graph: constrained clusters
 #' par(mar=c(3,6.5,2,2))
@@ -176,6 +176,6 @@
 #' ## plot_slice(3L,5L)
 #' ## plot_slice(4L,5L)
 #' ## plot_slice(5L,5L)
-#' 
+#' par(oldpar)
 NULL
 ##
diff --git a/R/WRperiodogram.R b/R/WRperiodogram.R
index 0f58ec7..4a26d2c 100644
--- a/R/WRperiodogram.R
+++ b/R/WRperiodogram.R
@@ -140,7 +140,7 @@
 #' # Periodogram of the lunar component only 
 #' res.moon.250 <- WRperiodogram(moon, nperm=0)  # T1=2, T2=n/2=250; no test
 #' res.moon.130 <- WRperiodogram(moon, T2=130, nperm=499)
-#' par(mfrow=c(1,2))
+#' oldpar <- par(mfrow=c(1,2))
 #' # Plot 2 moon cycles, n = 118*2 = 236 points
 #' plot(moon[1:236], xlab="One time unit = 6 hours") 
 #' plot(res.moon.130, prog=1) # Plot the periodogram
@@ -207,6 +207,7 @@
 #' # The yearly cycle and harmonics are significant
 #' # Compare the results of afc() to those of WRperiodogram above
 #' acf(ldeaths)   # lag=1.0 is one year; see ?acf
+#' par(oldpar)
 #' 
 #' @useDynLib adespatial, .registration = TRUE 
 #' @importFrom graphics lines points
diff --git a/R/aem.R b/R/aem.R
index 07e3e5f..ed97362 100755
--- a/R/aem.R
+++ b/R/aem.R
@@ -71,7 +71,7 @@
 #' # Illustrate 4 AEM eigenfunctions using bubble plots
 #' opal <- palette()
 #' palette(c("black","white"))
-#' par(mfrow=c(2,2))
+#' oldpar <- par(mfrow=c(2,2))
 #' symbols(x=xy[,2:3], circles=abs(res$vectors[,1]), inches=FALSE, asp=1,
 #'  fg=ifelse(sign(-res$vectors[,1])+1>0,1,0), 
 #'  bg=ifelse(sign(res$vectors[,1])+1>0,1,0), xlab="x", ylab="y")
@@ -113,6 +113,7 @@
 #' title("AEM 4")
 #' 
 #' palette(opal)
+#' par(oldpar)
 #' 
 #' # Construct AEM eigenfunctions with a function of the distance
 #' # as weights to put on the links
diff --git a/R/aem.build.binary.R b/R/aem.build.binary.R
index 40a119e..fb8bf7c 100755
--- a/R/aem.build.binary.R
+++ b/R/aem.build.binary.R
@@ -168,7 +168,7 @@
             xy.range.min<-xy.range[1,2]-((xy.range[2,2]-xy.range[1,2])/5)
             xy.range.max<-xy.range[2,2]
             
-            par(mar=c(1,1,1,1))
+            oldpar <- par(mar=c(1,1,1,1))
             plot(coords[, 2:3], pch = 20, asp = 1, cex = 0.5, axes = FALSE, 
                 xlab = "", ylab = "", ylim=c(xy.range.min,xy.range.max))
             segments(x0 = coords[link[, 1], 2], y0 = coords[link[, 
@@ -181,6 +181,7 @@
             points(site0[1],site0[2],pch=19,col="blue")
             segments(x0=site0[1],y0=site0[2],x1=coords[low.y,2],y1=coords[low.y,3],col="blue")
             par(mar=c(5,4,4,2))
+            par(oldpar)
         }
         n.low.y <- length(low.y)
         link <- rbind(cbind(rep(0, n.low.y), low.y), as.matrix(link))
diff --git a/R/aem.weight.edges.R b/R/aem.weight.edges.R
index 35a9160..23f04ae 100755
--- a/R/aem.weight.edges.R
+++ b/R/aem.weight.edges.R
@@ -195,7 +195,7 @@ aem.weight.edges <-
             xy.range.min<-xy.range[1,2]-((xy.range[2,2]-xy.range[1,2])/5)
             xy.range.max<-xy.range[2,2]
             
-            par(mar=c(1,1,1,1))
+            oldpar <- par(mar=c(1,1,1,1))
             plot(coords[, 2:3], pch = 20, asp = 1, cex = 0.5, axes = FALSE, xlab = "", ylab = "", ylim=c(xy.range.min,xy.range.max))
             segments(x0 = coords[links[, 1], 2], y0 = coords[links[, 1], 3], x1 = coords[links[, 2], 2], y1 = coords[links[, 2], 3], col = "red")
             text(coords[,2:3],labels=as.character(coords[,1]),pos=2,col="red")
@@ -205,6 +205,7 @@ aem.weight.edges <-
             points(site0[1],site0[2],pch=19,col="blue")
             segments(x0=site0[1],y0=site0[2],x1=coords[low.y,2],y1=coords[low.y,3],col="blue")
             par(mar=c(5,4,4,2))
+            par(oldpar)
         }
         
         ### Reorganize the edges 
diff --git a/R/chooseCN.R b/R/chooseCN.R
index 939c705..d1917ae 100644
--- a/R/chooseCN.R
+++ b/R/chooseCN.R
@@ -59,12 +59,12 @@
 #' if(require("ade4", quietly = TRUE)){
 #' data(mafragh)
 #'
-#' par(mfrow=c(2,2))
+#' oldpar <- par(mfrow=c(2,2))
 #' cn1 <- chooseCN(mafragh$xy,ask=FALSE,type=1)
 #' cn2 <- chooseCN(mafragh$xy,ask=FALSE,type=2)
 #' cn3 <- chooseCN(mafragh$xy,ask=FALSE,type=3)
 #' cn4 <- chooseCN(mafragh$xy,ask=FALSE,type=4)
-#' 
+#' par(oldpar)
 #' }
 #'
 #' @export chooseCN
diff --git a/R/constr.hclust.R b/R/constr.hclust.R
index 863b48d..22de1bb 100644
--- a/R/constr.hclust.R
+++ b/R/constr.hclust.R
@@ -277,7 +277,7 @@
 #' ## A fonction to plot the Whiskey clustering results:
 #' 
 #' plotWhiskey <- function(wh, k) {
-#'    par(fig=c(0,1,0,1))
+#'    oldpar <- par(fig=c(0,1,0,1))
 #'    plot(grpWD2cst_ScotchWhiskey[[wh]], k=k, links=TRUE, las=1,
 #'         xlab="Eastings (km)", ylab="Northings (km)", cex=0.1, lwd=3,
 #'         main=sprintf("Feature: %s",wh))
@@ -294,6 +294,7 @@
 #'    text(ScotchWhiskey$geo@coords/1000,labels=1:length(ScotchWhiskey$geo))
 #'    rect(xleft=SpeyZoom$xlim[1L], ybottom=SpeyZoom$ylim[1L],
 #'         xright=SpeyZoom$xlim[2L], ytop=SpeyZoom$ylim[2L], lwd=2, lty=1L)
+#'    par(oldpar)
 #' }
 #' 
 #' ## Plot the clustering results on the map of Scotland for 5 groups.
@@ -306,14 +307,13 @@
 #' 
 #' ## End of the Scotch Whiskey tasting data example
 #' 
-#' \dontrun{
+#' \donttest{
 #' 
 #' ## Third example: Fish community composition along the Doubs River,
 #' ## France. The sequence is analyzed as a case of chronological
 #' ## clustering, substituting space for time.
 #' 
-#' library(ade4)
-#' library(adespatial)
+#' if(require("ade4", quietly = TRUE)){
 #' data(doubs, package="ade4")
 #' Doubs.D <- dist.ldc(doubs$fish, method="hellinger")
 #' grpWD2cst_fish <- constr.hclust(Doubs.D, method="ward.D2", chron=TRUE,
@@ -349,12 +349,13 @@
 #'         coords=ex.xy      # File of geographic coordinates
 #'     )
 #' 
-#' par(mfrow=c(1,2))
+#' oldpar <- par(mfrow=c(1,2))
 #' ## Plot the map of the results for k = 3
 #' plot(test.out, k=3)
 #' ## Plot the dendrogram
 #' stats:::plot.hclust(test.out, hang=-1)
-#' 
+#' par(oldpar)
+#' }
 #'
 #' ## Same example modified: disjoint clusters
 #' ## Same ex.Y and ex.xy as in the previous example
@@ -381,12 +382,13 @@
 #'     )
 #' cutree(test.out2, k=2)
 #' 
-#' par(mfrow=c(1,2))
+#' oldpar <- par(mfrow=c(1,2))
 #' ## Plot the map of the results for k = 3
 #' plot(test.out2, k=3)
 #' ## Plot the dendrogram showing the disconnected groups
 #' stats:::plot.hclust(test.out2, hang=-1)
 #' axis(2,at=0:ceiling(max(test.out2$height,na.rm=TRUE)))
+#' par(oldpar)
 #' 
 #' ## End of the disjoint clusters example
 #' 
@@ -396,8 +398,7 @@
 #' ## Computing time grows with N at roughly the same speed as the memory
 #' ## storage requirements to store the dissimilarity matrices.
 #' 
-#' require(magrittr)
-#' require(pryr)
+#' if(require("magrittr", quietly = TRUE) & require("pryr", quietly = TRUE)){
 #' 
 #' benchmark <- function(nobj) {
 #'     # Argument -
@@ -436,7 +437,7 @@
 #' 
 #' ## Plotting the results:
 #' ok <- res %>% apply(1L, function(x) !x %>% is.na %>% any)
-#' par(mar=c(3,6,2,2),mfrow=c(2L,1L))
+#' oldpar <- par(mar=c(3,6,2,2),mfrow=c(2L,1L))
 #' barplot(height = res[ok,"Time (sec)"], names.arg= res[ok,"N.objects"],
 #'         ylab="Time (seconds)\n",xlab="",las=1L,log="y")
 #' par(mar=c(5,6,0,2))
@@ -446,9 +447,10 @@
 #' 
 #' ## Examine the output file
 #' res
-#' 
+#' par(oldpar)
 #' ## End of the benchmarking example
 #' }
+#' }
 #' ## End of examples
 #' 
 #' @useDynLib adespatial, .registration = TRUE
diff --git a/R/envspace.test.R b/R/envspace.test.R
index dd179d1..701bfa8 100644
--- a/R/envspace.test.R
+++ b/R/envspace.test.R
@@ -148,7 +148,7 @@
 #' @keywords spatial
 #' 
 #' @examples
-#' \dontrun{
+#' \donttest{
 #' if(require(vegan)) { 
 #' # Illustration of the test of the SSEF on the oribatid mite data
 #' # (Borcard et al. 1992, 1994 for details on the dataset):
diff --git a/R/listw.explore.R b/R/listw.explore.R
index d9fe827..0548284 100644
--- a/R/listw.explore.R
+++ b/R/listw.explore.R
@@ -6,7 +6,7 @@
 #' @seealso \code{\link{chooseCN}}
 #' @returns No return value
 #' @examples
-#' \dontrun{
+#' if(interactive()){
 #' ## a matrix or an object of class 'Spatial*' should be in the global environment
 #' xy <- matrix(rnorm(50), 25)
 #' listw.explore()
diff --git a/R/listw.select.R b/R/listw.select.R
index 1b7adda..ffee85d 100644
--- a/R/listw.select.R
+++ b/R/listw.select.R
@@ -125,7 +125,7 @@
 #'@keywords spatial
 #'
 #' @examples
-#' \dontrun{
+#' \donttest{
 #' if(require(spdep)) {
 #' ### Create a grid of 15 x 15:
 #' grid <- expand.grid(x = seq(1, 15, 1), y = seq(1, 15, 1))
diff --git a/R/multispati.R b/R/multispati.R
index a09064f..0d787a0 100644
--- a/R/multispati.R
+++ b/R/multispati.R
@@ -108,7 +108,7 @@
 #'       labels <- lapply(1:length(obj$eig),function(i) bquote(lambda[.(i)]))
 #'       # draw the plot
 #'       xmax <- eval(as.list(obj$call)$dudi)$eig[1]*1.1
-#'       par(las=1)
+#'       oldpar <- par(las=1)
 #'       var <- sum.obj[,2]
 #'       moran <- sum.obj[,3]
 #'       plot(x=var,y=moran,type='n',xlab='Inertia',ylab="Spatial autocorrelation (I)",
@@ -123,6 +123,7 @@
 #'       segments(0,I0,xmax,I0,lty=2)
 #'       abline(v=0)
 #'       title("Spatial and inertia components of the eigenvalues")
+#'       par(oldpar)
 #'     }
 #'     fgraph(maf.coa.ms)
 #'     ## end eigenvalues details
@@ -134,11 +135,11 @@
 #'       g4 <- s.corcircle(maf.coa.ms$as, plot = FALSE)
 #'       G1 <- ADEgS(list(g2, g3, g4), layout = c(1, 3))
 #'     } else {
-#'       par(mfrow = c(1, 3))
+#'       oldpar <- par(mfrow = c(1, 3))
 #'       barplot(maf.coa$eig)
 #'       barplot(maf.coa.ms$eig) 
 #'       s.corcircle(maf.coa.ms$as)
-#'       par(mfrow = c(1, 1))
+#'       par(oldpar)
 #'     }
 #'  
 #'  
@@ -149,12 +150,12 @@
 #'       g8 <- s.value(maf.xy, maf.coa.ms$li[, 2], plot = FALSE)
 #'       G2 <- ADEgS(list(g5, g6, g7, g8), layout = c(2, 2))
 #'     } else {
-#'       par(mfrow = c(2, 2))
+#'       oldpar <- par(mfrow = c(2, 2))
 #'       s.value(maf.xy, -maf.coa$li[, 1])
 #'       s.value(maf.xy, -maf.coa$li[, 2])
 #'       s.value(maf.xy, maf.coa.ms$li[, 1])
 #'       s.value(maf.xy, maf.coa.ms$li[, 2])
-#'       par(mfrow = c(1, 1))
+#'       par(oldpar)
 #'     }
 #' 
 #' 
@@ -167,10 +168,10 @@
 #'       g10 <- s.match(w1.ms, w1.msm, plab.cex = 0.75, plot = FALSE)
 #'       G3 <- cbindADEg(g9, g10, plot = TRUE)
 #'     } else {
-#'       par(mfrow = c(1,2))
+#'       oldpar <- par(mfrow = c(1,2))
 #'       s.match(w1, w1m, clab = 0.75)
 #'       s.match(w1.ms, w1.msm, clab = 0.75)
-#'       par(mfrow = c(1, 1))
+#'       par(oldpar)
 #'     }
 #' 
 #'     maf.pca <- dudi.pca(mafragh$env, scannf = FALSE)
diff --git a/R/plot.TBI.R b/R/plot.TBI.R
index d92e09a..ae9226b 100644
--- a/R/plot.TBI.R
+++ b/R/plot.TBI.R
@@ -109,7 +109,7 @@
 #'
 #' ## Draw BC plots
 #' 
-#' par(mfrow=c(1,2))
+#' oldpar <- par(mfrow=c(1,2))
 #' 
 #' s.names <- paste("Surv",1:12,sep=".")
 #' 
@@ -125,7 +125,7 @@
 #' 
 #' ## In the 3rd plot, draw all symbols small and of the same size, using cex.symb=NULL
 #' 
-#' dev.off()
+#' par(oldpar)
 #' 
 #' plot(res1, s.names=1:12, col.bg="gold", pch.loss=23, pch.gain=24, 
 #' main="B-C plot, Pyrifos, surveys 4 & 11", cex.symb=NULL)
@@ -145,7 +145,7 @@
 #' 
 #' ## Draw BC plots
 #' 
-#' par(mfrow=c(1,2))
+#' oldpar <- par(mfrow=c(1,2))
 #' 
 #' s.names=rownames(res2$BCD.mat)
 #' 
@@ -158,7 +158,7 @@
 #' 
 #' plot(res2, s.names=1:20, col.rim="coral2", pch.loss=19, pch.gain=15, 
 #' main="B-C plot, Mite data", xlim=c(0,0.6), ylim=c(0,0.6))
-#' 
+#' par(oldpar)
 #' }
 #' 
 #' @importFrom graphics abline
diff --git a/R/stimodels.R b/R/stimodels.R
index 43250d5..4bdabc3 100644
--- a/R/stimodels.R
+++ b/R/stimodels.R
@@ -194,7 +194,7 @@
 #' # of the times shows a significant spatial structure can be done with model 6a. For a  
 #' # global test of temporal structure at the various sites, use model 6b.
 #' 
-#' \dontrun{     # Code not run during CRAN software tests
+#' \donttest{     # Code not run during CRAN software tests
 #'
 #' # Example 2. Run space-time analysis with global tests for main effects after testing
 #' # the interaction, which is significant in this example 
diff --git a/man/LCBD.comp.Rd b/man/LCBD.comp.Rd
index 18ade67..9fc9d4d 100644
--- a/man/LCBD.comp.Rd
+++ b/man/LCBD.comp.Rd
@@ -99,7 +99,8 @@ s.value(fish.xy, out.fish.AbDiff$LCBD, method="size", symbol = "circle",
 col = c("white", "brown"), main = "Doubs fish abundance diff. LCBD")
 }
 
-\dontrun{
+\donttest{
+if(require("ade4", quietly = TRUE) & require("betapart", quietly = TRUE)){
 ### Example 3
 ### This example requires packages \code{"betapart"} and \code{"ade4"} for data. 
 ### For the Baselga-family indices, the same partitioning results are obtained using
@@ -132,6 +133,7 @@ out3.AbDiff = LCBD.comp(out3$beta.bray.gra, sqrt.D=TRUE)
 out3.AbDiff$beta
 # Compare BDtotal here to RichDiff/Nes in out.comp$part (above)
 }
+}
 
 }
 \references{
diff --git a/man/Tiahura.Rd b/man/Tiahura.Rd
index fc87052..4417f7f 100644
--- a/man/Tiahura.Rd
+++ b/man/Tiahura.Rd
@@ -97,7 +97,7 @@ tiah.chclust <- constr.hclust(tiah.jac, coords=Tiahura$habitat[,"distance"],
                               chron=TRUE)
 
 ## Plotting the results
-par(mfrow=c(3,1))
+oldpar <- par(mfrow=c(3,1))
 
 ## First graph: constrained clusters
 par(mar=c(3,6.5,2,2))
@@ -165,7 +165,7 @@ plot_slice(1L,5L)
 ## plot_slice(3L,5L)
 ## plot_slice(4L,5L)
 ## plot_slice(5L,5L)
-
+par(oldpar)
 }
 \references{
 Galzin, R. & P. Legendre. 1987. The fish communities of a coral
diff --git a/man/WRperiodogram.Rd b/man/WRperiodogram.Rd
index c3506df..7606541 100755
--- a/man/WRperiodogram.Rd
+++ b/man/WRperiodogram.Rd
@@ -152,7 +152,7 @@ tide <- periodic.component(1:n, 2, 0)
 # Periodogram of the lunar component only 
 res.moon.250 <- WRperiodogram(moon, nperm=0)  # T1=2, T2=n/2=250; no test
 res.moon.130 <- WRperiodogram(moon, T2=130, nperm=499)
-par(mfrow=c(1,2))
+oldpar <- par(mfrow=c(1,2))
 # Plot 2 moon cycles, n = 118*2 = 236 points
 plot(moon[1:236], xlab="One time unit = 6 hours") 
 plot(res.moon.130, prog=1) # Plot the periodogram
@@ -218,6 +218,7 @@ plot(ld.res.perio, prog=3) # Progressive correction for multiple tests
 # The yearly cycle and harmonics are significant
 # Compare the results of afc() to those of WRperiodogram above
 acf(ldeaths)   # lag=1.0 is one year; see ?acf
+par(oldpar)
 
 }
 \references{
diff --git a/man/aem.Rd b/man/aem.Rd
index 878301b..52f1b4e 100644
--- a/man/aem.Rd
+++ b/man/aem.Rd
@@ -75,7 +75,7 @@ res$values
 # Illustrate 4 AEM eigenfunctions using bubble plots
 opal <- palette()
 palette(c("black","white"))
-par(mfrow=c(2,2))
+oldpar <- par(mfrow=c(2,2))
 symbols(x=xy[,2:3], circles=abs(res$vectors[,1]), inches=FALSE, asp=1,
  fg=ifelse(sign(-res$vectors[,1])+1>0,1,0), 
  bg=ifelse(sign(res$vectors[,1])+1>0,1,0), xlab="x", ylab="y")
@@ -117,6 +117,7 @@ symbols(x=xy[,2:3], circles=abs(res2$vectors[,4]), inches=FALSE,asp=1,
 title("AEM 4")
 
 palette(opal)
+par(oldpar)
 
 # Construct AEM eigenfunctions with a function of the distance
 # as weights to put on the links
diff --git a/man/chooseCN.Rd b/man/chooseCN.Rd
index 5c40cef..215b126 100644
--- a/man/chooseCN.Rd
+++ b/man/chooseCN.Rd
@@ -86,12 +86,12 @@ using 'graph' as keyword.
 if(require("ade4", quietly = TRUE)){
 data(mafragh)
 
-par(mfrow=c(2,2))
+oldpar <- par(mfrow=c(2,2))
 cn1 <- chooseCN(mafragh$xy,ask=FALSE,type=1)
 cn2 <- chooseCN(mafragh$xy,ask=FALSE,type=2)
 cn3 <- chooseCN(mafragh$xy,ask=FALSE,type=3)
 cn4 <- chooseCN(mafragh$xy,ask=FALSE,type=4)
-
+par(oldpar)
 }
 
 }
diff --git a/man/constr.hclust-class.Rd b/man/constr.hclust-class.Rd
index cb33db3..6963c44 100644
--- a/man/constr.hclust-class.Rd
+++ b/man/constr.hclust-class.Rd
@@ -39,6 +39,9 @@ to the function).}
 \item{coords}{A copy of the coordinates (if coordinates were provided to the
 function).} }
 }
+\value{
+No return value
+}
 \description{
 Files belonging to this class hold information about the constrained
 agglomerative clustering and allows one to display results graphically.
diff --git a/man/constr.hclust.Rd b/man/constr.hclust.Rd
index 06d8770..d926a68 100644
--- a/man/constr.hclust.Rd
+++ b/man/constr.hclust.Rd
@@ -244,7 +244,7 @@ grpWD2cst_ScotchWhiskey[["norm"]] <-
 ## A fonction to plot the Whiskey clustering results:
 
 plotWhiskey <- function(wh, k) {
-   par(fig=c(0,1,0,1))
+   oldpar <- par(fig=c(0,1,0,1))
    plot(grpWD2cst_ScotchWhiskey[[wh]], k=k, links=TRUE, las=1,
         xlab="Eastings (km)", ylab="Northings (km)", cex=0.1, lwd=3,
         main=sprintf("Feature: \%s",wh))
@@ -261,6 +261,7 @@ plotWhiskey <- function(wh, k) {
    text(ScotchWhiskey$geo@coords/1000,labels=1:length(ScotchWhiskey$geo))
    rect(xleft=SpeyZoom$xlim[1L], ybottom=SpeyZoom$ylim[1L],
         xright=SpeyZoom$xlim[2L], ytop=SpeyZoom$ylim[2L], lwd=2, lty=1L)
+   par(oldpar)
 }
 
 ## Plot the clustering results on the map of Scotland for 5 groups.
@@ -273,14 +274,13 @@ plotWhiskey("norm", 5L)
 
 ## End of the Scotch Whiskey tasting data example
 
-\dontrun{
+\donttest{
 
 ## Third example: Fish community composition along the Doubs River,
 ## France. The sequence is analyzed as a case of chronological
 ## clustering, substituting space for time.
 
-library(ade4)
-library(adespatial)
+if(require("ade4", quietly = TRUE)){
 data(doubs, package="ade4")
 Doubs.D <- dist.ldc(doubs$fish, method="hellinger")
 grpWD2cst_fish <- constr.hclust(Doubs.D, method="ward.D2", chron=TRUE,
@@ -316,12 +316,13 @@ test.out <-
         coords=ex.xy      # File of geographic coordinates
     )
 
-par(mfrow=c(1,2))
+oldpar <- par(mfrow=c(1,2))
 ## Plot the map of the results for k = 3
 plot(test.out, k=3)
 ## Plot the dendrogram
 stats:::plot.hclust(test.out, hang=-1)
-
+par(oldpar)
+}
 
 ## Same example modified: disjoint clusters
 ## Same ex.Y and ex.xy as in the previous example
@@ -348,12 +349,13 @@ test.out2 <-
     )
 cutree(test.out2, k=2)
 
-par(mfrow=c(1,2))
+oldpar <- par(mfrow=c(1,2))
 ## Plot the map of the results for k = 3
 plot(test.out2, k=3)
 ## Plot the dendrogram showing the disconnected groups
 stats:::plot.hclust(test.out2, hang=-1)
 axis(2,at=0:ceiling(max(test.out2$height,na.rm=TRUE)))
+par(oldpar)
 
 ## End of the disjoint clusters example
 
@@ -363,8 +365,7 @@ axis(2,at=0:ceiling(max(test.out2$height,na.rm=TRUE)))
 ## Computing time grows with N at roughly the same speed as the memory
 ## storage requirements to store the dissimilarity matrices.
 
-require(magrittr)
-require(pryr)
+if(require("magrittr", quietly = TRUE) & require("pryr", quietly = TRUE)){
 
 benchmark <- function(nobj) {
     # Argument -
@@ -403,7 +404,7 @@ res <- benchmark(nobj=c(1000,2000,5000,10000,20000,50000,100000))
 
 ## Plotting the results:
 ok <- res \%>\% apply(1L, function(x) !x \%>\% is.na \%>\% any)
-par(mar=c(3,6,2,2),mfrow=c(2L,1L))
+oldpar <- par(mar=c(3,6,2,2),mfrow=c(2L,1L))
 barplot(height = res[ok,"Time (sec)"], names.arg= res[ok,"N.objects"],
         ylab="Time (seconds)\n",xlab="",las=1L,log="y")
 par(mar=c(5,6,0,2))
@@ -413,9 +414,10 @@ barplot(height = res[ok,"Storage (MiB)"], names.arg= res[ok,"N.objects"],
 
 ## Examine the output file
 res
-
+par(oldpar)
 ## End of the benchmarking example
 }
+}
 ## End of examples
 
 }
diff --git a/man/envspace.test.Rd b/man/envspace.test.Rd
index 3ad9f13..ce3e7ac 100644
--- a/man/envspace.test.Rd
+++ b/man/envspace.test.Rd
@@ -140,7 +140,7 @@ environmental variables. This does not exclude that this effect may exist with r
 to some unmeasured variables.
 }
 \examples{
-\dontrun{
+\donttest{
 if(require(vegan)) { 
 # Illustration of the test of the SSEF on the oribatid mite data
 # (Borcard et al. 1992, 1994 for details on the dataset):
diff --git a/man/listw.explore.Rd b/man/listw.explore.Rd
index 0ed0336..523f9b9 100644
--- a/man/listw.explore.Rd
+++ b/man/listw.explore.Rd
@@ -6,11 +6,14 @@
 \usage{
 listw.explore()
 }
+\value{
+No return value
+}
 \description{
 Interactive tool to generate R code that creates a spatial weighting matrix
 }
 \examples{
-\dontrun{
+if(interactive()){
 ## a matrix or an object of class 'Spatial*' should be in the global environment
 xy <- matrix(rnorm(50), 25)
 listw.explore()
diff --git a/man/listw.select.Rd b/man/listw.select.Rd
index 208afbe..dd0205d 100644
--- a/man/listw.select.Rd
+++ b/man/listw.select.Rd
@@ -131,7 +131,7 @@ While the selection of the SWM is the most critical step of the
  negative spatial predictors.
 }
 \examples{
-\dontrun{
+\donttest{
 if(require(spdep)) {
 ### Create a grid of 15 x 15:
 grid <- expand.grid(x = seq(1, 15, 1), y = seq(1, 15, 1))
diff --git a/man/multispati.Rd b/man/multispati.Rd
index 4bb2a92..2c3484d 100644
--- a/man/multispati.Rd
+++ b/man/multispati.Rd
@@ -113,7 +113,7 @@ if (require(spdep, quiet = TRUE) & require(ade4, quiet = TRUE)) {
       labels <- lapply(1:length(obj$eig),function(i) bquote(lambda[.(i)]))
       # draw the plot
       xmax <- eval(as.list(obj$call)$dudi)$eig[1]*1.1
-      par(las=1)
+      oldpar <- par(las=1)
       var <- sum.obj[,2]
       moran <- sum.obj[,3]
       plot(x=var,y=moran,type='n',xlab='Inertia',ylab="Spatial autocorrelation (I)",
@@ -128,6 +128,7 @@ if (require(spdep, quiet = TRUE) & require(ade4, quiet = TRUE)) {
       segments(0,I0,xmax,I0,lty=2)
       abline(v=0)
       title("Spatial and inertia components of the eigenvalues")
+      par(oldpar)
     }
     fgraph(maf.coa.ms)
     ## end eigenvalues details
@@ -139,11 +140,11 @@ if (require(spdep, quiet = TRUE) & require(ade4, quiet = TRUE)) {
       g4 <- s.corcircle(maf.coa.ms$as, plot = FALSE)
       G1 <- ADEgS(list(g2, g3, g4), layout = c(1, 3))
     } else {
-      par(mfrow = c(1, 3))
+      oldpar <- par(mfrow = c(1, 3))
       barplot(maf.coa$eig)
       barplot(maf.coa.ms$eig) 
       s.corcircle(maf.coa.ms$as)
-      par(mfrow = c(1, 1))
+      par(oldpar)
     }
  
  
@@ -154,12 +155,12 @@ if (require(spdep, quiet = TRUE) & require(ade4, quiet = TRUE)) {
       g8 <- s.value(maf.xy, maf.coa.ms$li[, 2], plot = FALSE)
       G2 <- ADEgS(list(g5, g6, g7, g8), layout = c(2, 2))
     } else {
-      par(mfrow = c(2, 2))
+      oldpar <- par(mfrow = c(2, 2))
       s.value(maf.xy, -maf.coa$li[, 1])
       s.value(maf.xy, -maf.coa$li[, 2])
       s.value(maf.xy, maf.coa.ms$li[, 1])
       s.value(maf.xy, maf.coa.ms$li[, 2])
-      par(mfrow = c(1, 1))
+      par(oldpar)
     }
 
 
@@ -172,10 +173,10 @@ if (require(spdep, quiet = TRUE) & require(ade4, quiet = TRUE)) {
       g10 <- s.match(w1.ms, w1.msm, plab.cex = 0.75, plot = FALSE)
       G3 <- cbindADEg(g9, g10, plot = TRUE)
     } else {
-      par(mfrow = c(1,2))
+      oldpar <- par(mfrow = c(1,2))
       s.match(w1, w1m, clab = 0.75)
       s.match(w1.ms, w1.msm, clab = 0.75)
-      par(mfrow = c(1, 1))
+      par(oldpar)
     }
 
     maf.pca <- dudi.pca(mafragh$env, scannf = FALSE)
diff --git a/man/plot.TBI.Rd b/man/plot.TBI.Rd
index 7039033..d064b53 100644
--- a/man/plot.TBI.Rd
+++ b/man/plot.TBI.Rd
@@ -125,7 +125,7 @@ res1$BCD.mat
 
 ## Draw BC plots
 
-par(mfrow=c(1,2))
+oldpar <- par(mfrow=c(1,2))
 
 s.names <- paste("Surv",1:12,sep=".")
 
@@ -141,7 +141,7 @@ main="B-C plot, Pyrifos, surveys 4 & 11", xlim=c(0,0.5), ylim=c(0.1,0.6))
 
 ## In the 3rd plot, draw all symbols small and of the same size, using cex.symb=NULL
 
-dev.off()
+par(oldpar)
 
 plot(res1, s.names=1:12, col.bg="gold", pch.loss=23, pch.gain=24, 
 main="B-C plot, Pyrifos, surveys 4 & 11", cex.symb=NULL)
@@ -161,7 +161,7 @@ res2$BCD.mat
 
 ## Draw BC plots
 
-par(mfrow=c(1,2))
+oldpar <- par(mfrow=c(1,2))
 
 s.names=rownames(res2$BCD.mat)
 
@@ -174,7 +174,7 @@ main="B-C plot, Mite data")
 
 plot(res2, s.names=1:20, col.rim="coral2", pch.loss=19, pch.gain=15, 
 main="B-C plot, Mite data", xlim=c(0,0.6), ylim=c(0,0.6))
-
+par(oldpar)
 }
 
 }
diff --git a/man/plot.constr.hclust.Rd b/man/plot.constr.hclust.Rd
index d4e17af..9adb3b0 100644
--- a/man/plot.constr.hclust.Rd
+++ b/man/plot.constr.hclust.Rd
@@ -60,6 +60,9 @@ details).}
 
 \item{...}{Other \link{graphical parameters}.}
 }
+\value{
+No return value
+}
 \description{
 Method \code{plot.constr.hclust} displays the results of space-constrained or
 time-constrained agglomerative cluster analyses obtained from multivariate
diff --git a/man/stimodels.Rd b/man/stimodels.Rd
index 3760673..d9c0b9e 100644
--- a/man/stimodels.Rd
+++ b/man/stimodels.Rd
@@ -220,7 +220,7 @@ out.1 <- stimodels(trich.log, S=22, Ti=10, model="5", nperm=99)
 # of the times shows a significant spatial structure can be done with model 6a. For a  
 # global test of temporal structure at the various sites, use model 6b.
 
-\dontrun{     # Code not run during CRAN software tests
+\donttest{     # Code not run during CRAN software tests
 
 # Example 2. Run space-time analysis with global tests for main effects after testing
 # the interaction, which is significant in this example 
diff --git a/vignettes/tutorial.Rmd b/vignettes/tutorial.Rmd
index 32c71cf..49deade 100644
--- a/vignettes/tutorial.Rmd
+++ b/vignettes/tutorial.Rmd
@@ -311,9 +311,10 @@ names(attributes(mem.gab))
 
 The eigenvalues associated to MEMs are stored in the attribute called `values`:
 ```{r, echo = -1}
-par(mar = c(0, 2, 3, 0))
+oldpar <- par(mar = c(0, 2, 3, 0))
     barplot(attr(mem.gab, "values"), 
         main = "Eigenvalues of the spatial weighting matrix", cex.main = 0.7)
+par(oldpar)
 ```
 
 A `plot` method is provided to represent MEMs. By default, eigenvectors are represented as a table (sites as rows, MEMs as columns). This representation is usually not informative and it is better to map MEMs in the geographical space by documenting the argument `SpORcoords`: