Форум врачей-аспирантов

Здравствуйте, гость ( Вход | Регистрация )

> Построение регрессионной модели с нелинейной зависимостью от времени, различающийся в зависимости от категориального предиктора, Кривые должны стартовать из одной точки
ИНО
сообщение 26.04.2025 - 19:23
Сообщение #1





Группа: Пользователи
Сообщений: 290
Регистрация: 1.06.2022
Из: Донецк
Пользователь №: 39632



Столкнулся с проблемой. Имеются две группы животных, делая реверанс в сторону основной тематики форума, назовем их "лечение" и "контроль". В определенные дни из каждой группы отбиралось по одной особи и провождились замеры некого параметра. При этом особь изымалась безвозвратно (секир башка smile.gif), т. е. можно с легким сердцем считать, что все наблюдения независимы друг от друга. Из теоретических соображений, что в о временной точке 0, перед началом "лечения", различий между группами не было, но со временем появились. Интерес заключен в различии динамик некого параметра между "леченными" и "нелеченными".

Создадим тестовый набор данных, "приближенных к боевым":

Код
testdata<-data.frame(group=factor(c(rep("A", 7),rep("B", 6))), time=c(15, 17, 18, 19, 20, 21, 22, 15, 18, 19, 20, 21, 22), response=c(1.078, 0.949, 0.793, -2.201, -1.181, -1.95, -2.391, 2.628, 2.525, 2.852, 2.811, 3.484, 2.206))


Отобразим их на диаграмме рассеяния:

Код
plot(response~time, data=testdata, col=group, xlim=c(0,22))


Такие границы оси абсцисс выбраны, дабы нулевая временная точка уместилась, хотя вблизи нее данные собраны не были. Несмотря на малый объем выборки, различия между группами на глаз вполне явственно. Но хочется подогнать модель, которая и кривые с ДИ нарисует, и гипотезы проверит. Важное условие - кривые для обеих групп должны выходить из одной точки с абсциссой 0. Увы, я смог додуматься только до того, как этого добиться с обычной линейной регрессией:

Код
time<-seq(0, 22, 0.1)
groupA<-data.frame(group=rep(testdata$group[1], 221), time=time)
groupB<-data.frame(group=rep(testdata$group[8], 221), time=time)

predict(m1, newdata=groupA)[1]==predict(m1, newdata=groupB)[1]


Цитата
TRUE


Таким образом, вышеозначенное условие выполнено: прямые выходят из одной точки. Проблема в том, что они прямые, а зависимость явно криволинейная. Поэтому пересечение по оси ординат лежит у черта на рогах:

Код
predict(m1, newdata=groupA)[1]


Цитата
6.533422


На диаграмме не уместилась:

Код
lines(time, predict(m1, newdata=groupA))
lines(time, predict(m1, newdata=groupB), col="red")


См. сплошные линии на прикрепленном рисунке. Из теоретических соображений относительно изучаемого процесса, это полный нонсенс. Предположительно, истинная точка старта лежит по оси ординат где-то между 0 и 1. Я накарябал от руки пунктирные линии, которые примерно отображают характер зависимости, какой она мне видится. Увы повторить подобное с помощью, танцев вокруг lm() с полиномами и различными сплайнами не смог - кривые пересекаются где угодно, но только не в нулевой временной точке. В настоящей нелинейной регрессии типа nls() я ни в зуб ногой, да и сомнительно, что в принципе возможно определить истинный математический закон по по столь малым данным. И не нужно оно мне в принципе. Подойдут любые кривые б. м. похожие на мои пунктиры, самое главное - чтоб они плавно расходилась из нулевой временной точки, и сама эта точка отстояла по оси абсцисс не шибко далеко от интервала (0; 1). Также требуется возможность проверки гипотезы о взаимодействии фактора принадлежности группе со временем и построения ДИ для пересказанных значений.

Сообщение отредактировал ИНО - 26.04.2025 - 19:24
Эскизы прикрепленных изображений
Прикрепленное изображение
 
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
 
Открыть тему
Ответов
ИНО
сообщение 28.04.2025 - 18:29
Сообщение #2





Группа: Пользователи
Сообщений: 290
Регистрация: 1.06.2022
Из: Донецк
Пользователь №: 39632



Кажись, додумался. Надо построить две модели: одну полную, а вторую такую же, но вообще без участия фактора группы, после чего сделать anova(модель1, модель 2).

Код
m_1<-lm(response~time:group+I(time^3):group,data=testdata)
m_2<-lm(response~time+I(time^3),data=testdata)
anova(m_1, m_2)


Верно?
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 

Сообщений в этой теме


Добавить ответ в эту темуОткрыть тему