After your code below you could add...
x1g <- X*(1-Group)
x2g <- X*Group
## different moderator effects in the different groups...
mgu <- gam(Y~s(M,by=x1g)+s(M,by=x2g))
## same moderator effects, but different intercepts in the groups...
mgc <- gam(Y~s(M,by=X)+Group)
## Alternative way of specifying different moderator effects
## in the 2 groups...
mgu.alt <- gam(Y~s(M,by=X)+s(M,by=x2g))
## compare alternatives via generalized AIC...
AIC(mgu,mgc,mgu.alt)
## ... in general you could use an approximate GLRT as well.
## e.g. anova(mgc, mgu), but in this case it makes no sense
## as mgc is estimated to have *higher* edf than mgu, clearly
## indicating that mgu is preferable by any sensible measure.
Note that mgu and mgu.alt will not generally give identical fitted
values, despite superficially appearing to be the same model
parameterized in two different ways. This is because the smoothing
penalties mean that the models are actually different (In mgu the
effects in the 2 groups are assumed to be smooth, but in mgu.alt the
difference in effects is assumed smooth, as is the overall effect...)
best,
Simon
On 27/06/11 13:53, Denes Toth wrote:>
>
> Dear UseRs,
>
>
> I built varying coefficient models (in mgcv) for two groups
> separately, with one explanatory and one moderator variable (see the
> example below).
>
> # ------- # Example:
> # ------ # generate moderator variable (can the> same for both groups)
> modvar<- c(1:1000)>
> # generate group1 values
> x1<- rnorm(1000)
> y1<- scale(cbind(1,poly(modvar,2))%*%c(1,2,1)*x1 +
rnorm(1000,0,0.3))>
> # generate group2 values
> x2<- rnorm(1000)
> y2<-scale(cbind(1,poly(modvar,2))%*%c(-1,0.5,-1)*x2 +
rnorm(1000,0,0.3))>
> # separate models for each group
> mg1<- gam(y1~s(modvar,by=x1))
> mg2<-gam(y2~s(modvar,by=x2))>
>
> Having done this, the next step would be to test if there is a
> significant difference between the modvar-dependency of the
> coefficient of the X values (i.e. to check if the coefficients vary
> in a different manner in the two groups):
>
> # unified model Y<- c(y1,y2)
> X<- c(x1,x2)
> M<- c(modvar,modvar)> Group<- rep(c(0,1),each=1000) ???
>
>
> And at this point I do not know how to proceed. Maybe a permutation
> approach would be helpful, but I would be more comfortable with a
> built-in solution. I tried to find similar threads but had no
> success. Any help is highly appreciated.
>
> Regards, Denes [[alternative HTML version deleted]]
>
>
>
>
> ______________________________________________ R-help at r-project.org
> mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do
> read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
--
Simon Wood, Mathematical Science, University of Bath BA2 7AY UK
+44 (0)1225 386603 http://people.bath.ac.uk/sw283