Support for NaNStatistics

It would be very handy if there was support for NaNStatistics. I couldn't find an easy way to rewrap the output array so I did this ugly workaround:

NaNStatistics.nansum(A::KeyedArray; dim=:, dims=:) = rewrap_keyedarray(nansum(A.data.data; dim, dims), A, dims isa Colon ? dim : dims)

function rewrap_keyedarray(new, original, dims)
    axs = NamedTuple{propertynames(original)}(original.keys) |> pairs
    if ndims(new) < ndims(original)
        axs = collect(axs)[1:end .!= dims]
    else
        axs = Dict(axs)
        for dim in dims
            axs[propertynames(original)[dim]] = axes(new, dim)
        end
        axs = [k => axs[k] for k in propertynames(original)]
    end
    KeyedArray(new; axs...)
end

Thanks

Start by reproducing the NaNStatistics.nansum example from the issue and inspect how AxisKeys currently wraps array operations and preserves KeyedArray axes. Done means NaNStatistics operations work on KeyedArray inputs with the expected keys and dimensions, without requiring the shown rewrap_keyedarray workaround.