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@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 ______________________________________________ R-help@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.
