As this answer points out, if you use sf, you can use st_area to get the area of each polygon, then order the polygons based on area with forcats::fct_reorder. Something like this should work:
mutate(area = st_area(geometry)) %>%
mutate(region = region %>% forcats::fct_reorder(-area)) %>%
Working with an outline of Pakistan from the GADM data:
p = raster::getData("GADM",country="PAK",level=0)
p is a SpatialPolygons data frame in lat-long coordinates. You should convert this to a proper equal-area projection but EPSG:3857 should be close enough:
pp = sp::spTransform(p, "+init=epsg:3857")
Now chop using the code you pasted in your question. ...