Jan Vanhove
2013-Jan-10 14:26 UTC
[R] mgcv: Plotting probabilities for binomial GAM with crossed random intercepts and factor by variable
mgcv: Constructing probabilities for binomial GAM with crossed random intercepts and factor by variable Hello, (I'm sorry if this has been discussed elsewhere; I may not have been looking in the right places.) I ran a binomial GAM in which "Correct" is modelled in terms of the participant's age and the modality in which the stimulus is presented (written vs spoken). Participants ("Subject") and stimuli are also included as crossed random intercepts. age.gam <- bam(Correct ~ Modality + s(Age, by=Modality) + s(Subject, bs="re") + s(Stimulus, bs="re"), data = dat, family="binomial") Parametric coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.242 1.121 -2.000 0.0455 * ModalityWritten -1.043 1.263 -0.826 0.4087 Approximate significance of smooth terms: edf Ref.df Chi.sq p-value s(Age):ModalitySpoken 4.883 5.477 98.45 <2e-16 *** s(Age):ModalityWritten 5.675 6.363 126.81 <2e-16 *** s(Subject) 133.296 161.000 1472.35 <2e-16 *** s(Stimulus) 83.376 85.000 5100.59 <2e-16 *** I would now like to plot the predicted probabilities for "Correct" == "yes", say for written stimuli: plot(age.gam, select=2, shift=(-2.242-1.043), trans=plogis) Though the overall shape of the curve is about what I expected based on a earlier visual exploration, the plotted probabilities are much too low (about 5% for values of Age where I'd expected 40%). Here are the raw counts: xtabs(~ Modality + Correct, dat) Correct Modality 0 1 Spoken 4146 2693 Written 4323 3011 Is it possible that I need to substract the values of s(Subject).1 and s(Stimulus).1 from coef(age.gam) to the "trans" argument, too? Those are 0.55 and -2.86, respectively, which would provide a much better match between the plotted and the expected probabilities. (But then what does the "(Intercept)" represent in this model?) Thanks, Jan