Could skipmissing help if you use missing instead of NaN? For some reason it doesn’t seem to work on your MWE but it promises to be the efficient way to filter missing values in calculations.

Maybe you are asking for a package that does this out-of-the-box, but I think using LazyArrays.@~ to create Broadcasted objects is a neat approach:

xs = [[NaN 1; 2 3], [1 NaN; 2 3]]
add!!(a, b) = a .+ b
add!!(a::AbstractArray, b) = a .+= b
sumpw(f, xs) = mapreduce(f, add!!, xs)
using LazyArrays: @~
fillnan(x) = isnan(x) ? zero(x) : x
sumpw(x -> (@~ fillnan.(x)), xs) ./
sumpw(x -> (@~ .!isnan.(x)), xs)

With this approach, you don’t have to allocate intermediate arrays like fillnan.(x). The array allocated by materializing fillnan.(xs[1]) +. fillnan.(xs[2]) is used for a entire base case of the reduction. You can also use mapfoldl instead of mapreduce to minimize the allocation.

(I’ll probably add add!! to BangBang.jl so that above code would be more generic.)

This is also easily parallelizable by using Transducers; sumpw(f, xs) = reduce(add!!, Map(f), xs).

That’s right. The reason why I thought it was neat was that @~-based implementation allocates arrays only when necessary; i.e., thanks to the clever design of the broadcasting machinery, fillnan.(x) and a .+ b are fused even though those expressions are not syntactically in the same “dot call expression.”