Question for factor levels in gam

[Nataliannash]

Dear Gavin,

I used gam to simulate and then predict Y by X, Z, and Xt variables, but I am not sure about how to use these variables in gam. I appreciate it if you could help me about it.

X and Z are numbers and I assume there is an interaction between them.
Xt is treatment as factor (with 2 levels, 1 = low and 2 = high).
Xt has a direct effect on X variable. The increase in Xt should decreases the X value.
I used the gam function as follow:

gm1 = gam(Y ~ s(X) + s(X, by = Xt) + ti(X, Z), data = Input_data , family = Gamma(link ="log") ,
                   method="REML" ,na.action = na.exclude)  
Predict_data$Y1 = predict(gm1, exclude = ‘s(X, by = Xt)’, new_data , type = "response")

1. If I want to find the effects of treatment Xt on X and then Y variable which form of these smooth functions are correct?

   1. s(X, by = Xt)
   2. s(X) + s(X, by = Xt)
   3. s(X , Xt, bs = 'tp' )
   4. s(X , Xt, bs = c( 'tp', ‘re’) )

2. Is it right if I just use s(X) to measure the effects of Xt and X on Y? because the effect of Xt treatment on Y is already reflected on X variable.

3. Can I use the predict function by excluding s(X, by = Xt) for a new dataset without Xt variable? (the new dataset contains dummy data of all variables except Xt)

Thank you
Nat

[gavinsimpson]

Q1

Models iii and iv are wrong - you can't smooth a factor and with the odd exception (the bs = "fs" being one) you can't specify variables in different units in s() (iii) and you can't specify separate marginal bases with s() (iv).

Models i and ii could work but you must include Xt as a parametric effect also:

• Xt + s(X, by = Xt)
• Xt + s(X) + s(X, by = Xt)

which one you want depends on several things. The first would model two entirely separate smooths, one per level of Xt. The second would model some average effect of X plus separate smooth effects of X per level of Xt; this model assumes there is some common relationship (that there is some common pattern to the effect of X on Y), but as you have ti currently there is a potential identifiability issue as there are now multiple smooths of X in the model. To address the identifiability problem, you could add m = 1 to the last smooth (s(X, by = Xt, m = 1)) which has the effect of penalising those smooths in terms of deviations from the constant function (which implies deviations from the common or global smooth) or you could recode Xt as an ordered factor, in which case s(X) would estimate the smooth effect of X on Y for the reference level of Xt, while s(X, by = Xt) would model the smooth difference between the effect of the reference level and the other level of Xt.

Q2

No, I don't think this is right. You have a potential (or known) interaction between the effect of X and the factor Xt; the effect of X on Y depends on Xt. Hence i would model this as a factor-by as described above.

Q3

If you fitted the model as any of these

• Xt + s(X, by = Xt)
• Xt + s(X) + s(X, by = Xt, m = 1)
• oXt + s(X) + s(X, by = oXt)

then no, you can't exclude the factor by term - well you could but it wouldn't you because you have to include the group means via the parametric Xt term. Note that I have indicated the ordered version of the model by inclusion of oXt in the formula in place of Xt.

What I would do instead would be to generate predictions under the two scenarios of Xt == "low" and Xt == "high" for the dummy values of X and Z.

You could perhaps treat Xt as a random effect but the justification for this is limited if Xt is, as you say, some fixed treatment effect.