Regresja liniowa z powtarzanymi pomiarami w R

12

Nie byłem w stanie wymyślić, jak przeprowadzić regresję liniową w R dla projektu z powtarzanymi pomiarami. W poprzednim pytaniu (wciąż bez odpowiedzi) zasugerowano mi, aby nie używać, lmale raczej używać modeli mieszanych. Użyłem lmw następujący sposób:

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

(więcej szczegółów na temat zestawu danych można znaleźć pod linkiem powyżej)

Nie udało mi się jednak znaleźć w Internecie żadnego przykładu z kodem R pokazującym, jak przeprowadzić analizę regresji liniowej.

Co chcę jest z jednej strony przedstawia wykres danych z linii armatury dane, az drugiej strony wartość wraz z wartością p dla testu istotności dla modelu.R2

Czy jest ktoś, kto może podać jakieś sugestie? Każdy przykład kodu R może być bardzo pomocny.


Edytuj
Zgodnie z sugestią, którą otrzymałem do tej pory, rozwiązanie mojej analizy moich danych w celu zrozumienia, czy istnieje liniowa zależność między dwiema zmiennymi Velocity_response (pochodzącymi z kwestionariusza) i Velocity (pochodzącymi z wyników) powinno być następujące:

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

Wynik podsumowania daje to:

    > 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
> 

Teraz nie rozumiem, gdzie mogę uzyskać R ^ 2 i odpowiadające im wartości p wskazujące, że nie ma liniowej zależności między dwiema zmiennymi, ani nie rozumiem, jak moje dane można wykreślić z linią pasującą do regresja.

Czy ktoś może być tak miły, aby mnie oświecić? Naprawdę potrzebuję waszej pomocy chłopaki ...

L_T
źródło
„Modele i rozszerzenia efektów mieszanych w ekologii z R” Zuur i in. to miłe wprowadzenie do liniowych modeli efektów mieszanych, które mniej koncentrują się na teorii, a bardziej na zastosowaniu metodologii.
Roland
Drogi Roland, uważam, że ta książka jest przydatna, ale raczej szukam czegoś w Internecie ... czy masz jakąś stronę internetową do zasugerowania?
L_T
1
Jak powiedziałem w twoim wcześniejszym poście, lm () ma z tym powiązany wykres. Jeśli więc masz model M1, możesz użyć wykresu (M1).
Peter Flom - Przywróć Monikę
drodzy @PeterFlom tak, ale powiedziałeś mi również, aby unikać używania lm do projektowania powtarzanych pomiarów. Moje pytanie brzmi więc, czy muszę użyć lm do analizy moich danych lub innej funkcji. Jakieś sugestie?
L_T
1
Tak jak powiedziałem, spójrz na modele wielopoziomowe. W R możesz spojrzeć na nlmepakiet. Przeszukaj również tę stronę, aby znaleźć temat, wiele o nim napisano.
Peter Flom - Przywróć Monikę

Odpowiedzi:

17

To, co robisz, zależy od celów analizy. Nie jestem do końca pewien, jakie są cele twojej analizy, ale przejdę przez kilka przykładów i mam nadzieję, że jeden z nich będzie miał zastosowanie w twojej sytuacji.

Przypadek 1 : Jedna zmienna ilościowa zmierzona dwukrotnie

Powiedzmy, że przeprowadziłeś badanie na ludziach, w którym uczestnicy dwukrotnie wzięli test statystyczny i chciałeś dowiedzieć się, czy średnie wyniki drugiego pomiaru różniły się od pierwszego pomiaru (w celu ustalenia, czy nastąpiło uczenie się). Jeśli wyniki test1 i test2 są przechowywane w ramce danych d, możesz to zrobić całkowicie za pomocą funkcji lm (), jak w:

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

Test przechwytywania jest testem różnicy między testem 1 i testem 2. Zauważ, że nie będziesz miał delta-R ^ 2 dla różnicy między test1 i test2 - zamiast tego, twoja miara wielkości efektu powinna być podobna do d cohena d.

Przypadek 2a : Jedna zmienna ilościowa zmierzona dwukrotnie, jedna zmienna dychotomiczna, zmierzona całkowicie między osobnikami

Powiedzmy, że mamy taki sam plan badań, ale chcemy wiedzieć, czy dla mężczyzn i kobiet występują różne wskaźniki uczenia się. Mamy więc jedną zmienną ilościową (wydajność testu), która jest mierzona dwukrotnie, i jedną zmienną dychotomiczną, mierzoną raz. Zakładając, że test1, test2 i płeć są zawarte w ramce danych d, moglibyśmy również przetestować ten model tylko za pomocą lm (), jak w:

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

