Nice.
An alternative to the 3rd dimension could be to use this tip to use a tuple shield to broadcast over the array of matrices:
using DataFrames, Statistics
dfn = [DataFrame(i*ones(4,3),:auto) for i in 1:4]
m₀, = mean.((map(Matrix, dfn),))
σₘ, = std.((map(Matrix, dfn),))