hendersontrent / gam.jl Goto Github PK

Could be due to the fact it is an integer variable.

This is looking awesome. Super well done. I am just putting this out there to keep track of the idea. On page 177 of his text, Wood describes how how to build the penalty matrix for a GAM with more than one term. Described in R terms, (whats comfortable for me), basically it consists of a cbind of the basis matricies and then rbinding a block diagonal matrix of the penalized difference matrices. I provided a code example below. This will allow for a single optimization run for a vector of penalty terms. Thankfully due to Julia's broadcasting capabilities, I think there should be minimal refactoring necessary. I think it might also help with the fit, but I am not 100% certain if that is the case. I am mor than happy to look into building this out, but I know development has been fast as furious so don't want to get in the way and break things. Below is code example assuming the trees dataset you are using in your tests.

using SparseArrays

# response variable
y = trees.Volume

# create vector if input variables
x = [trees.Girth, trees.Height]

# broadcast basis functions
Basis = QuantileBasis.(x,10,4)
X = BasisMatrix.(Basis, x) # Basis Matrix
D = DifferenceMatrix.(Basis) # D penalty matrix

# vector of penalties
λ = [10, 10]

# penalty design matrix
X_p = Matrix(
    vcat(
        # cbind the Basis Matricies
        hcat(X...), 
        # create a block diagonal matrix of penalized differences
        blockdiag((sqrt.(λ).*sparse.(D))...)
    )
)

# augmented penalty response
y_p = vcat(y, repeat([0],sum(first.(size.(D)))))

# fit model
lm(X_p, y_p)

Again, awesome work!

Black background with white writing!

Currently, GradientDescent is used, but this could be flexibly specified.

Currently, the code is as per the below, but should we really be using RSS for say a Poisson problem?

function GCV(param::AbstractVector, Basis::BSplineBasis{Vector{Float64}}, x::AbstractVector, y::AbstractVector)
n = length(Basis.breakpoints)
Xp, yp = PenaltyMatrix(Basis, param[1], x, y)
β = coef(lm(Xp,yp))
H = Xpinv(Xp'Xp)Xp' # Hat matrix
trF = sum(diag(H)[1:n])
y_hat = Xpβ
rss = sum((yp-y_hat)[1:n].^2) # Residual sums of squares
gcv = n*rss/(n-trF)^2
return gcv
end

When a heavy package like StatsPlots is a required dependency, then GAM is more difficult to use since StatsPlots or it's many dependencies might conflict with other packages. It's probably better to make it an optional dependency or move the plotting functions to a separate package

Is there any location where documentation can be found for this package?

This should then make it a proper generalised additive model once the intercept and all the smooth terms can be added.

I am aware this issue might be out of scope with respect to this package, but I couldn't think of a better place to ask ☺️

Would you know by any chance if it's feasible to specify a GAM in Turing, and if so how would such a mode look like, assuming the simple case y ~ s(x) where y and x are two continuous variables.

Hello,

Sorry if this is the wrong format to reach out. I am interested in contributing to the project. My background: I have an MPH in epidemiology, and I have been working with excess mortality time series models for the past 1+ year (R, GLM and mgcv packages). I have a master's in epidemiology but am taking coursework to apply for biostatistics PhD programs. I've developed proprietary packages for my organization, but I would be new to contributing to modeling packages. Thanks for your work on this package, and I would like to contribute in any way that I can.

-Jon

@yahrMason my current approach on main doesn't make use of the mutable struct you made. Our next iteration should!