Spørsmål:
Hvorfor gir lme og aov forskjellige resultater for gjentatte tiltak ANOVA i R?
Henrik
2011-08-10 20:16:04 UTC
view on stackexchange narkive permalink

Jeg prøver å gå fra å bruke ez -pakken til lme for gjentatte tiltak ANOVA (som jeg håper jeg vil kunne bruke egendefinerte kontraster på med lme ).

Etter råd fra dette blogginnlegget klarte jeg å sette opp den samme modellen ved å bruke begge aov (som gjør ez , når det blir bedt om det) og lme . Imidlertid er F -verdiene i eksemplet gitt i det innlegget helt enige mellom aov og lme ( Jeg sjekket det, og de gjør det), dette er ikke tilfelle for dataene mine. Selv om F -verdiene er like, er de ikke de samme.

aov returnerer en f-verdi på 1.3399, lme returnerer 1.36264. Jeg er villig til å godta resultatet aov som det "riktige", ettersom dette også er det SPSS returnerer (og dette er det som teller for mitt felt / veileder).

Spørsmål:

  1. Det ville være flott om noen kunne forklare hvorfor denne forskjellen eksisterer, og hvordan jeg kan bruke lme for å gi troverdige resultater. (Jeg vil også være villig til å bruke lmer i stedet for lme for denne typen ting, hvis det gir "riktig" resultat. Imidlertid har jeg ikke brukt det så langt.)

  2. Etter å ha løst dette problemet vil jeg kjøre en kontrastanalyse. Spesielt ville jeg være interessert i kontrasten med å slå sammen de to første faktornivåene (dvs. c ("MP", "MT") ) og sammenligne dette med det tredje faktornivået (dvs. "AC" ). Videre testing av det tredje versus det fjerde faktornivået (dvs. "AC" versus "DA").

Data: 

  tau.base <- struktur (liste (id = struktur (c (1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L , 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L , 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 1L, 2L, 3L, 4L,
5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L), .Label = c ("A18K", "D21C", "F25E", "G25D "," H05M "," H07A "," H08H "," H25C "," H28E "," H30D "," J10G "," J22J "," K20U "," M09M "," P20E "," P26G ", "P28G", "R03C", "U21S", "W08A", "W15V", "W18R"), klasse = "faktor"), faktor = struktur (c (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, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 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, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L), .Label = c ("MP", "MT", "AC", "DA"), klasse = "faktor"), verdi = c (0,9648092876, 0,2128662077, 1, 0,0607615485, 0,9912814024, 3,22e-08, 0,8073856412, 0,1465590332, 0,9981672618, 1, 1, 1, 0,9794401938, 0,6102546108, 0,428651501, 1, 0,1 710644881, 1, 0.7639763913, 1, 0.5298989196, 1, 1, 0.7162733447, 0.7871177434, 1, 1, 1, 0.8560509327, 0.3096989662, 1, 8.51e-08, 0.3278862311, 0.0953598576, 1, 1.38e-08, , 0.545290432, 0.1305621416, 2.61e-08, 1, 0.9834051136, 0.8044114935, 0.7938839461, 0.9910112678, 2.58e-08, 0.5762677121, 0.4750002288, 1e-08, 0.8584252623, 1, 1, 0.6020383727.5 0.263377904, 1, 1.07e-08, 0.3160751898, 5.8e-08, 0.3460325565, 0.6842217296, 1.01e-08, 0.9438301877, 0.5578367224, 2.18e-08, 1, 0.9161424562, 0.2924856039, 1e-08, 0.8642664992 0.9988463913, 0.2960361777, 0.0285680426, 0.0969063841, 0.6947998266, 0.0138254805, 1, 0.3494775301, 1, 2.61e-08, 1.52e-08, 0.5393467752, 1, 0.9069223275)), .Names = "factor" verdi "), klasse =" data.frame ", rad.names = c (1L, 6L, 10L, 13L, 16L, 17L, 18L, 22L, 23L, 24L, 27L, 29L, 31L, 33L, 42L, 43L, 44L, 45L, 54L, 56L, 58L, 61L, 64L, 69L, 73L, 76L, 79L,
80L, 81L, 85L, 86L, ​​87L, 90L, 92L, 94L, 96L, 105L, 106L, 107L, 108L, 117L, 119L, 121L, 124L, 127L, 132L, 136L, 139L, 142L, 143L, 144L, 148L, 149L, 150L, ​​153L, 155L, 157L, 159L, 168L, 169L, 170L, 171L, 180L, 182L, 184L, 187L, 190L, 195L, 199L, 202L, 205L, 206L, 207L, 211L, 212L, 213L, 216L, 218L, 220L, 222L, 231L, 232L, 233L, 234L, 243L, 245L, 247L, 250L))

