Draft of fweights, pweights, aweights for glms

src/GLM.jl

@@ -14,6 +14,7 @@ module GLM
     import Base: (\), cholfact, convert, cor, show, size
     import StatsBase: coef, coeftable, confint, deviance, nulldeviance, dof, dof_residual, loglikelihood, nullloglikelihood, nobs, stderr, vcov, residuals, predict, fit, model_response, r2, r², adjr2, adjr², PValue
     import StatsFuns: xlogy
+    import StatsBase: FrequencyWeights, ProbabilityWeights, AnalyticWeights
     export coef, coeftable, confint, deviance, nulldeviance, dof, dof_residual, loglikelihood, nobs, stderr, vcov, residuals, predict, fit, fit!, model_response, r2, r², adjr2, adjr²
 
     export                              # types

src/glmfit.jl

@@ -1,3 +1,4 @@
+# TODO: expand weights to include fweights, pweights, aweights
 """
     GlmResp
 
@@ -17,13 +18,19 @@ immutable GlmResp{V<:FPVector,D<:UnivariateDistribution,L<:Link} <: ModResp
     offset::V
     "`wts:` prior case weights.  Can be of length 0."
     wts::V
+
+    fweights::FrequencyWeights
+    pweights::ProbabilityWeights
+    aweights::AnalyticWeights
+
     "`wrkwt`: working case weights for the Iteratively Reweighted Least Squares (IRLS) algorithm"
     wrkwt::V
     "`wrkresid`: working residuals for IRLS"
     wrkresid::V
 end
 
-function GlmResp{V<:FPVector, D, L}(y::V, d::D, l::L, η::V, μ::V, off::V, wts::V)
+function GlmResp{V<:FPVector, D, L}(y::V, d::D, l::L, η::V, μ::V, off::V, wts::V,
+    fweights::FrequencyWeights, pweights::ProbabilityWeights, aweights::AnalyticWeights)
     if d == Binomial()
         for yy in y
             0. <= yy <= 1. || throw(ArgumentError("$yy in y is not in [0,1]"))
@@ -38,7 +45,10 @@ function GlmResp{V<:FPVector, D, L}(y::V, d::D, l::L, η::V, μ::V, off::V, wts:
         "lengths of η, μ, y and wts ($nη, $nμ, $(length(wts)), $n) are not equal"))
     lo = length(off)
     lo == 0 || lo == n || error("offset must have length $n or length 0")
-    res = GlmResp{V,D,L}(y, d, similar(y), η, μ, off, wts, similar(y), similar(y))
+    res = GlmResp{V,D,L}(y, d, similar(y), η, μ, off,
+                         normalizeWeights(fweights, pweights, aweights),
+                         similar(y), similar(y),
+                         fweights, pweights, aweights)
     updateμ!(res, η)
     res
 end

src/glmtools.jl

@@ -98,3 +98,22 @@ loglik_obs(::Binomial, y, μ, wt, ϕ) = logpdf(Binomial(Int(wt), μ), Int(y*wt))
 loglik_obs(::Gamma, y, μ, wt, ϕ) = wt*logpdf(Gamma(1/ϕ, μ*ϕ), y)
 loglik_obs(::Normal, y, μ, wt, ϕ) = wt*logpdf(Normal(μ, sqrt(ϕ)), y)
 loglik_obs(::Poisson, y, μ, wt, ϕ) = wt*logpdf(Poisson(μ), y)
+
+# TODO: combine fweights, pweights, aweights
+function normalizeWeights(fweights::FrequencyWeights, pweights::ProbabilityWeights,
+                          aweights::AnalyticWeights)
+    nf = length(fweights)
+    np = length(pweights)
+    na = length(aweights)
+    if (nf != 0 && np == 0 && na == 0)
+        return fweights |> collect
+    elseif (nf == 0 && np != 0 && na == 0)
+        return pweights |> collect
+    elseif (nf == 0 && np == 0 && na != 0)
+        return aweights |> collect
+    elseif (nf != 0 && np != 0)
+        return (fweights .* pweights) |> collect
+    else
+        return zeros(Float64, 0)
+    end
+end

src/linpred.jl

@@ -116,6 +116,8 @@ Base.cholfact!{T}(p::SparsePredChol{T}) = p.chol
 
 invchol(x::DensePred) = inv(cholfact!(x))
 invchol(x::SparsePredChol) = cholfact!(x) \ eye(size(x.X, 2))
+# TODO: support fweights, pweights, and aweights
+# This vcov represents σ^2  X'X
 vcov(x::LinPredModel) = scale!(invchol(x.pp), dispersion(x, true))
 
 function cor(x::LinPredModel)
@@ -150,10 +152,11 @@ For linear and generalized linear models, returns the number of rows, or,
 when prior weights are specified, the sum of weights.
 """
 function nobs(obj::LinPredModel)
-    if isempty(obj.rr.wts)
+    if (length(obj.rr.fweights) == 0)
         oftype(sum(one(eltype(obj.rr.wts))), length(obj.rr.y))
     else
-        sum(obj.rr.wts)
+        sum(obj.rr.fweights) # This is the behavior for fweights
+        # TODO: add support for pweights and aweights
     end
 end
 

src/lm.jl

@@ -1,20 +1,35 @@
+# TODO: expand weights to include fweights, pweights, aweights
+# Example LmResp{Array{Float64,1}}(randn(5), randn(5), randn(5), randn(5))
 type LmResp{V<:FPVector} <: ModResp  # response in a linear model
     mu::V                                  # mean response
     offset::V                              # offset added to linear predictor (may have length 0)
     wts::V                                 # prior weights (may have length 0)
+    fweights::FrequencyWeights
+    pweights::ProbabilityWeights
+    aweights::AnalyticWeights
     y::V                                   # response
-    function (::Type{LmResp{V}}){V}(mu::V, off::V, wts::V, y::V)
+    function (::Type{LmResp{V}}){V}(mu::V, off::V, wts::V, y::V,
+              fweights::FrequencyWeights, pweights::ProbabilityWeights,
+              aweights::AnalyticWeights)
         n = length(y)
         length(mu) == n || error("mismatched lengths of mu and y")
         ll = length(off)
         ll == 0 || ll == n || error("length of offset is $ll, must be $n or 0")
         ll = length(wts)
         ll == 0 || ll == n || error("length of wts is $ll, must be $n or 0")
-        new{V}(mu, off, wts, y)
+        ll = length(fweights)
+        ll == 0 || ll == n || error("length of fweights is $ll, must be $n or 0")
+        ll = length(pweights)
+        ll == 0 || ll == n || error("length of pweights is $ll, must be $n or 0")
+        ll = length(aweights)
+        ll == 0 || ll == n || error("length of aweights is $ll, must be $n or 0")
+        new{V}(mu, off, normalizeWeights(fweights, pweights, aweights), y,
+               fweights, pweights, aweights)
     end
 end
 convert{V<:FPVector}(::Type{LmResp{V}}, y::V) =
-    LmResp{V}(zeros(y), similar(y, 0), similar(y, 0), y)
+    LmResp{V}(zeros(y), similar(y, 0), similar(y, 0), y,
+    similar(y, 0), similar(y, 0), similar(y, 0))
 
 function convert{T<:Real}(::Type{LmResp}, y::AbstractVector{T})
     yy = float(y)