Оптимизация степеней свободы в сплайн-регрессии

У меня есть два набора данных о времени экспрессии генов:

Во-первых, экспрессию генов измеряли в течение 14 моментов времени из 4 групп:

df1 <- structure(list(val = c(-0.1, -0.13, -0.4, -0.3, -0.3, -0.2, -0.24, 
                            0.1, 0.2, 0.13, 0, 0.63, 0.83, 0.85, -0.07, -0.07, -0.27, -0.2, 
                            -0.2, -0.1, 0.2, 0.1, 0.07, 0.17, 0.6, 0.75, 1.1, 1.1, -0.13, 
                            -0.15, -0.26, -0.25, -0.14, 0.04, 0.2, 0.24, 0.23, 0.2, 0.1, 
                            0.73, 1, 1.3, 0, 0.06, -0.24, -0.17, -0.17, -0.04, 0.16, 0.1, 
                            0.14, 0.27, 0.34, 0.9, 0.97, 1.04), 
                    time = c(-1, 0, 1, 1.58,2, 2.58, 3, 3.32, 3.58, 4.17, 4.58, 5.58, 6.17, 7.39, 
                             -1, 0, 1, 1.58, 2, 2.58, 3, 3.32, 3.58, 4.17, 4.58, 5.58, 6.17, 7.39, 
                             -1, 0, 1, 1.58, 2, 2.58, 3, 3.32, 3.58, 4.17, 4.58, 5.58, 6.17,7.39, 
                             -1, 0, 1, 1.58, 2, 2.58, 3, 3.32, 3.58, 4.17, 4.58, 5.58,6.17, 7.39), 
                    group = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L,1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                                        2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,2L, 2L, 2L, 2L, 2L, 
                                        3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,3L, 3L, 3L, 
                                        4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,4L), 
                                      .Label = c("a", "b", "c", "d"), class = "factor")), .Names = c("val","time", "group"), 
               row.names = c(NA, -56L), class = "data.frame")


df1$group <- factor(df1$group,levels=c("a","b","c","d"))

который выглядит так (добавление сглаженной линии тренда loess):

library(ggplot2)
ggplot(df1,aes(x=time,y=val,color=group))+geom_point()+theme_minimal()+geom_smooth(se=F)+theme(legend.position="top",legend.title=element_blank())

Во-вторых, экспрессия генов измерялась в течение 14 одинаковых моментов времени, но теперь в двух разных группах, каждая из которых представлена представителями двух полов:

