Линейная регрессия с повторными измерениями в R

Я не смог понять, как выполнить линейную регрессию в R для повторного измерения проекта. В предыдущем вопросе (все еще без ответа) мне было предложено не использовать, lmа использовать смешанные модели. Я использовал lmследующим образом:

lm.velocity_vs_Velocity_response <- lm(Velocity_response~Velocity*Subject, data=mydata)

(более подробную информацию о наборе данных можно найти по ссылке выше)

Однако я не смог найти в Интернете ни одного примера с кодом R, показывающего, как выполнять анализ линейной регрессии.

То, что я хочу, - это, с одной стороны, график данных с линией, соответствующей данным, а с другой стороны, значение вместе со значением p для критерия значимости для модели. $R^2$

Есть ли кто-нибудь, кто может дать некоторые предложения? Любой пример кода R может быть очень полезным.

Редактировать В
соответствии с предложением, которое я получил до сих пор, решение для моего анализа моих данных, чтобы понять, существует ли линейная связь между двумя переменными Velocity_response (производными от вопросника) и Velocity (производными от производительности), должно быть следующим:

library(nlme)
summary(lme(Velocity_response ~ Velocity*Subject, data=scrd, random= ~1|Subject))

Результат резюме дает это:

    > summary(lme(Velocity_response ~ Velocity*Subject, data=scrd, random= ~1|Subject))
    Linear mixed-effects model fit by REML
     Data: scrd 
           AIC      BIC   logLik
      104.2542 126.1603 -30.1271

    Random effects:
     Formula: ~1 | Subject
            (Intercept) Residual
    StdDev:    2.833804 2.125353

Fixed effects: Velocity_response ~ Velocity * Subject 
                              Value Std.Error DF    t-value p-value
(Intercept)               -26.99558  25.82249 20 -1.0454288  0.3083
Velocity                   24.52675  19.28159 20  1.2720292  0.2180
SubjectSubject10           21.69377  27.18904  0  0.7978865     NaN
SubjectSubject11           11.31468  33.51749  0  0.3375754     NaN
SubjectSubject13           52.45966  53.96342  0  0.9721337     NaN
SubjectSubject2           -14.90571  34.16940  0 -0.4362299     NaN
SubjectSubject3            26.65853  29.41574  0  0.9062674     NaN
SubjectSubject6            37.28252  50.06033  0  0.7447517     NaN
SubjectSubject7            12.66581  26.58159  0  0.4764880     NaN
SubjectSubject8            14.28029  31.88142  0  0.4479188     NaN
SubjectSubject9             5.65504  34.54357  0  0.1637076     NaN
Velocity:SubjectSubject10 -11.89464  21.07070 20 -0.5645111  0.5787
Velocity:SubjectSubject11  -5.22544  27.68192 20 -0.1887672  0.8522
Velocity:SubjectSubject13 -41.06777  44.43318 20 -0.9242591  0.3664
Velocity:SubjectSubject2   11.53397  25.41780 20  0.4537754  0.6549
Velocity:SubjectSubject3  -19.47392  23.26966 20 -0.8368804  0.4125
Velocity:SubjectSubject6  -29.60138  41.47500 20 -0.7137162  0.4836
Velocity:SubjectSubject7   -6.85539  19.92271 20 -0.3440992  0.7344
Velocity:SubjectSubject8  -12.51390  22.54724 20 -0.5550080  0.5850
Velocity:SubjectSubject9   -2.22888  27.49938 20 -0.0810519  0.9362
 Correlation: 
                          (Intr) Velcty SbjS10 SbjS11 SbjS13 SbjcS2 SbjcS3 SbjcS6 SbjcS7 SbjcS8 SbjcS9 V:SS10 V:SS11 V:SS13 Vl:SS2 Vl:SS3
