GLAM_2d_covariates function

Fit the 2d GLAM with covariates