CODEWITHSUNDEEP
rdrawhierarchical-clusteringgammgcv
David Smith
David Smith

Reputation: 383

Using gratia::draw() in R to display partial effect plots within an HGAM that are not relative to the global smooth

I have a data set that looks like this:

df <- data.frame(
  Lake = structure(c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 
    2L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 
    1L, 2L, 2L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 
    1L, 2L, 1L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 1L, 1L, 
    1L, 1L, 2L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 2L, 2L, 
    2L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 
    2L, 2L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 2L, 2L, 1L, 2L, 1L, 2L, 1L, 
    1L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 
    1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 
    2L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 
    2L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 
    1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 2L, 
    2L, 2L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 1L, 
    1L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 
    2L, 2L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 1L, 
    2L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 
    1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 
    2L, 2L, 1L, 1L, 2L, 2L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 
    2L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 
    2L, 2L, 2L, 2L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 
    2L, 1L, 2L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 2L, 1L, 1L, 1L, 
    1L, 2L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 
    2L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 
    2L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 
    2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 
    2L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 
    1L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 
    1L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 
    1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 
    1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 2L, 
    2L, 2L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 
    1L, 1L, 1L, 1L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 
    2L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 1L, 2L, 1L, 
    2L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 
    2L, 2L, 2L, 1L, 2L, 2L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 
    1L, 2L, 1L, 2L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 
    1L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 
    2L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 1L, 1L, 
    1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 
    1L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 1L, 
    1L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 
    1L, 2L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 
    1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 1L, 2L, 
    2L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 
    1L, 2L, 1L, 1L, 2L), .Label = c("Fork", "TBend"), class = "factor"),
  Depth = c(1.63, 3, 10, 3, 10, 4, 13, 17, 14, 2.81, 20, 3, 28, 24, 6, 
    1.81999999999999, 7, 25, 2.20999999999998, 10, 15, 7.25999999999999, 
    4, 4, 6.64999999999998, 8.83999999999997, 6, 2.20999999999998, 
    22.96, 5.63, 11, 30, 32.31, 25, 1, 3, 4, 7.41000000000003, 2, 
    6, 17, 7, 5, 4.20999999999998, 3, 22, 5, 4.74000000000001, 7, 
    10, 3, 11, 14, 2, 24, 1, 7, 15, 16, 2.68000000000001, 12, 11, 
    5, 10, 10, 6, 12, 4, 4, 4.64999999999998, 18, 7.5, 13, 3, 15, 
    10, 22, 19, 4, 12, 2, 3, 5.41000000000003, 6, 19, 6, 3, 3, 34, 
    3.63, 11, 6, 7, 25, 4, 2.81, 4.70999999999998, 3, 12.31, 5, 17, 
    28, 3.63, 8, 9, 3, 30, 20, 11, 4, 12, 3, 4, 16, 5, 10, 2, 14, 
    58, 10, 2.06, 15, 2.74000000000001, 7, 10.74, 2.81, 11, 6, 5, 
    7.25999999999999, 10, 2.68000000000001, 9, 2.83999999999997, 
    5.5, 15, 7, 6.56, 14, 6, 3.25999999999999, 2.31, 1, 7, 3, 4, 
    2, 3, 9, 28, 18.84, 5, 5, 2.75999999999999, 7.63, 8.20999999999998, 
    18, 3, 11, 1, 24, 4, 22, 2, 3, 4.20999999999998, 14.65, 16, 9, 
    5, 3, 7, 1, 2, 4.5, 2, 20, 1, 10, 17, 4, 2, 1, 23, 5, 11, 12, 
    17, 10, 3, 18, 6, 7, 5, 3, 32, 16, 5, 7, 9, 29, 2, 12, 4, 23, 
    14, 4, 5, 11, 11.82, 6.20999999999998, 7, 12, 3, 6, 4, 17, 4, 
    24, 6, 12, 11.63, 4, 2, 25, 2, 54, 7, 8, 9.25999999999999, 14, 
    15, 11, 6, 21, 1, 3, 8, 1, 2.83999999999997, 19, 6, 19, 2.06, 
    3, 3, 4, 8, 6, 9.41000000000003, 4, 8.64999999999998, 3, 3, 2.5, 
    30, 12, 14, 15, 16, 10.56, 24, 12, 16.71, 25, 1, 10, 17, 1, 1.25999999999999, 
    12, 4, 24, 15, 8.68000000000001, 8, 3, 15.82, 17, 5, 3, 6.70999999999998, 
    5.63, 10, 10.68, 8, 3, 8.81, 5.25999999999999, 22, 12, 5.81999999999999, 
    6, 6, 3.5, 1.52999999999997, 4, 22, 15, 4, 23, 12, 25, 4, 22, 
    5.41000000000003, 9, 19, 8, 4, 8.56, 20, 10.21, 24, 1, 6, 3, 
    10, 3, 28, 12, 6, 17, 1, 3.41000000000003, 6.16000000000003, 
    4, 20.68, 4, 2.74000000000001, 5, 12, 1, 45, 4.74000000000001, 
    18, 15, 1, 8, 20, 21, 3, 16, 1, 3, 30, 10, 6.06, 4, 10.84, 25, 
    26, 12, 2.56, 2, 6, 10.56, 10.31, 16, 29.26, 5, 6, 3.81999999999999, 
    15, 1, 8, 3, 2, 22, 5, 2.95999999999998, 4.5, 1, 18, 2.66000000000003, 
    19, 12, 4, 14, 3, 7, 28, 4, 23, 6, 5, 3, 22, 1, 4, 12, 7, 1.63, 
    12.21, 15, 4, 3, 9, 20.65, 4.74000000000001, 22, 8.81, 5.81999999999999, 
    4.16000000000003, 7, 10, 24, 4.95999999999998, 30, 2, 10, 5, 
    9, 5, 12, 29.82, 2, 6.5, 6.20999999999998, 1, 1, 22, 22, 6.64999999999998, 
    32, 11, 15, 1, 18, 1.81999999999999, 4, 8, 20, 15, 4, 7, 22, 
    2, 2, 1, 1, 15, 20, 3, 5, 1.63, 4.66000000000003, 22, 6, 2, 31, 
    20, 5, 9.5, 30, 18, 13, 12, 12, 4.20999999999998, 12, 10.06, 
    2.68000000000001, 2, 1, 5, 2, 9, 2, 4, 1, 6, 1, 1, 2.16000000000003, 
    7, 8.95999999999998, 2.74000000000001, 5, 4, 5, 15, 20, 5.41000000000003, 
    29.41, 7, 32, 4, 14, 2.74000000000001, 4, 15, 8, 21, 32, 13.41, 
    3, 14, 4, 3, 18, 2.31, 25, 3.5, 4.74000000000001, 19, 21, 5.25999999999999, 
    10.21, 12.84, 2.95999999999998, 2, 4.31, 7, 7, 2.31, 17, 10.71, 
    23.41, 3, 3.41000000000003, 4.68000000000001, 22, 3, 13, 15, 
    8.74000000000001, 14.81, 5, 1, 4, 16, 1.41000000000003, 13, 3, 
    2, 6.06, 7, 3, 22, 4.83999999999997, 7, 2.81, 21, 3, 19, 6, 14, 
    2, 1, 10, 7.5, 8.70999999999998, 30, 14, 20, 1, 18, 30, 28, 1.41000000000003, 
    20, 5, 1.41000000000003, 3.5, 4.64999999999998, 5, 9.5, 3, 1.63, 
    11, 21, 2.66000000000003, 20.74, 15, 15, 14, 5, 14, 4.5, 4, 6.06, 
    4.20999999999998, 12, 18, 10.16, 7.81999999999999, 1, 2.95999999999998, 
    15, 2.5, 2.70999999999998, 11, 13.63, 18, 6, 18, 11, 6, 12, 7.5, 
    4.56, 1.38, 2.95999999999998, 17, 4, 1, 15, 4.74000000000001, 
    5.5, 11, 4, 1, 3, 25, 3, 9, 15, 11, 29, 8.56, 23, 14.65, 1, 7, 
    8, 14.06, 2, 3, 26, 2.56, 2.5, 25, 2.74000000000001, 1, 3, 8.56, 
    9.38, 2, 18, 3, 30, 16.96, 4, 22, 11, 6, 4, 3, 8.83999999999997, 
    22, 18, 7, 2.68000000000001, 6, 14.76, 7, 5, 1, 21, 3.81999999999999, 
    10, 3, 5, 7, 6, 20, 6.81, 7, 19, 24, 5, 1, 21.41, 3, 1.81999999999999, 
    10, 11.41, 6, 30, 3, 4, 4, 4, 1.5, 10.5, 18, 10, 2, 25, 14, 4, 
    5.63, 4.20999999999998, 2, 10.84, 10, 7, 30, 1, 17, 3, 3, 22, 
    2.74000000000001, 1, 8, 7, 32.65, 4, 3, 5, 4, 5, 1, 5, 10.76, 
    4, 2, 3.41000000000003, 4, 17),
  OrdDay = c(254, 271, 286, 88, 181, 209, 246, 259, 218, 324, 230, 181, 
    271, 351, 364, 224, 268, 232, 210, 215, 260, 281, 286, 351, 195, 
    167, 248, 54, 308, 254, 322, 125, 33, 248, 336, 319, 322, 238, 
    181, 304, 195, 181, 273, 210, 153, 230, 28, 349, 195, 78, 286, 
    41, 355, 109, 78, 187, 31, 286, 41, 336, 187, 146, 305, 70, 290, 
    129, 290, 160, 83, 195, 147, 7, 159, 195, 146, 195, 181, 11, 
    349, 230, 140, 146, 268, 305, 181, 244, 299, 124, 155, 254, 232, 
    218, 12, 78, 286, 324, 177, 131, 33, 304, 56, 211, 254, 218, 
    60, 167, 147, 167, 322, 181, 299, 167, 215, 351, 230, 334, 25, 
    63, 11, 246, 5, 281, 349, 209, 91, 324, 246, 63, 203, 281, 167, 
    336, 63, 167, 88, 248, 153, 184, 237, 28, 281, 33, 195, 167, 
    109, 260, 56, 268, 248, 259, 187, 11, 124, 75, 254, 54, 218, 
    319, 322, 91, 12, 204, 195, 211, 125, 54, 195, 271, 364, 83, 
    335, 75, 75, 109, 75, 299, 160, 124, 334, 7, 146, 153, 184, 129, 
    146, 181, 131, 364, 31, 124, 11, 304, 290, 181, 204, 195, 322, 
    290, 305, 28, 336, 101, 174, 335, 109, 322, 273, 304, 364, 224, 
    210, 246, 25, 305, 349, 319, 83, 160, 28, 224, 187, 254, 124, 
    7, 167, 195, 12, 12, 187, 281, 101, 336, 304, 195, 244, 75, 232, 
    322, 246, 167, 237, 167, 336, 5, 125, 232, 187, 204, 286, 268, 
    131, 195, 322, 155, 104, 325, 28, 215, 195, 224, 184, 224, 174, 
    177, 167, 21, 363, 244, 268, 281, 286, 286, 335, 286, 336, 286, 
    109, 224, 181, 322, 299, 177, 254, 124, 336, 268, 218, 324, 281, 
    12, 119, 224, 248, 187, 215, 234, 159, 7, 204, 167, 78, 167, 
    325, 244, 290, 238, 305, 322, 246, 334, 184, 195, 210, 335, 160, 
    248, 218, 299, 78, 322, 167, 41, 211, 184, 238, 21, 281, 336, 
    322, 349, 268, 363, 273, 334, 349, 83, 78, 75, 204, 25, 237, 
    104, 232, 195, 319, 363, 355, 5, 335, 167, 237, 349, 286, 184, 
    75, 91, 184, 33, 215, 281, 28, 78, 224, 215, 116, 268, 124, 248, 
    7, 70, 308, 160, 336, 237, 105, 195, 273, 305, 273, 155, 248, 
    281, 160, 209, 259, 63, 101, 143, 67, 187, 203, 11, 254, 210, 
    31, 167, 363, 70, 195, 91, 41, 324, 224, 21, 351, 146, 268, 308, 
    28, 334, 259, 56, 12, 232, 174, 224, 101, 335, 54, 195, 143, 
    25, 171, 195, 167, 336, 281, 203, 25, 224, 75, 218, 248, 160, 
    181, 237, 195, 133, 172, 146, 75, 143, 260, 215, 56, 254, 105, 
    271, 319, 88, 364, 12, 230, 271, 125, 203, 248, 211, 286, 54, 
    63, 5, 336, 259, 105, 28, 299, 224, 172, 125, 75, 299, 177, 105, 
    21, 28, 308, 91, 88, 63, 281, 167, 349, 238, 238, 204, 12, 237, 
    349, 91, 364, 174, 237, 63, 363, 268, 167, 28, 181, 155, 160, 
    33, 304, 244, 349, 248, 28, 281, 54, 167, 308, 116, 33, 224, 
    181, 33, 364, 177, 268, 268, 238, 336, 281, 181, 299, 246, 349, 
    324, 56, 75, 273, 271, 268, 195, 246, 181, 5, 248, 146, 322, 
    167, 140, 324, 286, 286, 174, 322, 60, 187, 260, 335, 104, 177, 
    167, 203, 304, 177, 232, 336, 209, 238, 125, 260, 268, 203, 195, 
    363, 88, 232, 254, 203, 246, 105, 349, 268, 160, 336, 336, 260, 
    88, 56, 5, 54, 363, 31, 21, 224, 260, 308, 355, 25, 177, 167, 
    254, 224, 70, 349, 281, 119, 7, 75, 184, 124, 308, 273, 146, 
    202, 167, 349, 88, 218, 70, 210, 160, 147, 155, 181, 244, 195, 
    56, 184, 41, 195, 160, 260, 101, 5, 116, 230, 351, 184, 25, 224, 
    349, 91, 67, 184, 124, 355, 237, 167, 209, 308, 167, 268, 31, 
    218, 101, 155, 167, 12, 125, 143, 336, 286, 75, 167, 187, 260, 
    304, 224, 203, 290, 125, 195, 290, 355, 324, 153, 187, 349, 355, 
    324, 238, 260, 224, 281, 238, 140, 290, 273, 119, 181, 153, 129, 
    271, 75, 230, 116, 41, 91, 167, 254, 54, 290, 167, 11, 237, 336, 
    105, 181, 11, 286, 244, 349, 91, 230, 336, 195, 119, 230, 349, 
    349, 203, 238, 63, 75, 335, 91, 268, 322, 83),
  stringsAsFactors = FALSE)