Zakładając, że płeć jest wyśrodkowana (tj. Zakodowana, na przykład: mężczyzna = -5, a kobieta = +.5), punkt przecięcia w tym modelu jest testem różnicy między testem 1 a testem 2, uśrednionym dla mężczyzn i kobiet. Współczynnikiem dla płci jest interakcja między czasem a płcią. Aby uzyskać efekt płci uśredniony w czasie, musisz wykonać:

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

Przypadek 2b : Jedna zmienna ilościowa zmierzona dwukrotnie, jedna zmienna ilościowa, zmierzona tylko raz

Załóżmy, że znowu mamy jedną zmienną ilościową zmierzoną dwukrotnie i jedną zmienną ilościową zmierzoną raz. Załóżmy na przykład, że mieliśmy miarę podstawowego zainteresowania statystykami i chcieliśmy ustalić, czy ludzie, którzy mieli wyższy poziom podstawowego zainteresowania, uczyli się więcej od 1 do 2 czasu. Najpierw musielibyśmy skoncentrować zainteresowanie, jak w :

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

Zakładając, że test1, test2 i odsetki znajdują się w ramce danych d, pytanie to można następnie przetestować bardzo podobnie do przypadku 1a:

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

Po raz kolejny punkt przecięcia w tym modelu sprawdza, czy uśrednione dla odsetka wyniki testu zmieniły się od czasu 1 do czasu 2. Jednak ta interpretacja obowiązuje tylko wtedy, gdy odsetki są wyśrodkowane. Współczynnik oprocentowania zależy od tego, czy wpływ czasu zależy od oprocentowania podstawowego. Możemy uzyskać efekt zainteresowania, uśredniony w czasie, poprzez uśrednienie razem testu 1 i testu 2, jak powyżej, i przetestowanie wpływu zainteresowania na tę zmienną złożoną.

Przypadek 2c : Jedna zmienna ilościowa zmierzona dwukrotnie, jedna zmienna kategorialna, zmierzona tylko raz

Załóżmy, że twoja zmienna między podmiotami była kategorią mierzoną tylko raz. Załóżmy na przykład, że interesowałeś się tym, czy ludzie różnych ras (biały kontra azjatycki kontra czarny kontra hiszpanie) mieli różną wiedzę w zakresie od 1 do 2. Zakładając, że test1, test2 i rasa znajdują się w ramce danych d , najpierw trzeba by kontrastować wyścig kodu. Można to zrobić za pomocą planowanych kontrastów ortogonalnych, kodów fikcyjnych lub kodów efektów, w zależności od konkretnych hipotez / pytań, które chcesz przetestować (polecam skorzystanie z lm.setContrasts (), jeśli szukasz funkcji pomocnika, aby to zrobić) . Zakładając, że zmienna wyścigowa jest już zakodowana kontrastowo, użyłbyś lm () bardzo podobnie do powyższych dwóch przypadków, jak w:

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

Zakładając, że kontrasty rasowe są wyśrodkowane, przecięcie w tym modelu jest po raz kolejny „głównym efektem” czasu. Współczynniki dla kontrastów rasowych to interakcje między tymi kontrastami i czasem. Aby uzyskać omnibus efekty rasy, użyj następującego kodu:

Anova(mod, type = 3)

Przypadek 3 : Jedna zmienna ilościowa zmierzona 3 razy (tj. Trzypoziomowa manipulacja wewnątrz osobników)

Załóżmy, że dodałeś trzeci projekt do projektu od pierwszego przypadku. Tak więc twoi uczestnicy wzięli test statystyczny trzy razy zamiast dwa razy. Tutaj masz kilka możliwości, w zależności od tego, czy chcesz przeprowadzić test zbiorczy różnic między punktami czasowymi (czasami nie).

Załóżmy na przykład, że twoja główna hipoteza jest taka, że ​​wyniki testu będą liniowo rosły od czasu 1 do czasu 3. Zakładając, że test1, test2 i test3 znajdują się w ramce danych d, hipotezę tę można przetestować, najpierw tworząc następujący kompozyt:

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

Następnie sprawdziłbyś, czy model tylko przechwytujący wykorzystujący lin jako zmienną zależną ma przecięcie inne niż 0, jak w:

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

To da ci test, czy wyniki statystyk rosły z czasem. Możesz oczywiście tworzyć inne typy niestandardowych wyników różnic, w zależności od konkretnych hipotez.

Jeśli zależy Ci na testach zbiorczych o znaczeniu, musisz użyć funkcji Anova () z pakietu samochodowego. Konkretne wdrożenie jest nieco skomplikowane. Zasadniczo za pomocą lm () określasz, które zmienne są między podmiotami, a które między nimi. Następnie tworzysz część modelu w ramach podmiotów (tj. Określasz, które z testów1, test2 i test3 zostały zmierzone najpierw, drugie i trzecie), a następnie przekazujesz ten model do Anova (), tworząc ramkę danych o nazwie idata. Korzystając z mojego hipotetycznego przykładu:

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)