Velocity                  -0.993                                                                                                         
SubjectSubject10          -0.950  0.943                                                                                                  
SubjectSubject11          -0.770  0.765  0.732                                                                                           
SubjectSubject13          -0.479  0.475  0.454  0.369                                                                                    
SubjectSubject2           -0.756  0.751  0.718  0.582  0.362                                                                             
SubjectSubject3           -0.878  0.872  0.834  0.676  0.420  0.663                                                                      
SubjectSubject6           -0.516  0.512  0.490  0.397  0.247  0.390  0.453                                                               
SubjectSubject7           -0.971  0.965  0.923  0.748  0.465  0.734  0.853  0.501                                                        
SubjectSubject8           -0.810  0.804  0.769  0.624  0.388  0.612  0.711  0.418  0.787                                                 
SubjectSubject9           -0.748  0.742  0.710  0.576  0.358  0.565  0.656  0.386  0.726  0.605                                          
Velocity:SubjectSubject10  0.909 -0.915 -0.981 -0.700 -0.435 -0.687 -0.798 -0.469 -0.883 -0.736 -0.679                                   
Velocity:SubjectSubject11  0.692 -0.697 -0.657 -0.986 -0.331 -0.523 -0.607 -0.357 -0.672 -0.560 -0.517  0.637                            
Velocity:SubjectSubject13  0.431 -0.434 -0.409 -0.332 -0.996 -0.326 -0.378 -0.222 -0.419 -0.349 -0.322  0.397  0.302                     
Velocity:SubjectSubject2   0.753 -0.759 -0.715 -0.580 -0.360 -0.992 -0.661 -0.389 -0.732 -0.610 -0.563  0.694  0.528  0.329              
Velocity:SubjectSubject3   0.823 -0.829 -0.782 -0.634 -0.394 -0.622 -0.984 -0.424 -0.799 -0.667 -0.615  0.758  0.577  0.360  0.629       
Velocity:SubjectSubject6   0.462 -0.465 -0.438 -0.356 -0.221 -0.349 -0.405 -0.995 -0.449 -0.374 -0.345  0.425  0.324  0.202  0.353  0.385
Velocity:SubjectSubject7   0.961 -0.968 -0.913 -0.740 -0.460 -0.726 -0.844 -0.496 -0.986 -0.778 -0.718  0.886  0.674  0.420  0.734  0.802
Velocity:SubjectSubject8   0.849 -0.855 -0.807 -0.654 -0.406 -0.642 -0.746 -0.438 -0.825 -0.988 -0.635  0.783  0.596  0.371  0.649  0.709
Velocity:SubjectSubject9   0.696 -0.701 -0.661 -0.536 -0.333 -0.526 -0.611 -0.359 -0.676 -0.564 -0.990  0.642  0.488  0.304  0.532  0.581
                          Vl:SS6 Vl:SS7 Vl:SS8
Velocity                                      
SubjectSubject10                              
SubjectSubject11                              
SubjectSubject13                              
SubjectSubject2                               
SubjectSubject3                               
SubjectSubject6                               
SubjectSubject7                               
SubjectSubject8                               
SubjectSubject9                               
Velocity:SubjectSubject10                     
Velocity:SubjectSubject11                     
Velocity:SubjectSubject13                     
Velocity:SubjectSubject2                      
Velocity:SubjectSubject3                      
Velocity:SubjectSubject6                      
Velocity:SubjectSubject7   0.450              
Velocity:SubjectSubject8   0.398  0.828       
Velocity:SubjectSubject9   0.326  0.679  0.600

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.47194581 -0.46509026 -0.05537193  0.39069634  1.89436646 

Number of Observations: 40
Number of Groups: 10 
Warning message:
In pt(q, df, lower.tail, log.p) : NaNs produced
>

Теперь я не понимаю, где я могу получить R ^ 2 и соответствующие p-значения, указывающие на то, есть ли линейная зависимость между двумя переменными или нет, и я не понял, как мои данные могут быть нанесены с помощью линии, соответствующей регрессия.

