Adjusting p value with GAM in R?

Question

I am interested in understanding how X predicts three different Y values across four different age groups, controlling for Sex.

Some data for Y1:

structure(list(Y = c(0.2801485, 0.281114, 0.2556435, 0.2670325, 
0.2509095, 0.440258, 0.410253, 0.268027, 0.2641495, 0.287518, 
0.4374505, 0.370858, 0.3135055, 0.3216865, 0.2836315, 0.4176995, 
0.317726, 0.43616, 0.4236695, 0.4383975, 0.25266, 0.4136055, 
0.332767, 0.3132415, 0.3261405, 0.34355, 0.4005185, 0.3477285, 
0.3688955, 0.4653245, 0.3945735, 0.322322, 0.386704, 0.3627, 
0.4209505, 0.3291925, 0.3906035, 0.313632, 0.389333, 0.2887485, 
0.3224685, 0.4356595, 0.4760495, 0.4045595, 0.39168, 0.532351, 
0.3976815, 0.287972, 0.375473, 0.3411195, 0.352923, 0.4852365, 
0.413251, 0.5854595, 0.317848, 0.364348, 0.645034, 0.397519, 
0.339168, 0.3987795, 0.4516355, 0.29165, 0.3580975, 0.3960375, 
0.3760265, 0.345663, 0.3908205, 0.235979, 0.323781, 0.278155, 
0.3412205, 0.358712, 0.4360125, 0.3428525, 0.380136, 0.2429335, 
0.38803, 0.2917035, 0.450765, 0.3115525, 0.3234875, 0.4238855, 
0.4918315, 0.425098, 0.3816665, 0.32467, 0.300854, 0.333812, 
0.321521, 0.3045225, 0.398721, 0.2946885, 0.3402805, 0.3469465, 
0.5559045, 0.3365045, 0.413945, 0.3900765, 0.3724925, 0.3720195, 
0.441148, 0.5548215, 0.3664545, 0.307332, 0.305285, 0.348574, 
0.5294355, 0.434247, 0.3060095, 0.21842535, 0.323011, 0.319916, 
0.337996, 0.3816965, 0.4141275, 0.432481, 0.3867495, 0.4777715, 
0.39888, 0.328043, 0.290918, 0.205073, 0.376553, 0.3969375, 0.337341, 
0.2866875, 0.528681, 0.3649845, 0.4004325, 0.337676, 0.4870895, 
0.3818915, 0.30209, 0.5245945, 0.362414, 0.3280005, 0.5290235, 
0.384249, 0.453579, 0.466193, 0.2766955, 0.4699585, 0.4377835, 
0.5094885, 0.4556425, 0.251821, 0.3767055, 0.219233, 0.394665, 
0.402117, 0.346433, 0.424634, 0.3787175, 0.5182415, 0.357635, 
0.2883515, 0.38723, 0.32965, 0.474527, 0.3939005, 0.372789, 0.3428715, 
0.298425, 0.3394395, 0.486052, 0.3854875, 0.331162, 0.482598, 
0.3667065, 0.3231195, 0.293581, 0.2639825, 0.4304105, 0.4137305, 
0.3630545, 0.419433, 0.3816475, 0.311356, 0.440799, 0.368029, 
0.401345, 0.353489, 0.3324415, 0.394598, 0.2543955, 0.3698885, 
0.5063065, 0.3598435, 0.3013415, 0.2098135, 0.4068725, 0.3670155, 
0.4511215, 0.437022, 0.474631, 0.2720705, 0.4512285, 0.3314975, 
0.302051, 0.3499875, 0.487781, 0.4791155, 0.3049765, 0.389113, 
0.4181505, 0.373844, 0.291942, 0.439281, 0.2481589, 0.26381, 
0.3468945, 0.407657, 0.269458, 0.458682, 0.405109, 0.314457, 
0.401793, 0.4404885, 0.257801, 0.322888, 0.3306025, 0.4065515, 
0.4044915, 0.4647615, 0.45721, 0.409373, 0.2932925, 0.32179, 
0.366739, 0.3950615, 0.3352645, 0.3392295, 0.425147, 0.4731345, 
0.41028, 0.4798275, 0.388645, 0.27228195, 0.340866, 0.358758, 
0.3259775, 0.464307, 0.4323105, 0.56906, 0.2833945, 0.3703535, 
0.339217, 0.329044, 0.380571, 0.3088705, 0.21881395, 0.351079, 
0.4363515, 0.2193697, 0.4234945, 0.49115), X = c(44L, 37L, 35L, 
54L, 18L, 48L, 51L, 38L, 49L, 66L, 54L, 46L, 38L, 38L, 50L, 56L, 
56L, 61L, 59L, 56L, 41L, 25L, 37L, 62L, 55L, 32L, 49L, 48L, 50L, 
49L, 51L, 61L, 59L, 26L, 57L, 42L, 61L, 55L, 34L, 55L, 47L, 50L, 
55L, 65L, 50L, 41L, 58L, 26L, 52L, 52L, 49L, 50L, 58L, 70L, 48L, 
64L, 67L, 48L, 55L, 63L, 57L, 47L, 65L, 50L, 40L, 52L, 39L, 57L, 
49L, 39L, 59L, 52L, 61L, 66L, 52L, 59L, 53L, 54L, 41L, 17L, 45L, 
24L, 57L, 42L, 55L, 63L, 49L, 50L, 28L, 58L, 33L, 34L, 48L, 60L, 
68L, 52L, 46L, 62L, 38L, 49L, 35L, 33L, 54L, 55L, 58L, 45L, 62L, 
58L, 45L, 48L, 59L, 57L, 39L, 48L, 42L, 47L, 40L, 66L, 51L, 50L, 
52L, 52L, 39L, 58L, 46L, 48L, 55L, 64L, 67L, 45L, 47L, 61L, 40L, 
38L, 68L, 54L, 41L, 69L, 39L, 48L, 51L, 47L, 45L, 37L, 39L, 67L, 
41L, 61L, 56L, 49L, 52L, 40L, 45L, 57L, 36L, 47L, 49L, 43L, 67L, 
55L, 63L, 51L, 53L, 51L, 62L, 38L, 51L, 51L, 61L, 44L, 58L, 32L, 
56L, 49L, 41L, 33L, 55L, 45L, 32L, 51L, 65L, 60L, 56L, 55L, 49L, 
40L, 50L, 55L, 55L, 44L, 42L, 35L, 44L, 27L, 36L, 47L, 35L, 52L, 
41L, 48L, 45L, 45L, 55L, 51L, 43L, 61L, 62L, 49L, 32L, 36L, 56L, 
42L, 48L, 47L, 44L, 47L, 51L, 55L, 45L, 60L, 45L, 47L, 54L, 39L, 
50L, 52L, 51L, 41L, 36L, 58L, 49L, 64L, 40L, 38L, 55L, 55L, 49L, 
42L, 45L, 52L, 60L, 40L, 38L, 50L, 34L, 36L, 49L, 53L, 60L, 38L, 
43L, 57L, 32L, 56L, 50L, 65L), Age = c(77L, 65L, 79L, 9L, 34L, 
11L, 16L, 9L, 21L, 16L, 34L, 22L, 19L, 23L, 25L, 14L, 53L, 28L, 
22L, 22L, 21L, 82L, 81L, 16L, 19L, 77L, 15L, 18L, 15L, 78L, 24L, 
16L, 14L, 29L, 18L, 50L, 17L, 43L, 8L, 14L, 85L, 31L, 20L, 30L, 
23L, 78L, 29L, 61L, 14L, 22L, 10L, 83L, 15L, 13L, 15L, 15L, 29L, 
8L, 9L, 15L, 8L, 9L, 15L, 9L, 34L, 8L, 9L, 9L, 16L, 8L, 25L, 
21L, 23L, 13L, 56L, 10L, 27L, 8L, 8L, 8L, 8L, 80L, 80L, 15L, 
42L, 25L, 23L, 21L, 8L, 11L, 43L, 69L, 34L, 34L, 14L, 12L, 10L, 
22L, 78L, 16L, 76L, 12L, 10L, 16L, 13L, 66L, 11L, 26L, 12L, 16L, 
13L, 24L, 76L, 10L, 65L, 13L, 25L, 14L, 12L, 15L, 43L, 51L, 27L, 
15L, 24L, 34L, 63L, 17L, 15L, 9L, 12L, 17L, 82L, 75L, 24L, 44L, 
69L, 11L, 10L, 12L, 10L, 10L, 70L, 54L, 45L, 42L, 84L, 54L, 23L, 
23L, 14L, 81L, 17L, 42L, 44L, 16L, 15L, 43L, 45L, 50L, 53L, 23L, 
13L, 69L, 14L, 65L, 14L, 13L, 22L, 67L, 59L, 52L, 54L, 44L, 78L, 
62L, 69L, 10L, 63L, 57L, 22L, 12L, 62L, 9L, 82L, 53L, 54L, 66L, 
49L, 63L, 51L, 9L, 45L, 49L, 77L, 49L, 61L, 62L, 57L, 16L, 65L, 
75L, 45L, 16L, 55L, 17L, 64L, 67L, 56L, 52L, 63L, 10L, 62L, 14L, 
66L, 68L, 15L, 13L, 43L, 47L, 55L, 69L, 21L, 67L, 34L, 52L, 15L, 
31L, 64L, 55L, 44L, 13L, 71L, 64L, 13L, 25L, 34L, 50L, 61L, 70L, 
33L, 57L, 51L, 46L, 57L, 69L, 46L, 8L, 11L, 46L, 71L, 33L, 38L, 
17L, 29L, 28L), Age_group = structure(c(4L, 4L, 4L, 1L, 3L, 1L, 
1L, 1L, 2L, 1L, 3L, 2L, 2L, 2L, 2L, 1L, 3L, 2L, 2L, 2L, 2L, 4L, 
4L, 1L, 2L, 4L, 1L, 1L, 1L, 4L, 2L, 1L, 1L, 2L, 1L, 3L, 1L, 3L, 
1L, 1L, 4L, 3L, 2L, 3L, 2L, 4L, 2L, 4L, 1L, 2L, 1L, 4L, 1L, 1L, 
1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 1L, 1L, 1L, 1L, 1L, 
2L, 2L, 2L, 1L, 3L, 1L, 2L, 1L, 1L, 1L, 1L, 4L, 4L, 1L, 3L, 2L, 
2L, 2L, 1L, 1L, 3L, 4L, 3L, 3L, 1L, 1L, 1L, 2L, 4L, 1L, 4L, 1L, 
1L, 1L, 1L, 4L, 1L, 2L, 1L, 1L, 1L, 2L, 4L, 1L, 4L, 1L, 2L, 1L, 
1L, 1L, 3L, 3L, 2L, 1L, 2L, 3L, 4L, 1L, 1L, 1L, 1L, 1L, 4L, 4L, 
2L, 3L, 4L, 1L, 1L, 1L, 1L, 1L, 4L, 3L, 3L, 3L, 4L, 3L, 2L, 2L, 
1L, 4L, 1L, 3L, 3L, 1L, 1L, 3L, 3L, 3L, 3L, 2L, 1L, 4L, 1L, 4L, 
1L, 1L, 2L, 4L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 1L, 4L, 3L, 2L, 1L, 
4L, 1L, 4L, 3L, 3L, 4L, 3L, 4L, 3L, 1L, 3L, 3L, 4L, 3L, 4L, 4L, 
3L, 1L, 4L, 4L, 3L, 1L, 3L, 1L, 4L, 4L, 3L, 3L, 4L, 1L, 4L, 1L, 
4L, 4L, 1L, 1L, 3L, 3L, 3L, 4L, 2L, 4L, 3L, 3L, 1L, 3L, 4L, 3L, 
3L, 1L, 4L, 4L, 1L, 2L, 3L, 3L, 4L, 4L, 3L, 3L, 3L, 3L, 3L, 4L, 
3L, 1L, 1L, 3L, 4L, 3L, 3L, 1L, 2L, 2L), .Label = c("Adolescent", 
"Young", "Middle", "Older"), class = "factor"), Sex = c("Female", 
"Female", "Male", "Male", "Female", "Male", "Male", "Male", "Male", 
"Female", "Female", "Female", "Male", "Male", "Male", "Male", 
"Female", "Male", "Female", "Male", "Male", "Female", "Female", 
"Male", "Female", "Female", "Female", "Male", "Female", "Male", 
"Male", "Male", "Male", "Male", "Male", "Female", "Female", "Female", 
"Male", "Female", "Female", "Male", "Male", "Male", "Female", 
"Female", "Female", "Male", "Male", "Female", "Male", "Female", 
"Male", "Male", "Male", "Male", "Female", "Female", "Female", 
"Male", "Female", "Female", "Female", "Male", "Female", "Female", 
"Male", "Female", "Male", "Male", "Male", "Male", "Male", "Female", 
"Female", "Male", "Female", "Female", "Male", "Male", "Female", 
"Male", "Female", "Male", "Female", "Female", "Female", "Female", 
"Male", "Female", "Male", "Male", "Female", "Female", "Female", 
"Male", "Male", "Female", "Female", "Male", "Male", "Female", 
"Female", "Male", "Male", "Female", "Male", "Male", "Female", 
"Female", "Female", "Female", "Male", "Male", "Male", "Female", 
"Female", "Male", "Male", "Female", "Female", "Female", "Female", 
"Female", "Female", "Male", "Male", "Male", "Female", "Female", 
"Female", "Female", "Male", "Male", "Male", "Female", "Female", 
"Female", "Male", "Male", "Male", "Male", "Female", "Female", 
"Female", "Female", "Male", "Female", "Male", "Male", "Male", 
"Male", "Male", "Female", "Female", "Male", "Male", "Female", 
"Male", "Female", "Female", "Male", "Male", "Female", "Male", 
"Male", "Female", "Female", "Male", "Male", "Female", "Female", 
"Female", "Female", "Female", "Female", "Female", "Female", "Female", 
"Female", "Male", "Female", "Female", "Female", "Female", "Female", 
"Female", "Male", "Female", "Male", "Female", "Male", "Female", 
"Female", "Male", "Female", "Female", "Female", "Male", "Male", 
"Female", "Female", "Female", "Male", "Female", "Male", "Male", 
"Female", "Female", "Male", "Female", "Female", "Female", "Male", 
"Female", "Male", "Male", "Male", "Female", "Female", "Male", 
"Male", "Female", "Male", "Female", "Female", "Female", "Female", 
"Male", "Female", "Female", "Male", "Female", "Male", "Female", 
"Female", "Male", "Male", "Male", "Male", "Female", "Male", "Female", 
"Male", "Male", "Female", "Male", "Male", "Male", "Male", "Male", 
"Female", "Male", "Female", "Male", "Male")), row.names = c(NA, 
-256L), class = c("tbl_df", "tbl", "data.frame"), na.action = structure(c(`1` = 1L, 
`2` = 2L, `3` = 3L, `4` = 4L, `5` = 5L, `6` = 6L, `7` = 7L, `8` = 8L, 
`9` = 9L, `10` = 10L, `11` = 11L, `12` = 12L, `13` = 13L, `14` = 14L, 
`15` = 15L, `16` = 16L, `17` = 17L, `18` = 18L, `19` = 19L, `20` = 20L, 
`21` = 21L, `22` = 22L, `23` = 23L, `24` = 24L, `25` = 25L, `26` = 26L, 
`27` = 27L, `28` = 28L, `29` = 29L, `30` = 30L, `31` = 31L, `32` = 32L, 
`33` = 33L, `34` = 34L, `35` = 35L, `36` = 36L, `37` = 37L, `38` = 38L, 
`39` = 39L, `40` = 40L, `41` = 41L, `42` = 42L, `43` = 43L, `44` = 44L, 
`45` = 45L, `46` = 46L, `47` = 47L, `48` = 48L, `49` = 49L, `50` = 50L, 
`51` = 51L, `52` = 52L, `53` = 53L, `54` = 54L, `55` = 55L, `56` = 56L, 
`57` = 57L, `76` = 76L, `106` = 106L, `136` = 136L, `144` = 144L, 
`166` = 166L, `178` = 178L, `226` = 226L, `227` = 227L, `265` = 265L, 
`299` = 299L, `321` = 321L, `325` = 325L), class = "omit"))

