1

I'm wondering if there is a way to plot a difference curve for two smooths involving the intercept.

Here is a summary of the model smooths.

s(percent)                 
  s(percent):x1x2voiced.D    
  s(percent):x1x2voiceless.f 
  s(percent):x1x2voiced.G    
  s(percent):x1x2voiceless.h 
  s(percent):x1x2voiceless.s 
  s(percent):x1x2voiceless.S 
  s(percent):x1x2voiceless.T 
  s(percent):x1x2voiceless.x 
  s(percent):x1x2voiceless.X-
  s(percent):x1x2voiced.z

The following code (from itsadug) works fine.

plot_diff(m1, view='percent', comp=list(x1x2=c('voiced.D', 'voiceless.f')))

But I'm interested also in other contrasts such as s(percent) vs. s(percent):x1x2voiced.G. Is there a way to include the first smooth s(percent) in the comp call?

I'm aware of this code (from mgcViz) which works fine:

plotDiff(s1 = sm(b, 1), s2 = sm(b, 2)) + l_ciPoly() + l_fitLine() + 
    geom_hline(yintercept = 0, linetype = 2)

But I'd love something following the same approach as itsadug in order to integrate the output with other functions such as get_smooths_difference() from tidymv.

kjetil b halvorsen
  • 77,844
Dallak
  • 21
  • What does your model structure look like? - show the gam() call please. Are these using the ordered factor trick to get a reference smooth plus smooth differences? – Gavin Simpson Sep 14 '22 at 09:38
  • Yes, I am using ordered factors. Here it is.

    bam(val ~ x1x2 + s(percent, bs= "cr", k=10) + s(percent, bs= "cr", k=10, by=x1x2) + s(percent, stim, bs="fs", m=1, k=10, by=x1x2)+ s(percent, stim, bs="fs", m=1, k=10)+ s(percent, word, bs="fs", m=1, k=10), data = kurt_gam_aa_ini, AR.start=kurt_gam_aa_ini$start.event, method = "fREML", discrete=T, family = "scat", nthread=8)))

     – Dallak Sep 14 '22 at 13:36
  • s(percent) refers to the reference level of x1x2, so can't you just include that in the comp definition? Something like comp = list(x1x2, = c(levels(kurt_gam_aa_ini$x1x2)[1L], "voiced.G")). The way the differencing works is that it generates data over percent for the 2 levels of the factor mentioned, all other covariates are held constant at typical values, then it predict()s using type = "lpmatrix" to get matrix Xp and then we take the rows of Xp that correspond to the first level, and rows of Xp that correspond to the other level we are comparing with, then... – Gavin Simpson Sep 14 '22 at 14:10
  • it subtracts these two sub matrices from one another to yield the difference Xp matrix, we then multiply %*% by the model coefficients to get the actual difference between the two smooths, and we compute a standard error for the difference using the vcov() of the fitted model. So the trick is just I think to know how to specify the data that you want, and in this case it's which levels of x1x2 you want to compare. – Gavin Simpson Sep 14 '22 at 14:13
  • I tried this code plot_diff(m1, view='percent', comp=list(x1x2=c(levels(m1$x1x2)[1L], "voiced.G"))) , but it produced this error Error in get_difference(model, comp = comp, cond = cond, se = ifelse(se > : Provide two levels for x1x2 to calculate difference. – Dallak Sep 14 '22 at 14:47
  • I guess the issue is just in specifying the data I want to use as your noted, which is something I'm not familiar with. – Dallak Sep 14 '22 at 14:59

1 Answers1

1

I guess I figured it out following Gavin's comment. Here is the answer for anybody else who might come later looking for the answer.

This code extracts the reference level for the ordered factors.

plot_diff(cog_ini__aa_m2, view='percent', comp=list(x1x2=c(cog_ini__aa_m2[["xlevels"]]$x1x2[1], "voiceless.X-"))) enter image description here

And here is how I fed it to tidymv for extra customizations.

get_smooths_difference(cog_ini__aa_m2, percent, list(x1x2=c(cog_ini__aa_m2[["xlevels"]]$x1x2[1], "voiceless.X-")))-> inter_diff

inter_diff %>% ggplot(aes(percent, difference, group = group)) + geom_hline(aes(yintercept = 0), colour = "#8f5f3f") + geom_ribbon(aes(ymin = CI_lower, ymax = CI_upper, fill = sig_diff), alpha = 0.3) + geom_line(aes(colour = sig_diff), size = 1) + scale_colour_manual(values = c("#e35760", "#6f849c")) + scale_fill_manual(values = c("#e35760", "#6f849c")) + labs(colour = "significant", fill = "significant") + theme(legend.position = "top")

This produces the plot I'm looking for. enter image description here

Thanks Gavin for your time and help!

Dallak
  • 21