df2 <- structure(list(val = c(-0.23, -0.01, -0.14, -0.01, -0.21, -0.16, 
                       -0.24, -0.11, 0.02, -0.11, -0.01, -0.25, -0.47, -1.25, 0.02, 
                       -0.3, -0.02, 0.14, 0.25, -0.05, 0.15, 0.11, -0.24, -0.18, -0.39, 
                       -0.49, -0.5, -0.65, -0.06, 0.09, 0.1, 0.15, 0.08, 0.15, 0.4, 
                       0.24, 0.07, 0.08, -0.18, -0.35, -0.19, -0.81, -0.16, 0.29, -0.05, 
                       0.14, 0.14, 0.48, 0.34, 0.11, -0.07, -0.13, -0.41, -0.22, -0.54, 
                       -0.76, 0.35, 0.34, -0.06, 0.21, 0.14, 0.14, 0.25, 0.22, 0.25, 
                       0.16, 0.3, 0.44, 0.08, 0.48, 0.1, 0.16, -0.03, -0.22, 0.2, 0.01, 
                       -0.09, -0.02, -0.01, 0.06, -0.13, 0.19, 0.11, -0.04, -0.39, 0.03, 
                       -0.01, 0.09, 0.1, -0.14, -0.12, -0.1, 0.36, 0.08, 0.09, 0.09, 
                       0.42, 0.37, -0.14, 0.12, 0.09, 0.03, 0.06, -0.25, 0.2, -0.06, 
                       -0.44, 0.23, 0.03, 0.16, 0.81, 0.83),
               time = c(-1, 0, 1, 1.58,2, 2.58, 3, 3.32, 3.58, 4.17, 4.58, 5.58, 6.17, 7.39, 
                        -1, 0,1, 1.58, 2, 2.58, 3, 3.32, 3.58, 4.17, 4.58, 5.58, 6.17, 7.39, 
                        -1, 0, 1, 1.58, 2, 2.58, 3, 3.32, 3.58, 4.17, 4.58, 5.58, 6.17,7.39, 
                        -1, 0, 1, 1.58, 2, 2.58, 3, 3.32, 3.58, 4.17, 4.58, 5.58,6.17, 7.39, 
                        -1, 0, 1, 1.58, 2, 2.58, 3, 3.32, 3.58, 4.17, 4.58,5.58, 6.17, 7.39, 
                        -1, 0, 1, 1.58, 2, 2.58, 3, 3.32, 3.58, 4.17,4.58, 5.58, 6.17, 7.39, 
                        -1, 0, 1, 1.58, 2, 2.58, 3, 3.32, 3.58, 4.17, 4.58, 5.58, 6.17, 7.39, 
                        -1, 0, 1, 1.58, 2, 2.58, 3, 3.32, 3.58, 4.17, 4.58, 5.58, 6.17, 7.39), 
               sex = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                                 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
                                 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), 
                               .Label = c("F", "M"), class = "factor"), group = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                                                                                            2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), 
                                                                                          .Label = c("a", "b"), class = "factor")), .Names = c("val", "time", "sex", "group"), row.names = c(NA, -112L), class = "data.frame")
df2$sex <- ordered(df2$sex,levels=c("M","F"))

df2$group <- ordered(df2$group,levels=c("a","b"))

df2$col <- factor(paste0(df2$group,":",df2$sex))

который выглядит так (добавление линии тренда, сглаженной лёссом):

ggplot(df2,aes(x=time,y=val,color=col))+geom_point()+theme_minimal()+geom_smooth(se=F)+theme(legend.position="top",legend.title=element_blank())

Для df1 я хотел бы оценить влияние time на val с поправкой на group.

Для df2 я хотел бы оценить влияние time:group на val с поправкой на sex.

Глядя на данные, я подумал, что использование spline regressions будет уместным, поэтому я использовал функцию gam из пакета mgcv, которая, насколько я понимаю, оптимизирует степени свободы splines, подогнанных к данным.

Вот что я подогнал для df1:

mgcv1.fit <- mgcv::gam(val ~ group+s(time),data=df1)

Который дает:

Family: gaussian 
Link function: identity 

Formula:
val ~ group + s(time)

Estimated degrees of freedom:
7.18  total = 11.18 

GCV score: 0.01258176

Но 7,18 степеней свободы кажутся слишком большими для этих данных.

Для df2:

mgcv2.fit <- mgcv::gam(val ~ sex+s(time,by=group),data=df2)

который дает:

Family: gaussian 
Link function: identity 

Formula:
val ~ sex + s(time, by = group)

Estimated degrees of freedom:
1.72  total = 3.72 

GCV score: 0.08522094

Я предполагаю, что в этом случае я бы предположил, что степени свободы будут немного выше.

Еще один момент. Построение соответствующих значений для этих двух наборов данных:

df1$mgcv <- mgcv1.fit$fitted.values
ggplot(df1,aes(x=time,y=mgcv,color=group))+geom_point()+theme_minimal()+geom_smooth(se=F)+theme(legend.position="top",legend.title=element_blank())

что выглядит нормально.

Но для df2

df2$mgcv <- mgcv2.fit$fitted.values
ggplot(df2,aes(x=time,y=mgcv,color=col))+geom_point()+theme_minimal()+geom_smooth(se=F)+theme(legend.position="top",legend.title=element_blank())

