Hello R users,
I have a puzzle with the VGAM package, on my first excursion into
generalized additive models, in that this very nice package seems to
want to do either more or less than what I want.
Precisely, I have a 4-component outcome, y, and am fitting multinomial
logistic regression with one predictor x. What I would like to find
out is, is there a single nonlinear function f(x) which acts in place
of the linear predictor x. There is a mechanistic reason to believe
this is sensible. So I'd like to fit a model
\eta_j = \beta_{ (j) 0 } + \beta_{ (j) x } f(x)
where both the function f(x) and its scaling coefficients \beta_{ (j)
x } are fit simultaneously. Here \eta_j is the linear predictor, the
logodds of outcome j vs the reference outcome. I cannot see how to fit
exactly this. Instead I seem to be able to do the following:
vgam(formula = y ~ s(x), family = multinomial)
fits the model
\eta_j = \beta_{ (j) 0 } + \beta_{ (j) x } f_j (x)
i.e. a different function f_j (x) is fit for each outcome.
vgam(formula = y ~ s(x), family = multinomial, constraints list(`(Intercept)`=
diag(1,3), 's(x)' = matrix(c(1,1,1),3,1)) )
fits the model
\eta_j = \beta_{ (j) 0 } + f (x)
i.e. a single function f (x) is fit, but scaled the same for each
outcome. I'd like one function, scaled differently.
Of course,
vgam(formula = y ~ x, family = multinomial)
fits the model
\eta_j = \beta_{ (j) 0 } + \beta_{ (j) x } x
which has the scaling, but not the nonlinear function.
Perhaps this is achievable using bs(), xij, and vglm, or even via the
constraint matrix, but I did not succeed.
Any help appreciated!
Edward
--
Edward Wallace, PhD
Postdoctoral fellow, Drummond lab
Harvard FAS center for Systems Biology
ewallace at cgr.harvard.edu
773-517-4009
Edward Wallace <ewjwallace <at> gmail.com> writes:> > Hello R users, > I have a puzzle with the VGAM package, on my first excursion into > generalized additive models, in that this very nice package seems to > want to do either more or less than what I want. > > Precisely, I have a 4-component outcome, y, and am fitting multinomial > logistic regression with one predictor x. What I would like to find > out is, is there a single nonlinear function f(x) which acts in place > of the linear predictor x. There is a mechanistic reason to believe > this is sensible. So I'd like to fit a model > \eta_j = \beta_{ (j) 0 } + \beta_{ (j) x } f(x) > where both the function f(x) and its scaling coefficients \beta_{ (j) > x } are fit simultaneously. Here \eta_j is the linear predictor, the > logodds of outcome j vs the reference outcome. I cannot see how to fit > exactly this. Instead I seem to be able to do the following: >Hello, try rrvglm(y ~ 1 + bs(x), fam = multinomial, trace = TRUE) It seems what you want is a stereotype model with a smooth function. Unfortunately rrvglm() is restricted to regression splines. You could extract out the scaling coeffs and feed them into vgam() using the constraints argument, but that would not be optimal in any strict sense. cheers Thomas
Edward Wallace <ewjwallace <at> gmail.com> writes:> > Hello R users, > I have a puzzle with the VGAM package, on my first excursion into > generalized additive models, in that this very nice package seems to > want to do either more or less than what I want. > > Precisely, I have a 4-component outcome, y, and am fitting multinomial > logistic regression with one predictor x. What I would like to find > out is, is there a single nonlinear function f(x) which acts in place > of the linear predictor x. There is a mechanistic reason to believe > this is sensible. So I'd like to fit a model > \eta_j = \beta_{ (j) 0 } + \beta_{ (j) x } f(x) > where both the function f(x) and its scaling coefficients \beta_{ (j) > x } are fit simultaneously. Here \eta_j is the linear predictor, the > logodds of outcome j vs the reference outcome. I cannot see how to fit > exactly this.Thomas Yee wrote :> Hello, > > try > > rrvglm(y ~ 1 + bs(x), fam = multinomial, trace = TRUE) > > It seems what you want is a stereotype model with > a smooth function. > Unfortunately rrvglm() is restricted to regression splines.Thank you very much. This seems to work, but occasionally quits producing the cryptic error "Error in devmu[smallmu] = smy * log(smu) : NAs are not allowed in subscripted assignments" Any ideas?> You could extract out the scaling coeffs and feed them > into vgam() using the constraints argument, but that > would not be optimal in any strict sense.Really? I thought I had tried adding constraint matrices such as [1 ; 2] and vglm raised an error saying it needed entries to be 1 or 0. I can check that if you'd like. Edward -- Edward Wallace, PhD Harvard FAS center for Systems Biology +1-773-517-4009