[R] compare gam fits

Mike Lawrence Mike.Lawrence at dal.ca
Thu Aug 5 18:49:38 CEST 2010


Hi folks,

I originally tried R-SIG-Mixed-Models for this one
(https://stat.ethz.ch/pipermail/r-sig-mixed-models/2010q3/004170.html),
but I think that the final steps to a solution aren't mixed-model
specific, so I thought I'd ask my final questions here.

I used gamm4 to fit a generalized additive mixed model to data from a
AxBxC design, where A is a random effect (human participants in an
experiment), B is a 2-level factor predictor variable, and C is a
continuous variable that is likely non-linear. I tell gamm4 to fit a
smooth across C to each level of B independently, and I can use
predict.gam(...,se.fit=T) to obtain predictions from the fitted model
as well as the standard error for the prediction. I'd like to
visualize the BxC interaction to see if smoothing C within each level
of B was really necessary, and if so, where it is (along the C
dimension) that B affects the smooth. It's easy enough to obtain the
predicted B1-B2 difference function, but I'm stuck on how to convey
the uncertainty of this function (e.g. computing the confidence
interval of the difference at each value of C).

One thought is that predict.gam(...,se.fit=T) returns SE values, so if
I could find out the N on which these SE values are computed, I could
compute the difference CI as
sqrt( ( (SD_B1)^2 + (SD_B2)^2 ) / N ) * qt( .975, df=N-1 )

However, I can't seem to figure out what value of N was used to
compute the SEs that predict.gam(...,se.fit=T) produces. Can anyone
point me to where I might find N?

Further, is N-1 the proper df for the call to qt()?

Finally, with a smooth function and 95% confidence intervals computed
at each of a large number of points, don't I run into a problem of an
inflated Type I error rate? Or does the fact that each point is not
independent from those next to it make this an inappropriate concern?

Cheers,

Mike

-- 
Mike Lawrence
Graduate Student
Department of Psychology
Dalhousie University

Looking to arrange a meeting? Check my public calendar:
http://tr.im/mikes_public_calendar

~ Certainty is folly... I think. ~



More information about the R-help mailing list