Compute weighted mean from a dataframe?

[floswald]

i can only compute a weighted mean from a dataframe object after converting the column x to a float array. I don't understand why, since DataArray <: AbstractArray is true?

julia> df = DataFrame(x = rand(3),w=rand(3))
3x2 DataFrame
| Row | x        | w         |
|-----|----------|-----------|
| 1   | 0.226874 | 0.266516  |
| 2   | 0.346735 | 0.57437   |
| 3   | 0.232796 | 0.0512731 |

julia> mean(df[:x],WeightVec(df[:w]))
ERROR: `start` has no method matching start(::WeightVec{Float64,DataArray{Float64,1}})
 in reduced_dims at reducedim.jl:17
 in mean at /Users/florianoswald/.julia/v0.3/DataArrays/src/reducedim.jl:335

julia> mean(array(df[:x]),WeightVec(df[:w]))
0.30438070685782537

julia> methods(mean)
...
mean(v::AbstractArray{T,N},w::WeightVec{W<:Real,Vec<:AbstractArray{T<:Real,1}}) at /Users/florianoswald/.julia/v0.3/StatsBase/src/weights.jl:234
...

julia> typeof(df[:x])
DataArray{Float64,1} (constructor with 1 method)

julia> DataArray <: AbstractArray
true

[johnmyleswhite]

DataArray does inherit from AbstractArray, but it's probably more precise to say that DataArray{T} is a kind of AbstractArray{Union(T, NAtype)}.

As such, you should expect that none of the functions in StatsBase work on DataArrays. Those that do only work because we've redefined the functions to behave differently when given DataArrays as inputs.

[floswald]

ok, I see. I guess that's not a design choice rather than something that's being implemented? I mean DataFrame being the prime data vehicle, it's clearly desirable to be able to compute some stats from it.

If that's true, can you give me a quick hint what needed to be done here? would one extend the Base.mean method overloaded in DataArrays to use wsum(a::DataArray,w::WeightVec;skipna::Bool), similar to what is done in StatsBase?

Base.mean(a::DataArray; skipna::Bool=false) =
    sum(a; skipna=skipna) / (length(a.na)-(skipna ? countnz(a.na) : 0))
Base.mean(a::DataArray,w::WeightVec; skipna::Bool=false) =
    wsum(a,w; skipna=skipna) / sum(w;skipna=skipna)

[simonster]

I think the method in StatsBase would actually work (albeit slowly and without support for skipna) if it were the method being dispatched, but mean{T}(A::DataArray{T,N},region) sorts before mean(v::AbstractArray{T,N},w::WeightVec{W<:Real,Vec<:AbstractArray{T<:Real,1}}). We could restrict the type of region, although mean in Base has no such restriction.

Computing mean quickly with a WeightVec is easy if you don't need skipna:

Base.mean(a::DataArray,w::WeightVec) =
    anyna(a) ? NA : mean(a.data, w)
Base.mean{V<:DataArray}(a::DataArray,w::WeightVec{W,V}) =
    anyna(a) || anyna(w.values) ? NA : wsum(a.data, w.values) / w.sum

(caveat: code is untested)

skipna is harder. There is not presently a variant of wsum for DataArrays, so you'd have to start by implementing that. There are several possible ways to do this that have different tradeoffs between performance, accuracy, and code complexity. A simple and accurate way would be to compute a.*w.values and sum it, although that allocates another vector where it isn't necessary. The fastest way would be to perform the element-wise summation and multiplication in a loop. In that case the question is how best to determine which values should be skipped due to NAs. The best compromise between speed and accuracy would be to use an algorithm similar to that used by sum for DataArrays, which is pairwise summation with a tweak to handle NAs, but this is also the most complicated option and wsum in StatsBase doesn't presently bother with pairwise summation.

And then there is the weighted mean across a dimension, which is another level of complexity.

[floswald]

yeah that Base.mean extension would work for me currently. Of course the whole point of DataArrays is the skipna switch, so i guess at some point wsum should take a DataArray. not sure how a.*w.values would work though, you have to subset somehow to the non - NA values there as well right?
so how do you guys feel about adding this in a PR

Base.mean{V<:DataArray}(a::DataArray,w::WeightVec{W,V}) =
    anyna(a) || anyna(w.values) ? NA : wsum(a.data, w.values) / w.sum

last question: how do you find out the order in which methods are looked up?

[johnmyleswhite]

What happens if you have non-zero weight on a missing value?

[simonster]

I think we'd want something like:

function Base.mean(a::DataArray, w::WeightVec; skipna::Bool=false)
    if skipna
        v = a .* w.values
        sum(v; skipna=true) / sum(DataArray(w.values, v.na); skipna=true)
    else
        anyna(a) ? NA : mean(a.data, w)
    end
end

function Base.mean{V<:DataArray}(a::DataArray, w::WeightVec{W,V}; skipna::Bool=false)
    if skipna
        v = a .* w.values
        sum(v; skipna=true) / sum(DataArray(w.values.data, v.na); skipna=true)
    else
        anyna(a) || anyna(w.values) ? NA : wsum(a.data, w.values.data) / w.sum
    end
end

which should be able to handle non-zero weights on missing values.

[floswald]

yes that looks right to me. the assumption is that an NA implicitly gets reweighted to 0, which makes sense I think, and the sum(DataArray(w.values.data, v.na); skipna=true) accounts for the changed sum of weights.

[simonster]

Fixed by JuliaStats/DataArrays.jl#121