Og koden: 

  krever ( nlme) sammendrag (aov (verdi ~ faktor + feil (id / faktor), data = tau.base)) anova (lme (verdi ~ faktor, data = tau.base, tilfeldig = ~ 1 | id))
Det ser ut til at du nettopp har svart på delen om kontrastene selv i svaret ditt [her] (http://stats.stackexchange.com/questions/14078/post-hoc-test-after-anova-with-repeated-measures-using -r / 14142 # 14142); Hvis ikke, kan du redigere dette spørsmålet slik at vi vet hvilken vanskelighet som gjenstår.
@Aaron, så lenge det er forskjeller i `lme`-resultatene fra standard lærebok ANOVA (gitt av` aov`, og det er det jeg trenger), er dette ikke et alternativ for meg. I papiret mitt vil jeg rapportere en ANOVA, ikke noe som en ANOVA. Interessant viser Venables & Ripley (2002, s. 285) at begge tilnærmingene fører til identiske estimater. Men forskjellene i * F * -verdier gir meg en dårlig følelse. Videre returnerer `Anova ()` (fra `bil`) bare Chi²-verdier for` lme`-objekter. Derfor blir ikke mitt første spørsmål besvart ennå.
Jeg forstår (men ikke del) din varighet med `lme '; men for kontraster, fungerer 'glht' på 'lm' passer også, ikke bare 'lme' passer. (Også, `lme`-resultatene er også standard lærebokresultater.)
Dessverre kan du ikke spesifisere 'lm' for en gjentatt måleanalyse. Bare "aov" kan håndtere gjentatte tiltak, men vil returnere et objekt av klasse "aovlist" som dessverre ikke håndteres av "glht".
Er dette bare et mindre eksempel da? Fordi du ikke trenger aov her, da det bare er innen emneeffekter.
@Aaron kan du utdype mer om dette. Hvorfor ingen `aov`? Generelt er jeg interessert i en generell løsning med både mellom og innen faktorer. Men hvis det er en annen fin måte for deler av dette problemet, vil jeg også være interessert.
`lm` bruker restfeilen som feiluttrykk for alle effekter; når det er effekter som skal bruke et annet feiluttrykk, er 'aov' nødvendig (eller i stedet for å bruke resultatene fra 'lm' for å beregne F-statistikken manuelt). I ditt eksempel er feiluttrykket for "faktor" interaksjonen "id: faktor", som er restfeiluttrykket i en additivmodell. Sammenlign resultatene dine med `anova (lm (verdi ~ faktor + id))`.
Jeg har svart på det første spørsmålet ditt nedenfor. Siden svaret er ganske langt, vil jeg foreslå at du starter et nytt spørsmål for det andre spørsmålet ditt.
Hvorfor skriver du "tilfeldig = ~ 1 | id" og ikke "tilfeldig = ~ 1 | (id / faktor)"?
To svar:
Aaron left Stack Overflow
2011-08-12 20:58:33 UTC
view on stackexchange narkive permalink

De er forskjellige fordi lme-modellen tvinger varianskomponenten til id til å være større enn null. Ser vi på den rå anovatabellen for alle termer, ser vi at den gjennomsnittlige kvadrerte feilen for id er mindre enn den for restene. 

  > anova (lm1 <- lm (verdi ~ faktor + id, data = tau.base)) Df Sum Sq Gjennomsnitt Sq F verdi Pr (>F) faktor 3 0.6484 0.21614 1.3399 0.2694id 21 3.1609 0.15052 0.9331 0.5526 Rester 63 10.1628 0.16131

Når vi beregner varianskomponentene, betyr dette at variansen på grunn av id vil være negativ. Minneet mitt om forventet gjennomsnittlig kvadratminne er skjelvende, men beregningen er omtrent som 

  (0.15052-0.16131) / 3 = -0.003597.

Dette høres ut rart, men kan skje. Hva det betyr er at gjennomsnittene for hver id er nærmere hverandre enn du forventer hverandre gitt mengden gjenværende variasjon i modellen.

Hvis du bruker lme, tvinger denne variansen til å være større enn null. 

  >-sammendrag (lme1 <- lme (verdi ~ faktor, data = tau.base, tilfeldig = ~ 1 | id)) ... Tilfeldige effekter: Formel: ~ 1 | id (Intercept) ResidualStdDev: 3.09076e-05 0.3982667

Dette rapporterer standardavvik, kvadrat for å få variansutbyttet 9.553e-10 for id-variansen og 0.1586164 for den gjenværende variansen.

Nå bør du vite at bruk av aov for gjentatte tiltak bare er passende hvis du mener at korrelasjonen mellom alle par av gjentatte tiltak er identiske; dette kalles sammensatt symmetri. (Teknisk sett er det nødvendig med sfærisitet, men dette er tilstrekkelig for nå.) En grunn til å bruke lme over aov er at den kan håndtere forskjellige typer korrelasjonsstrukturer.

I dette bestemte datasettet er estimatet for denne korrelasjonen negativt; Dette hjelper til med å forklare hvordan den gjennomsnittlige kvadrerte feilen for id var mindre enn den gjenværende kvadratiske feilen. En negativ korrelasjon betyr at hvis et individs første måling var under gjennomsnittet, i gjennomsnitt, ville deres andre være over gjennomsnittet, noe som gjør det totale gjennomsnittet for individene mindre variabelt enn vi ville forvente hvis det var en null korrelasjon eller en positiv korrelasjon. p>

Bruk av lme med en tilfeldig effekt tilsvarer å tilpasse en sammensatt symmetri-modell der korrelasjonen blir tvunget til å være ikke-negativ; vi kan passe på en modell der korrelasjonen får være negativ ved hjelp av gls:

  > anova (gls1 <- gls (verdi ~ faktor, korrelasjon = corCompSymm (form = ~ 1 | id), data = tau.base)) Denom. DF: 84 numDF F-verdi p-verdi (Intercept) 1 199.55223 <.0001factor 3 1.33985 0.267

Denne ANOVA-tabellen stemmer overens med tabellen fra aov fit og fra lm fit.

OK, så hva? Vel, hvis du mener at avviket fra id og korrelasjonen mellom observasjoner skal være ikke-negativ, er lme passformen faktisk mer passende enn passformen ved hjelp av aov eller lm da estimatet av restvariansen er litt bedre. Men hvis du mener at sammenhengen mellom observasjoner kan være negativ, er aov eller lm eller gls bedre.

Du kan også være interessert i å utforske korrelasjonsstrukturen nærmere; for å se på en generell korrelasjonsstruktur, vil du gjøre noe sånt som 

  gls2 <- gls (verdi ~ faktor, korrelasjon = corSymm (form = ~ uklass (faktor) | id), data = tau.base)

Her begrenser jeg bare produksjonen til korrelasjonsstrukturen. Verdiene 1 til 4 representerer de fire nivåene av faktor ; vi ser at faktor 1 og faktor 4 har en ganske sterk negativ korrelasjon: 

  >-sammendrag (gls2) ... Korrelasjonsstruktur: Generell formel: ~ uklass (faktor) | id Parameterestimat (er): Korrelasjon: 1 2 3 2 0,049 3 -0,127 0,208 4 -0,400 0,146 -0,024

En måte å velge mellom disse modellene er med en sannsynlighetsforholdstest; Dette viser at modellen for tilfeldige effekter og den generelle korrelasjonsstrukturmodellen ikke er statistisk signifikant forskjellige; når det skjer foretrekkes vanligvis den enklere modellen. 

  > anova (lme1, gls2) Model df AIC BIC logLik Test L.Ratio p-valuelme1 1 6 108.0794 122.6643 -48.03972 gls2 2 11 111.9787 138.7177 -44.98936 1 mot 2 6.100725 0.2965
Flott svar. Takk. En rask oppfølging. Er det en måte å tillate negative avvik i `lme`?
Ikke bruker `lme`, nei. Men `gls` med sammensatt symmetri er ekvivalent, som jeg viste ovenfor.
Det er faktisk mulig å bruke sammensatt symmetri med 'lme' for å oppnå de samme resultatene som med 'aov' (og derved muliggjøre 'lme' for alle ANOVAer), nemlig å bruke korrelasjonsargumentet i kallet til 'lme': 'anova ( lme (verdi ~ faktor, data = tau.base, tilfeldig = ~ 1 | id, korrelasjon = corCompSymm (form = ~ 1 | id)) `
Hyggelig funn. Men er det ikke en ekstra parameter som passer? Den har tre variansparametere; varians for id, restvarians og korrelasjon, mens gls bare har en restvarians og en korrelasjon.
Argumentet ditt høres plausibelt ut, men resultatene stemmer ikke overens. Alle anovatabellene (`aov`,` lme` uten sammensatt symmetri og `lme med sammensatt symmetri) har nøyaktig samme antall dfs.
Du må overbevise meg om at de tre parametrene virkelig er en overparametrisering av de to første. Har du funnet ut hvordan de er i slekt?
Nei. Jeg stoler på utdataene fra `anova.lme ()`. Fra svaret ditt fikk jeg at forholdet mellom ANOVA og de blandede modellene ligger i deres korrelasjonsstruktur. Jeg leste da at å pålegge en sammensatt symmetrisk korrelasjonsstruktur fører til likhet mellom de to tilnærmingene. Derfor innførte jeg det. Jeg aner ikke om dette spiser opp en annen df. Produksjonen er imidlertid uenig med denne tolkningen.
Vel, du er mer tillitsfull enn meg da. Jeg vet at Doug Bates ikke var fornøyd med df-beregningen (derav endringene i lme4), så jeg ville være forsiktig med mindre jeg visste nøyaktig hva som foregikk.
Hvorfor "tilfeldig = ~ 1 | id" og ikke "tilfeldig = ~ faktor | id"?
Fra hukommelsen ville det gi forskjellige avvik for hvert nivå av faktoren, og jeg tror også forskjellige kovarianter for hvert par.Det kan være verdt å utforske også, avhengig av detaljene i situasjonen din.En potensiell ulempe er at det kan være mange parametere å passe, spesielt hvis det er mange nivåer av faktoren.
Hva er betydningen av en negativ varians?Så vidt jeg vet er variansen summen av kvadrater.
Godt spørsmål!Det er mer enn det kan besvares i en kommentar.Dette svaret https://stats.stackexchange.com/a/253323/3601 antyder at det er mange spørsmål om det, men i min søken fant jeg ikke en som enten var bra eller med et godt klart kanonisk svar.Hvis du er interessert i å dra nytte av nettstedet, ville det være flott å gjøre et søk for å bekrefte det og deretter skrive et nytt klart spørsmål om det, med henvisning til de andre, og håpe at din blir det kanoniske spørsmålet om det.
Chris
2011-08-10 21:40:05 UTC
view on stackexchange narkive permalink

aov () passer til modellen via lm () med minst kvadrat, lme passer med maksimal sannsynlighet. Den forskjellen i hvordan parametrene til den lineære modellen er estimert, utgjør sannsynligvis (veldig liten) forskjell i f-verdiene.

I praksis (f.eks. for hypotesetesting) er disse estimatene de samme, så jeg kan ikke se hvordan den ene kan betraktes som 'mer troverdig' enn den andre. De kommer fra forskjellige modelltilpassede paradigmer.

For kontraster, må du sette opp en kontrastmatrise for faktorene dine. Venebles og Ripley viser hvordan du gjør dette på s 143, s.146 og s.293-294 i 4. utgave.

Hmm, men hvorfor er det noen ganger forskjeller, og noen ganger er resultatene nøyaktig like? Fortsatt ser det ut til å være umulig å bruke `lme` eller` lmer` for å beregne en ANOVA (strengt tatt), da den bruker en metode som er lik, men ikke identisk. Så er det ingen måte å beregne kontraster for ANOVAer med gjentatt mål i R?
Hvis systemet din modellering er virkelig lineær enn minste kvadrater, og ML skal gi den samme statistikken. Det er bare når det er annen struktur i dataene at de to metodene vil gi forskjellige resultater. Pinheiro og Bates dekker dette i modellboken med blandede effekter. De er sannsynligvis ikke "nøyaktig" like, hvis du skulle gå langt nok i sifre, er jeg sikker på at du vil finne en forskjell. Men for alle praktiske formål er de de samme.


ⓘ
Denne spørsmålet ble automatisk oversatt fra engelsk.Det opprinnelige innholdet er tilgjengelig på stackexchange, som vi takker for cc by-sa 3.0-lisensen den distribueres under.
Loading...