QGIS is out, it doesn't support raster attribute tables at all.
I don't know about GRASS.
GDAL would require you to (manually or in python code) handcraft VRTs with lookup tables from the attribute table fields for each "band", then composite the VRTs.
In ArcGIS, you can only use Lookup and Composite Bands. However, look at the versions of those functions ...
Note where the error is coming from and step through your code. Its not even getting to the Savitsky-Golay line.
> v = as.vector(r)
> z = na.spline(v)
> s1.ts2 = ts(z, start=1, end=nlayers(z), frequency=12)
Error in (function (classes, fdef, mtable) :
unable to find an inherited method for function ‘nlayers’ for signature ‘"numeric"’
The error ...
If all you need to do is visually represent "number of times an area has been burnt", there's a very simple solution with the polygon layer(s)
Simply make the polygon layers all the same color, and set their opacity to 30-50%. The more overlapping polygons, the darker they appear.
Here's what a layer looks like at 30% opacity.
Set the 'burnt' area of each of your 16 rasters equal to 1 and use r.null on each of your rasters to convert 'nodata' values to 0.
To use r.null, set "The value to replace the null value by [optional]" to 0 and set the "GRASS GIS 7 region extent" to anything larger than the total extent of all 16 rasters.
Then use raster calculator to add all 16 layers. ...
You normally would ST_ReClass() the useless values of your raster to nodata (NULL). Then you can convert the withvalue pixel zones to polygons with ST_DumpAsPolygons() or you could intersect them with other polygons with ST_Clip() or ST_Intersection().