library(tmap)

tm = tm_shape(NLD_dist) +
  tm_polygons("edu_appl_sci",
              fill.scale = tm_scale_intervals(values = "pu_gn", style = "kmeans", n = 7))

# default
tm +
  tm_facets_wrap(by = "province", byrow = TRUE)

# with byrow = FALSE -- not working
tm +
  tm_facets_wrap(by = "province", byrow = FALSE)