Кто-нибудь может быть так добр, чтобы просветить меня? Мне очень нужна ваша помощь, ребята ...

r regression mixed-model repeated-measures L_T
источник

«Модели смешанных эффектов и расширения в экологии с помощью R» Zuur et al. является хорошим введением в линейные модели смешанных эффектов, которые фокусируются не столько на теории, сколько на применении методологии.

Роланд

Дорогой Роланд, я считаю, что эта книга полезна, но я скорее ищу что-то в Интернете ... У вас есть какая-нибудь веб-страница, чтобы предложить?

L_T

Как я уже говорил в вашем предыдущем посте, с lm () связан сюжет. Итак, если ваша модель M1, вы можете использовать plot (M1).

Питер Флом - Восстановить Монику

Уважаемый @PeterFlom, да, но вы также сказали мне избегать использования lm для повторных измерений. Итак, мой вопрос, должен ли я использовать lm для анализа своих данных или другой функции. Любое предложение?

L_T

Как я уже сказал, посмотрите на многоуровневые модели. В R вы можете посмотреть на nlmeпакет. Также поищите эту тему на этом сайте, здесь много написано об этом.

Питер Флом - Восстановить Монику

Что вы делаете на самом деле, зависит от целей анализа. Я не уверен точно, каковы цели вашего анализа, но я рассмотрю несколько примеров, и, надеюсь, один из них будет применим к вашей ситуации.

Случай 1 : Одна количественная переменная, измеренная дважды

Допустим, вы провели исследование на людях, в котором участники дважды проходили статистический тест, и вы хотели выяснить, отличаются ли средние оценки по второму измерению от первого (чтобы определить, произошло ли обучение). Если оценки test1 и test2 хранятся во фрейме данных d, вы можете сделать это полностью, используя функцию lm (), например:

mod <- lm(test2 - test1 ~ 1, data = d)
summary(mod)

Тест перехвата - это тест на разницу между test1 и test2. Обратите внимание, что у вас не будет delta-R ^ 2 для разницы между test1 и test2 - вместо этого ваша мера величины эффекта должна быть чем-то вроде коэна d.

Случай 2a : одна количественная переменная, измеренная дважды, одна дихотомическая переменная, измеренная полностью между субъектами

Допустим, у нас одинаковый дизайн исследования, но мы хотим знать, имели ли место разные уровни обучения для мужчин и женщин. Итак, у нас есть одна количественная переменная (производительность теста), которая измеряется дважды, и одна дихотомическая переменная, измеряемая один раз. Предполагая, что test1, test2 и пол все содержатся во фрейме данных d, мы также можем протестировать эту модель только с помощью lm (), как в:

mod <- lm(test2 - test1 ~ gender, data = d)
summary(mod)
lm.sumSquares(mod) # lm.sumSquares() is located in the lmSupport package, and gives the change in R^2 due to the between-subjects part of the model

Предполагая, что пол центрирован (то есть закодирован, например, мужчина = -5 и женщина = +.5), перехват в этой модели является тестом разницы между тестом 1 и тестом 2, усредненным по мужчинам и женщинам. Коэффициент для пола - это взаимодействие между временем и полом. Чтобы получить усредненный по времени эффект пола, вам нужно сделать:

mod <- lm(rowMeans(cbind(test2, test1)) ~ gender, data = d)
summary(mod)

Случай 2b : Одна количественная переменная, измеренная дважды, одна количественная переменная, измеренная только один раз

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

d$interestc <- d$interest - mean(d$interest)

Предполагая, что test1, test2 и Interestc все находятся в кадре данных d, этот вопрос может быть протестирован очень похоже на Случай 1a:

mod <- lm(test2 - test1 ~ interestc, data = d)
summary(mod)
lm.sumSquares(mod)