Похоже, он перевернул ярлыки группы.

Итак, мои вопросы:

Правильно ли я использую mgcv::gam для оптимизации степеней свободы сплайна для моих вопросов?
Переупорядочивает ли mgcv образцы в своем fitted.values?

dan 18.08.2018 источник

comment

Я хочу изменить заголовок вопроса на что-то вроде: mgcv: by smooth, похоже, не возвращает правильное соответствие, потому что ваш вопрос в основном касается использования by. Кроме того, новое название, вероятно, упростит поиск для тех, кто борется с поиском. - Zheyuan Li 20.08.2018

Ответы (1)

arrow_upward
4
arrow_downward

Прежде всего, mgcv делает правильные вещи на уровнях факторов. Если вы отметите str(df2$sex), вы увидите, что «M» (мужской) — это первый уровень, а «F» (женский) — второй. Но из str(df2$col) кажется, что «F» стоит первым, поэтому при построении сюжета вы получаете неправильную маркировку.

Во-вторых, у вас не правильно указана вторая модель.

Сплайн s(time) находится под ограничением центрирования, когда нет переменной «by» или «by» является фактором. Таким образом, вы должны предоставить свою переменную "by" group в качестве отдельного термина в формуле вашей модели, чтобы уловить ее предельный эффект;
Поскольку переменная "by" group является упорядоченной переменной, mgcv применяет к ней контрасты, опуская первый уровень "a" при построении s(time, by = group). Поэтому вам нужно предоставить отдельный s(time) в качестве базового сглаживания.

Ваша текущая mgcv2.fit — довольно плохая модель (неудивительно), дающая объясненное отклонение 9%. Но если вы сделаете следующее, вы получите 64%.

gam(val ~ sex + s(time) + group + s(time, by = group), data = df2, method = "REML")

ggplot теперь выглядит правильно (я не менял df2$col, так что окраска все еще, вероятно, перевернута).

gam по умолчанию используется "GCV.Cp" в качестве метода выбора параметра сглаживания. Но рекомендуется использовать "REML", так как он менее подвержен переоснащению.

Примечание 1

Если переменная "by" group является (неупорядоченным) фактором, она не подлежит контрастам. Таким образом, формула модели должна быть:

val ~ sex + group + s(time, by = group)

Следующее цитируется из раздела переменные документа ?gam.models:

 If a ‘by’ variable is a ‘factor’ then it generates an indicator
 vector for each level of the factor, unless it is an ‘ordered’
 factor. In the non-ordered case, the model matrix for the smooth
 term is then replicated for each factor level, and each copy has
 its rows multiplied by the corresponding rows of its indicator
 variable. The smoothness penalties are also duplicated for each
 factor level.  In short a different smooth is generated for each
 factor level (the ‘id’ argument to ‘s’ and ‘te’ can be used to
 force all such smooths to have the same smoothing parameter).
 ‘ordered’ ‘by’ variables are handled in the same way, except that
 no smooth is generated for the first level of the ordered factor
 (see ‘b3’ example below).  This is useful for setting up
 identifiable models when the same smooth occurs more than once in
 a model, with different factor ‘by’ variables.

Примечание 2

Я не должен судить о вашей модели, но кажется, что между «F» и «M» существует четкая внутригрупповая разница. Из ваших данных мы видим, что «Ж» и «М» имеют большую разницу в группе «б», чем в группе «а». На данный момент эффект sex идентичен в обеих группах, это просто сдвиг по вертикали. Вы можете наблюдать это в приведенном выше ggplot в этом ответе. В конце концов, вам решать, какую модель выбрать, но на всякий случай, если вы хотите смоделировать это sex-group взаимодействие, вы можете сделать

df2$sex_group <- with(df2, interaction(sex, group))  ## the new variable is unordered
test <- gam(val ~ sex + group + s(time, by = sex_group), data = df2, method = "REML")