From 75bdaf0ca24901960c26c030c318fe9909f9fd29 Mon Sep 17 00:00:00 2001 From: Joshua Krause <52180639+JoKra1@users.noreply.github.com> Date: Tue, 16 Jan 2024 10:04:10 +0100 Subject: [PATCH] Changes to Difference constraints --- src/mssm/src/python/constraints.py | 3 ++- src/mssm/src/python/formula.py | 12 +++++++----- src/mssm/src/python/penalties.py | 7 +++++-- 3 files changed, 14 insertions(+), 8 deletions(-) diff --git a/src/mssm/src/python/constraints.py b/src/mssm/src/python/constraints.py index 452b509..0226761 100644 --- a/src/mssm/src/python/constraints.py +++ b/src/mssm/src/python/constraints.py @@ -38,7 +38,8 @@ class Constraint: Line 3 shows how the constraint can be absorbed for model fitting by first computing new bases (f_i(x) - f_{i-1}(x)) and then estimating the coefficients based on those (this is implemented in mgcv's ``smoothCon`` function when setting ``sparse.cons=2``). Note that the constrained model, just like when using the QR-based or dropping approach, requires dropping one of the k original coefficients for estimation since only k-1 - coefficients are allowed to vary. + coefficients are allowed to vary. The same sum-to-zero constraint can be achieved by fixing one of the central bases in the original model to it's neighbor (e.g., setting b_2 = -b_3) or + by setting b_1= -b_2 - b_3 - b_4 - b_5. mssm fixes one of the central bases to it's neighbor. With a B-splines basis, it would be necessary that b_1 = b_2 = b_3 = b_4 = b_5 for the model to fit a constant f(x) over all values of x (Eilers & Marx, 1996). In the constrained model this is no longer possible due to b_1 = -b_2, effectively removing the constant from the space of functions that f(x) can approximate. diff --git a/src/mssm/src/python/formula.py b/src/mssm/src/python/formula.py index b68a3ac..0040756 100644 --- a/src/mssm/src/python/formula.py +++ b/src/mssm/src/python/formula.py @@ -773,7 +773,7 @@ def __absorb_constraints(self): elif term_constraint == ConstType.DROP: sterm.Z.append(Constraint(int(matrix_term_v.shape[1]/2),ConstType.DROP)) elif term_constraint == ConstType.DIFF: - sterm.Z.append(Constraint(None,ConstType.DIFF)) + sterm.Z.append(Constraint(int(matrix_term_v.shape[1]/2),ConstType.DIFF)) else: if vi == 0: matrix_term = matrix_term_v @@ -789,7 +789,7 @@ def __absorb_constraints(self): elif term_constraint == ConstType.DROP: sterm.Z.append(Constraint(int(matrix_term.shape[1]/2),ConstType.DROP)) elif term_constraint == ConstType.DIFF: - sterm.Z.append(Constraint(None,ConstType.DIFF)) + sterm.Z.append(Constraint(int(matrix_term.shape[1]/2),ConstType.DIFF)) def build_penalties(self): """Builds the penalties required by ``self.__terms``. Called automatically whenever needed. Call manually only for testing.""" @@ -1401,8 +1401,9 @@ def build_smooth_term_matrix(ci,n_j,has_scale_split,sti,sterm,var_map,var_mins,v # Applies difference re-coding for sum-to-zero coefficients. # Based on smoothCon in mgcv(2017). See constraints.py # for more details. - matrix_term_v = np.diff(matrix_term_v) - + matrix_term_v = np.diff(np.concatenate((matrix_term_v[:,sterm.Z[vi].Z:matrix_term_v.shape[1]],matrix_term_v[:,:sterm.Z[vi].Z]),axis=1)) + matrix_term_v = np.concatenate((matrix_term_v[:,matrix_term_v.shape[1]-sterm.Z[vi].Z:],matrix_term_v[:,:matrix_term_v.shape[1]-sterm.Z[vi].Z]),axis=1) + if vi == 0: matrix_term = matrix_term_v else: @@ -1414,7 +1415,8 @@ def build_smooth_term_matrix(ci,n_j,has_scale_split,sti,sterm,var_map,var_mins,v elif sterm.Z[0].type == ConstType.DROP: matrix_term = np.delete(matrix_term,sterm.Z[0].Z,axis=1) elif sterm.Z[0].type == ConstType.DIFF: - matrix_term = np.diff(matrix_term) + matrix_term = np.diff(np.concatenate((matrix_term[:,sterm.Z[0].Z:matrix_term.shape[1]],matrix_term[:,:sterm.Z[0].Z]),axis=1)) + matrix_term = np.concatenate((matrix_term[:,matrix_term.shape[1]-sterm.Z[0].Z:],matrix_term[:,:matrix_term.shape[1]-sterm.Z[0].Z]),axis=1) m_rows, m_cols = matrix_term.shape #print(m_cols) diff --git a/src/mssm/src/python/penalties.py b/src/mssm/src/python/penalties.py index fec1137..87a8e1e 100644 --- a/src/mssm/src/python/penalties.py +++ b/src/mssm/src/python/penalties.py @@ -57,7 +57,8 @@ def diff_pen(n,constraint,m=2): S = np.delete(np.delete(S,Z,axis=1),Z,axis=0) D = np.delete(D,Z,axis=0) elif constraint.type == ConstType.DIFF: - D = np.diff(D,axis=0) # Correct for column differencing applied to X! See smoothCon help for mgcv (Wood, 2017) + D = np.diff(np.concatenate((D[Z:D.shape[0],:],D[:Z,:]),axis=0),axis=0) # Correct for column differencing applied to X! See smoothCon help for mgcv (Wood, 2017) + D = np.concatenate((D[D.shape[0]-Z:,:],D[:D.shape[0]-Z,:]),axis=0) S = D @ D.T S = scp.sparse.csc_array(S) @@ -113,7 +114,9 @@ def TP_pen(S_j,D_j,j,ks,constraint): S_TP = scp.sparse.csc_array(np.delete(np.delete(S_TP.toarray(),Z,axis=1),Z,axis=0)) D_TP = scp.sparse.csc_array(np.delete(D_TP.toarray(),Z,axis=0)) elif constraint.type == ConstType.DIFF: - D_TP = np.diff(D_TP.toarray(),axis=0) # Correct for column differencing applied to X! See smoothCon help for mgcv (Wood, 2017) + D_TP = D_TP.toarray() + D_TP = np.diff(np.concatenate((D_TP[Z:D_TP.shape[0],:],D_TP[:Z,:]),axis=0),axis=0) # Correct for column differencing applied to X! See smoothCon help for mgcv (Wood, 2017) + D_TP = np.concatenate((D_TP[D_TP.shape[0]-Z:,:],D_TP[:D_TP.shape[0]-Z,:]),axis=0) S_TP = D_TP @ D_TP.T S_TP = scp.sparse.csc_array(S_TP) D_TP = scp.sparse.csc_array(D_TP)