If I then applied these GAM model in R:

mgcv::gam(Y1 ~ s(X, by = Age_Group)  + Sex, data = DF, method = "REML") 
mgcv::gam(Y2 ~ s(X, by = Age_Group)  + Sex, data = DF, method = "REML") 
mgcv::gam(Y3 ~ s(X, by = Age_Group)  + Sex, data = DF, method = "REML")

My questions are: 1a. Is that the correct GAM model syntax? 1b. I would like to know how to correct the p-values for multiple comparisons. P.adjust isn't properly working.

I don't often see GAM p value correction in my field publications, are there any exceptions that allow this?

There are no `Y1, Y2, Y3` in your sample data, though. – jay.sf Dec 24 '21 at 10:04 — jay.sf, Dec 24 '21 at 10:04

jay.sf · Answer 1 · 2021-12-24T10:04:08.583

Here is one way for a single model. You could try p.adjust() (with a small "p"). To get the initial p-values, use summary() or mgcv::summary.gam().

fit <- mgcv::gam(Y ~ s(X, by=Age_group) + Sex, data=DF, method="REML") 
sfit <- mgcv::summary.gam(fit)

str() reveals the structure of the output,

str(sfit)

where s table is stored,

(S <- sfit$s.table)
#                               edf   Ref.df         F    p-value
# s(X):Age_groupAdolescent 3.243795 4.053323 1.9116622 0.11045294
# s(X):Age_groupYoung      1.402075 1.706386 4.2091583 0.01378258
# s(X):Age_groupMiddle     1.157060 1.299352 0.3545537 0.72443981
# s(X):Age_groupOlder      1.000195 1.000391 1.1435619 0.28589857

