[R] mgcv: Plotting probabilities for binomial GAM with crossed random intercepts and factor by variable

Jan Vanhove jan.vanhove at unifr.ch
Thu Jan 10 15:26:15 CET 2013

mgcv: Constructing probabilities for binomial GAM with crossed random 
intercepts and factor by variable


(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 

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 

xtabs(~ Modality + Correct, dat)
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?)


More information about the R-help mailing list