Instrukcja idesign mówi Anova, aby w modelu uwzględniła zmienną czasową (złożoną z test1, test2 i test3). Ten kod da ci wszechstronne testy wpływu czasu na wyniki testu.

Przypadek 4 : Jedna zmienna ilościowa zmierzona 3 razy, jedna międzyosobnicza zmienna ilościowa

Ten przypadek jest prostym rozszerzeniem Przypadku 3. Jak wyżej, jeśli zależy ci tylko na 1 stopniu swobody testów, możesz po prostu stworzyć niestandardowy wynik różnicy za pomocą zmiennej wewnątrzosobniczej. Zakładając, że test1, test2, test3 i odsetki znajdują się w ramce danych d, i zakładając, że jesteśmy zainteresowani liniowym wpływem czasu na wyniki testu (oraz tym, jak te efekty czasu różnią się w zależności od poziomu odniesienia), zrobiłbyś następujące:

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

Następnie wykonaj następujące czynności (ze średnim zainteresowaniem):

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

Jeśli chcesz przeprowadzić testy zbiorcze, wykonaj następujące czynności:

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)

Inne przypadki: pominię je dla zwięzłości, ale są to proste rozszerzenia tego, co już opisałem.

Należy pamiętać, że wszystkie (jednoczynnikowe) omnibusowe testy czasu, w których czas ma więcej niż 2 poziomy, zakładają sferyczność. To założenie staje się nie do utrzymania, gdy zwiększasz liczbę poziomów. Jeśli masz całkiem sporo punktów pomiarowych w swoim projekcie (powiedzmy 4+), zdecydowanie zalecamy użycie czegoś takiego jak modelowanie wielopoziomowe i przejście do pakietu specjalizującego się w tej technice (takiego jak nlme lub lme4 .

Mam nadzieję że to pomoże!

Patrick S. Forscher
źródło
Drogi Patricku @ user1188407, dziękuję bardzo za udzielenie takiej odpowiedzi. Niestety moja sprawa prawdopodobnie pasuje do tego, co napisałeś w ostatnich zdaniach ... więc potrzebowałbym przykładu kodu R, aby zrozumieć, jak traktować moje dane. Rzeczywiście, jeśli spojrzysz na projekt mojego eksperymentu opisany w poprzednim poście stackoverflow.com/questions/12182373/ ... możesz zobaczyć, że mam zmienną zmierzoną 4 razy (tj. Prędkość zmierzoną w 4 warunkach)
L_T
i chcę sprawdzić, czy istnieje relacja liniowa ze zmienną (odpowiedź_wysokości) wyrażającą postrzeganą prędkość w czterech warunkach. Tak więc każdy uczestnik przeszedł 4 warunki, a następnie ocenił ich postrzeganie. Chcę wiedzieć, czy przedstawienie jest powiązane z percepcją ...
L_T
Jeśli chcesz skorzystać z wielopoziomowego rozwiązania do modelowania, możesz skorzystać z wielu różnych zasobów online. Na początek powinieneś rzucić okiem na pakiet nlme i tę winietę . Winieta jest nieco nieaktualna (2002), uznałem ją za przydatną, kiedy uczyłem się o modelowaniu wielopoziomowym. Na koniec możesz sprawdzić książkę opublikowaną przez twórców pakietu nlme.
Patrick S. Forscher
Drogi Patricku @ user1188407 dzięki. Studiowałem modele wielopoziomowe i doszedłem do tej formuły, aby przeanalizować moje dane: lme (Velocity_response ~ Velocity * Subject, data = scrd, random = ~ 1 | Subject) Czy możesz mi potwierdzić, że ta formuła jest poprawna do analizy I chcesz wykonać na moich danych? Nie rozumiem jednak, w jaki sposób mogę uzyskać wartości R ^ 2 i p, ani jak wykreślić grafikę z punktami i linią pasującą do regresji. Czy mógłbyś mi pomóc? Nie jestem
statikiem
Formuła wydaje mi się poprawna w oparciu o moje rozumienie twojego badania (i zakładając, że sformatowałeś swoje dane w formacie okresu osobowego). Otrzymałbyś swoje wartości p, zapisując wyniki swojej analizy w obiekcie (tak jak robię to w moich przykładach) i uzyskując streszczenie tego obiektu. Jednak ze względu na różnice między modelami wielopoziomowymi a tradycyjną regresją (np. W metrykach wielkości efektu - nie ma standardowej metryki w modelach wielopoziomowych) zdecydowanie zalecamy przeczytanie więcej o tej technice przed jej użyciem. Wygląda na to, że inni użytkownicy zalecili kilka dobrych opcji.
Patrick S. Forscher,