and where p-values are stored.

(pv <- sfit$s.pv)  ## s p-values
# [1] 0.11045294 0.01378258 0.72443981 0.28589857

If you want to correct these p-values (where I'm not sure if it's right), they may be put into p.adjust() choosing one method (see ?p.adjust),

(p_adj <- p.adjust(pv, method='holm'))
# [1] 0.33135882 0.05513031 0.72443981 0.57179715

and perhaps cbinded to the summary.

cbind(S, p_adj)
#                               edf   Ref.df         F    p-value      p_adj
# s(X):Age_groupAdolescent 3.243795 4.053323 1.9116622 0.11045294 0.33135882
# s(X):Age_groupYoung      1.402075 1.706386 4.2091583 0.01378258 0.05513031
# s(X):Age_groupMiddle     1.157060 1.299352 0.3545537 0.72443981 0.72443981
# s(X):Age_groupOlder      1.000195 1.000391 1.1435619 0.28589857 0.57179715

However, as said, I'm not quite sure, whether p-values of a GAM should be corrected like this. The take-home-message of this answer is rather, how to correct any p-values for multiple comparisons. Concerning GAM, you maybe should take a look at this answer concerning GAM's on Cross Validated.

Thanks, yes I didnt add ''y2'' or ''y3'' as I didn't think it was necessary for getting the idea across, but can add them in if you think it's needed. I saw that GAM CV post before and it's the only topic I've found on it. Since that was pairwise adjustment I wasn't sure if it's a different situation that doesn't apply. — CanyonView, Dec 24 '21 at 16:37

Adjusting p value with GAM in R?

1 Answers1