Еще раз, перехват в этой модели проверяет, изменились ли в среднем по интересам результаты тестов со времени 1 на время 2. Однако эта интерпретация верна только тогда, когда интерес сосредоточен. Коэффициент для интереса - зависит ли влияние времени от базового интереса. Мы могли бы получить эффект интереса, усредненный по времени, усреднив вместе test1 и test 2, как указано выше, и протестировав эффект интереса для этой составной переменной.

Случай 2c : одна количественная переменная, измеренная дважды, одна категориальная переменная, измеренная только один раз

Давайте предположим, что ваша переменная между субъектами была категорией, измеренной только один раз. Так, например, предположим, что вас интересовало, имели ли люди разных рас (белые против азиатов против черных против испаноязычных) разные уровни обучения от времени 1 к времени 2. Предположим, что test1, test2 и race находятся в кадре данных d Вам сначала нужно противопоставить код гонке. Это можно сделать, используя запланированные ортогональные контрасты, фиктивные коды или коды эффектов, в зависимости от конкретных гипотез / вопросов, которые вы хотите проверить (я рекомендую обратиться к lm.setContrasts (), если вы ищете вспомогательную функцию для этого) , Предполагая, что переменная гонки уже контрастно закодирована, вы должны использовать lm () очень похоже на два вышеупомянутых случая, как в:

mod <- lm(test2 - test1 ~ race, data = d)
summary(mod)
lm.sumSquares(mod)

Если предположить, что контрасты расы центрированы, то пересечение в этой модели снова является «главным эффектом» времени. Коэффициенты для контрастов расы - это взаимодействия между этими контрастами и временем. Чтобы получить омнибусные эффекты расы, используйте следующий код:

Anova(mod, type = 3)

Случай 3 : Одна количественная переменная, измеренная 3 раза (т.е. трехуровневая манипуляция внутри субъекта)

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

Например, предположим, что ваша основная гипотеза состоит в том, что результаты тестов будут линейно увеличиваться со времени 1 во время 3. Предполагая, что test1, test2 и test3 находятся во фрейме данных d, эту гипотезу можно проверить, сначала создав следующую композицию:

d$lin <- d[, paste("test", sep = "", 1:3)] %*% c(-1, 0, 1)

Затем вы проверите, имеет ли модель только перехват, использующая lin в качестве зависимой переменной, перехват, отличный от 0, как в:

mod <- lm(lin ~ 1, data = d)
summary(mod)

Это даст вам возможность проверить, увеличивались ли показатели статистики с течением времени. Конечно, вы можете создавать другие типы пользовательских оценок различий, в зависимости от ваших конкретных гипотез.

Если вы заботитесь о значимых омнибусных тестах, вам нужно использовать функцию Anova () из автомобильного пакета. Конкретная реализация немного запутана. По сути, вы указываете, какие переменные находятся внутри субъекта, а какие - между субъектами, используя lm (). Затем вы создаете внутреннюю часть модели (то есть указываете, какие из test1, test2 и test3 были измерены первыми, вторыми и третьими) и затем передаете эту модель в Anova (), создавая фрейм данных с именем idata. Используя мой гипотетический пример:

mod <- lm(cbind(test1, test2, test3) ~ 1, data = d) # No between-subjects portion of the model
idata <- data.frame(time = c("time1", "time2", "time3")) # Specify the within-subjects portion of the model
mod.A <- Anova(mod, idata = idata, idesign = ~time, type = 3) # Gives multivariate tests.  For univariate tests, add multivariate = FALSE
summary(mod.A)

Оператор idesign указывает Anova включить в модель переменную времени (составленную из test1, test2 и test3). Этот код даст вам ваши сводные тесты влияния времени на результаты тестов.

Случай 4 : Одна количественная переменная, измеренная 3 раза, одна количественная переменная между субъектами

Этот случай является простым расширением варианта 3. Как и выше, если вы заботитесь только о тестах на 1 степень свободы, вы можете просто создать пользовательский показатель разницы с помощью переменной внутри объекта. Итак, если предположить, что test1, test2, test3 и интерес все находятся во фрейме данных d, и предположим, что нас интересует линейное влияние времени на результаты тестов (и как эти эффекты времени меняются в зависимости от базового интереса), вы должны сделать следующее:

