I aim to plot the gam() model fit back on the original data. My impression was that gam() has special treatment for ordered factors, such that a smooth is calculated for the reference level, and then difference smooths are calculated with respect to the reference smooth. I therefore expect that the parametric components will correspond to the mean of the reference level and to the mean difference between the reference level and the smooth.

In testing this behaviour, I found that if the factor is not ordered, this expectation holds. For an ordered factor, the global intercept instead goes through the mean of all the data, so does not pertain to the reference level directly. For toy data with global and difference-intercepts of +5 and -2, respectively, the unordered factor returns reasonable estimates of 4.75 and -2.00, whereas the ordered factor estimates are 3.76 and -1.41. I do not know what these intercepts signify. The resulting plots show the smooths are correct, but the parametric components miss the data entirely. What am I missing in the ordered factor parametrisation that needs to be accounted for?example by-smooth



glob.int <- 5
diff.int <- -2

x <- seq(0,10, by=0.25)
y1 <- glob.int-1/(x+1)
y2 <- y1 + rnorm(length(x), 0, 0.05) + diff.int
y1 <- y1 + rnorm(length(x), 0, 0.05)

foo <- data.frame( x=rep(x,2), 
                   fac1_=factor(rep( c('0','1'), each=length(x))),
                   fac2_=factor(rep( c('0','1'), each=length(x)),ordered=TRUE))

Model1 <- gam( y~s(x, bs='ad', k=length(x)) +
                 s(x,by=fac1_, m=1) +
Model2 <- gam( y~s(x, bs='ad', k=length(x)) +
                 s(x,by=fac2_, m=1) +
               method='ML' )


smdat1 <- plot(Model1, all.terms=FALSE, seWithMean=TRUE, pages=1)
smdat2 <- plot(Model2, all.terms=FALSE, seWithMean=TRUE, pages=1)

ref_fit1 <- data.frame( x=smdat1[[1]]$x, fit=smdat1[[1]]$fit )
diff_fit1 <- data.frame( x=smdat1[[2]]$x, fit=smdat1[[2]]$fit )
int.glob1 <- coefficients(Model1)[1]
int.fac1 <- coefficients(Model1)[2]

ref_fit2 <- data.frame( x=smdat2[[1]]$x, fit=smdat2[[1]]$fit )
diff_fit2 <- data.frame( x=smdat2[[2]]$x, fit=smdat2[[2]]$fit )
int.glob2 <- coefficients(Model2)[1]
int.fac2 <- coefficients(Model2)[2]

ref_fit1$fit_trans <- (ref_fit1$fit+int.glob1)
diff_fit1$fit_trans <- (ref_fit1$fit+diff_fit1$fit+int.glob1+int.fac1)
ref_fit2$fit_trans <- (ref_fit2$fit+int.glob2)
diff_fit2$fit_trans <- (ref_fit2$fit+diff_fit2$fit+int.glob2+int.fac2)

plot(y~x, data=foo, col=fac1_, pch=16, ylim=c(0,8))
lines(fit_trans~x, data=ref_fit1, col='black')
lines(fit_trans~x, data=diff_fit1, col='red')
lines(fit_trans~x, data=ref_fit2, col='black', lty='dashed')
lines(fit_trans~x, data=diff_fit2, col='red', lty='dashed')

legend( x = 0, y=8, c('unordered reference','unordered difference','ordered reference','ordered difference'), 
        lty=rep(c('solid','dashed'),each=2), col=rep(c('black','red'),2) )

The ordered factor model uses polynomial contrasts instead of the usual (default) treatment contrasts for unordered factors. As such, the groups means are not parameterised in terms of the mean for a reference group and deviations of the other group means from the reference.

In the unordered factor, the intercept is the estimated mean of the reference level:

> with(foo, mean(y[fac1_ == 0]))
[1] 4.757463
> coef(Model1)[1]

and the second coefficient is the difference between the mean for level 1 and the reference level

> with(foo, mean(y[fac1_ == 1]) - mean(y[fac1_ == 0]))
[1] -2.000921
> coef(Model1)[2]

In the ordered factor representation, polynomial contrasts are used, which do not result in binary (0,1) indicators for groups. For the unordered factor the contrasts are:

> contr.treatment(2)
1 0
2 1

So the model matrix for this set-up contains a vector of 1s for the intercept and vector of 0s or 1a indexing if the $i$th observation is from level 0 (the reference) or level 1 of fac1_

> model.matrix(Model1)[c(1:3, nrow(foo)- (0:2)), 1:2]
   (Intercept) fac1_1
1            1      0
2            1      0
3            1      0
82           1      1
81           1      1
80           1      1

The first three rows here are from observations from the reference level (fac1_ = 0) and the last three rows are for observations from the other level. When you are predicting, the coefficients are multiplied by these columns.

With polynomial contrasts, the "data" representing different levels, the parameterisation, is different. Now the intercept represents the overall mean of the response

> coef(Model2)[1]
> with(foo, mean(y))
[1] 3.757003

The second coefficient is likewise not the difference between the overall mean and the second group.

The second coefficient in Model2 need to be multiplied by the values of the orthonormal polynomial contrasts

> contr.poly(2)
[1,] -0.7071068
[2,]  0.7071068

and then we add the intercept to recover the group means

> coef(Model2)[1] + (contr.poly(2) * coef(Model2)[2])
[1,] 4.757463
[2,] 2.756543
> aggregate(y ~ fac1_, data = foo, FUN = mean)
  fac1_        y
1     0 4.757463
2     1 2.756543

These are the values you need to add to data from plot.gam() to recover the fits on the scale of the link function (which here is also the scale of the response).

  • 1
    $\begingroup$ Thank you, I assumed the issue was a result of my not understanding something fundamental about the contrasts. $\endgroup$
    – noname
    Jun 4 '19 at 21:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.