Implementing Generalized Linear Models in Julia

glm.jl

type Link
    name::String
    linkFun::Function
    linkInv::Function
    muEta::Function
end

type Dist
    name::String
    variance::Function
    devResid::Function
    deviance::Function
end

type GlmResp
    dist::Dist
    link::Link
    eta::Vector{Float64}
    mu::Vector{Float64}
    offset::Vector{Float64}
    wts::Vector{Float64}
    y::Vector{Float64}
end

GlmResp(dist::Dist, link::Link, y::Vector{Float64}) = (n=length(y);zz=zeros(Float64,(n,));GlmResp(dist,link,zz,zz,zz,ones(Float64,(n,)), y))

logitLink    = Link("logit", (mu)-> log(mu ./ (1 - mu)), (eta)-> 1 ./ (1 + exp(-eta)), (eta)-> (e = exp(-eta); f = 1. + e; e ./ (f .* f)))
logLink      = Link("log", (mu)-> log(mu), (eta)-> exp(eta), (eta)-> exp(-eta))
identityLink = Link("identity", (mu)-> mu, (eta)-> eta, (eta)-> ones(eltype(eta), size(eta)))
inverseLink  = Link("identity", (mu)-> 1 ./ mu, (eta)-> 1 ./ eta, (eta)-> -1. ./ (eta .* eta))

## some utility functions
function logN0(x)
    ret = [x]
    ret[ret > 0] = log(ret[ret > 0])
    ret
end
y_log_y(y, mu) = y .* logN0(y ./ mu)
 
GaussianDist =
    Dist("Gaussian",
         (mu)-> ones(typeof(mu), size(mu)),
         (y, mu, wt)-> (r = y - mu; wt .* r .* r),
         (y, mu, wt, drsum)-> (n = length(mu); n * (log(2*pi*drsum/n) + 1) + 2 - sum(log(wt))))

gammaDist =
    Dist("gamma",
         (mu)-> mu .* mu,
         (y, mu, wt)-> -2 * wt .* (logN0(y ./ mu) - (y - mu) ./ mu),
         (y, mu, wt, drsum)-> (n=sum(wt); disp=drsum/n; invdisp(1/disp); sum(wt .* dgamma(y, invdisp, mu * disp, true))))

BernoulliDist =
    Dist("Bernoulli",
         (mu)-> max(eps(Float64), mu .* (1. - mu)),
         (y, mu, wt)-> 2 * wt .* (y_log_y(y, mu) - y_log_y(1 - y, 1 - mu)),
         (y, mu, wt, drsum)-> -2. * sum(y * log(mu) + (1 - y) * log(1 - mu)))

PoissonDist =
    Dist("Poisson",
         (mu)-> mu,
         (y, mu, wt)-> 2 * wt .* (y .* logN0(y ./ mu) - (y - mu)),
         (y, mu, wt, drsum)-> -2 * sum(dpois(y, mu, true) * wt))

updateMu(r::GlmResp, linPr::Vector{Float64}) = (r.eta = linPr + r.offset; r.mu = r.link.linkInv(r.eta); nothing)