I am running an HGAM (GI model from Pederson et al. 2019 https://peerj.com/articles/6876/) that looks like this:

library(mgcv)
hgam_gi<-gam(Depth~Lake+s(OrdDay,bs="cc")+s(OrdDay,by=Lake,bs='cc')+s(Lake,bs="re"),data=df,family=nb)

My partial effect plots look like this:

library(gratia)
draw(hgam_gi)

enter image description here

My understanding of these partial effect plots is that the individual smooth LakeTBend doesn't differ much from the global smooth s(OrdDay) thus leadings to no effect in the plot and LakeFork has a stronger effect somewhere around OrdDay 200-250. I am wondering if there's a way to use an HGAM with a global smooth, and extract the individual smooths from that model similar to what you get with an HGAM without a global smooth.

With a model without a global smooth (I model according to Pederson et al. 2019), I can get a better understanding of the individual trends.

hgam_i<-gam(Depth~Lake+s(OrdDay,by=Lake,bs='cc'),data=df,family=nb)
draw(hgam_i)

enter image description here

My question is if there is a call I can make in gratia::draw() with an HGAM that can pull the partial effect plots independent of the global smooth? I know I can use predict.gam() to get the the individual predictions for each lake, which is great, but I really like the interpretability of partial effect plots and would love to find a way to run a single parsimonious model that displays global trends and individual trends.

Thanks!

Upvotes: 1

Views: 751

Answers (1)

Gavin Simpson
Gavin Simpson

Reputation: 174908

Your models are somewhat mis-specified - you have the group (lake) means in the model twice:

  1. the parametric Lake effect, and
  2. the random effect s(Lake, bs = "re")

you shouldn't include this effect twice; use one or the other.

Also, your first model really should have something for the by = Lake smooths to help the model identify those lake-specific differences from the global smooth. Typically we would add m = 1 to those by = Lake smooths when the model already contains a "global" smooth. You should also set knots when fitting a bs = "cc" smooth of day (assuming OrdDay is day of year?):

library("gratia")
library("mgcv")
library("dplyr")

knots <- list(OrdDay = c(0.5, 366.5))

m <- gam(Depth ~ Lake +
           s(OrdDay, bs = "cc") +
           s(OrdDay, by = Lake, bs = "cc", m = 1),
         data = df, family = nb, knots = knots)

another way of specifying this model is using an ordered factor instead of a normal factor. We have to set the contrasts differently so that the parametric term is easier to understand

df <- df |>
  mutate(oLake = ordered(Lake))
contrasts(df$oLake) <- contr.treatment

Now when you fit the model we use oLake in place of Lake and the first smooth of OrdDay represents the smooth for the reference level of oLake (here the reference is: Fork). The the by = oLake smooth (there is only one here as you only have two levels/lakes) represents the smooth difference between the two lakes, such that the smooth for TBend is the linear combination of the two smooths in the model

m2 <- gam(Depth ~ oLake +
            s(OrdDay, bs = "cc") +
            s(OrdDay, by = oLake, bs = "cc"),
          data = df, family = nb, knots = knots)

With either model, to get estimated smooths for your two lakes, you need to add combinations of smooths. With m (the global smooth plus m = 1 smooths for each lake) we need to form the following combinations:

  1. For Fork we want s(OrdDay) + s(OrdDay):LakeFork
  2. For TBend we want s(OrdDay) + s(OrdDay):LakeTBend

For m2 where we used ordered factors, we want:

  1. For Fork we want s(OrdDay)
  2. For TBend we want s(OrdDay) + s(OrdDay):oLakeTBend

You can choose to include the group effects or not. Including them would mean you are actually generating fitted values on the response scale which is not what draw() shows (it shows partial effects).

To get fitted values we just predict over the range of OrdDay for both lakes. Using the development version of {gratia} on GitHub we create a data slice and then use fitted_values() on the data slice to yield predicted/fitted values:

ds <- data_slice(m, OrdDay = evenly(OrdDay, n = 100), Lake = evenly(Lake))
fv <- fitted_values(m, data = ds, scale = "response")
fv

library("ggplot2")

fv |>
ggplot(aes(x = OrdDay, y = .fitted, group = Lake)) +
  geom_ribbon(aes(ymin = .lower_ci, ymax = .upper_ci, fill = Lake), 
    alpha = 0.2) +
  geom_line(aes(colour = Lake))

while for m2 we use

ds2 <- data_slice(m2, OrdDay = evenly(OrdDay, n = 100), oLake = evenly(oLake))
fv2 <- fitted_values(m2, data = ds2, scale = "response")
fv2

library("ggplot2")

fv2 |>
ggplot(aes(x = OrdDay, y = .fitted, group = oLake)) +
  geom_ribbon(aes(ymin = .lower_ci, ymax = .upper_ci, fill = oLake),
    alpha = 0.2) +
  geom_line(aes(colour = oLake)) +
  labs(fill = "Lake", colour = "Lake") +
  scale_colour_discrete() +
  scale_fill_discrete()

If you don't want the group means, so you are seeing something similar to that draw() produces (partial effects), you'll need to currently switch back to predict.gam(), using type = "terms", and then sum the rows over the relevant columns:

p <- predict(m2, newdata = ds2, type = "terms")

p contains

r$> head(p)                                                                     
      oLake s(OrdDay) s(OrdDay):oLakeTBend
1 0.0000000 0.2274295           0.00000000
2 0.3846254 0.2274295           0.07137302
3 0.0000000 0.2286423           0.00000000
4 0.3846254 0.2286423           0.07287456
5 0.0000000 0.2263583           0.00000000
6 0.3846254 0.2263583           0.07405621

so we only want to sum the rows over columns 2 and 3:

fv3 <- rowSums(p[, -1])

Getting the standard errors is a bit more difficult in this case, so if including the group means is OK with you, I'll omit the code to get the correct standard errors for the type = "terms" idea. Comment if you need that.

"terms" but with standard errors / credible intervals

If you want the "terms" idea but with standard errors/credible intervals, it's not easy (possible?) to do it directly with type = "terms" because we miss the covariances between terms. It is simpler to modify the fitted_values() approach (or use predict.gam()) to use the exclude or terms arguments to exclude specific terms or include specific terms in the predictions.

For m we just want the three smooths (1 global plus the two from the factor by smooth), so we would use:

fv <- fitted_values(m, data = ds, scale = "link",
  terms = c("s(OrdDay)", "s(OrdDay):LakeFork", "s(OrdDay):LakeTBend"))

and

fv2 <- fitted_values(m2, data = ds2, scale = "link",
  terms = c("s(OrdDay)", "s(OrdDay):oLakeTBend"))

noting that we are now predicting on the link (linear predictor) scale and because we are not using any parametric terms, these smooths are centred about the overall mean, as they would if we'd followed the "terms" example, but we get standard errors.

Note that I have updated the code above so that it uses the latest development version of gratia, where i have modified quite a few of the variable names. So be sure to install from GitHub, or you'll have to adjust the names for the fitted value (.fitted -> fitted) and credible interval (.lower_ci and .upper_ci -> lower and upper respectively).

Alternative way of fitting the GI model

m is a GI model, to use the terminology of Pedersen et al (2019). A better way of fitting this model, which wasn't available until recently in mgcv, is through the use of a new basis "sz":

# Alternative GI model
m3 <- gam(Depth ~ s(OrdDay, bs = "cc") +
           s(OrdDay, Lake, bs = "sz", xt = list(bs = "cc")),
         data = df, family = nb, knots = knots)

In this model formulation we note:

  1. There is no need to specify the parametric factor term for Lake as this is actually in the "sz" smooth,
  2. The first smooth s(OrdDay) is the average or global smooth, 3, The second smooth in the formula is the smooth differences between the global smooth and the individual effects of OrdDay in each Lake, including the group means,
  3. We have to pass the basis we want for the specific smooths, "cc" here, through to the constructor using the xt argument, because we need to use bs = "sz" to specify we want one of these constrained factor smooths,
  4. We are using Lake here, not oLake, so a regular factor,
  5. Partial effects plots that combine the global and the individual/group-specific smooths would be created using fitted_values() using exclude or terms as shown below

To get the partial effect plots, we can't now (as far as I know) exclude the difference in group means, but we can produce plots without the model constant term:

fv3 <- fitted_values(m3, data = ds, scale = "link",
  terms = c("s(OrdDay)", "s(OrdDay,Lake)"))

fv3 |>
ggplot(aes(x = OrdDay, y = .fitted, group = Lake)) +
  geom_ribbon(aes(ymin = .lower_ci, ymax = .upper_ci, fill = Lake), 
    alpha = 0.2) +
  geom_line(aes(colour = Lake))

Upvotes: 1

Related Questions