d$lin <- d[, paste("test", sep = "", 1:3)] %*% c(-1, 0, 1)

Затем сделайте следующее (со средним интересом):

mod <- lm(lin ~ interestc, data = d)
summary(mod)
lm.sumSquares(mod)

Если вы хотите провести комплексный тест, сделайте следующее:

mod <- lm(cbind(test1, test2, test3) ~ interest, data = d) # We now have a between-subjects portion of the model
idata <- data.frame(time = c("time1", "time2", "time3"))
mod.A <- Anova(mod, idata = idata, idesign = ~time * interest, type = 3) # The idesign statement assumes that we're interested in the interaction between time and interest
summary(mod.A)

Другие случаи: я опущу их для краткости, но они являются простыми расширениями того, что я уже описал.

Обратите внимание, что (одномерные) омнибусные тесты времени, когда время имеет более 2 уровней, предполагают сферичность. Это предположение становится довольно несостоятельным по мере увеличения количества уровней. Если у вас есть несколько точек измерения в вашем дизайне (скажем, 4+), я настоятельно рекомендую вам использовать что-то вроде многоуровневого моделирования и перейти к пакету, специализированному для этой техники (например, nlme или lme4) .

Надеюсь это поможет!

Патрик С. Форшер
источник

Уважаемый Патрик @ user1188407, спасибо вам большое за такой ответ. К сожалению, мой случай, вероятно, соответствует тому, что вы написали в последних предложениях ... поэтому мне понадобится пример кода R, чтобы понять, как обрабатывать мои данные. Действительно, если вы посмотрите на схему моего эксперимента, описанную в предыдущем посте stackoverflow.com/questions/12182373/…, вы увидите, что у меня есть переменная, измеренная 4 раза (т.е. скорость, измеренная в 4 условиях)

L_T

и я хочу выяснить, существует ли линейная зависимость с переменной (speed_response), выражающей воспринимаемую скорость в четырех условиях. Таким образом, каждый участник прошел 4 условия, а затем оценил восприятие этих условий. Я хочу знать, связано ли представление с восприятием ...

L_T

Хорошо, если вы хотите использовать многоуровневое решение для моделирования, вы можете использовать множество различных онлайн-ресурсов. Для начала взгляните на пакет nlme и эту виньетку . Виньетка немного устарела (2002), я нашел ее полезной, когда изучал многоуровневое моделирование. Наконец, вы можете проверить книгу, опубликованную создателями пакета nlme.

Патрик С. Форшер

Уважаемый Патрик @ user1188407 спасибо. Я изучил многоуровневые модели и пришел к этой формуле, чтобы проанализировать свои данные: lme (Velocity_response ~ Velocity * Subject, data = scrd, random = ~ 1 | Subject) Не могли бы вы подтвердить, что эта формула верна для анализа? хотите выполнить на моих данных? Однако я не понимаю, как я могу получить значения R ^ 2 и p, а также как построить график с точками и линией, соответствующей регрессии. Не могли бы вы мне помочь? Я не статик ...

L_T

Формула кажется мне правильной, основываясь на моем понимании вашего исследования (и при условии, что вы отформатировали свои данные в формате периода личности). Вы получили бы свои p-значения, сохранив результаты вашего анализа в объект (как я это делаю в моих примерах) и получив сводку этого объекта. Однако из-за различий между многоуровневыми моделями и традиционной регрессией (например, в метриках размера эффекта - в многоуровневых моделях нет стандартной метрики), я настоятельно рекомендую вам подробнее ознакомиться с этой техникой перед ее использованием. Похоже, что другие пользователи рекомендовали несколько хороших вариантов.

Патрик С. Форшер

Линейная регрессия с повторными измерениями в R

Ответы: