12 Korrelerede målinger

12.1 Omstrukturering af data

Data kan optræde i langt format, hvor hver enkelt måling på et individ er givet i en række for sig, eller i bredt format, hvor hvert individ har én række med hver måling angivet i en søjle.

I analysen af korrelerede målinger skal vi oftest have data på langt format. Dvs hvert individ skal have lige så mange linjer, som det har målinger af outcome. Ofte får vi dog data på bredt format, men det er nemt nok at lave om til langt format:

12.1.1 Bredt til langt format

Eksempel: Kanin-data uge 10

Er data i bredt format, og vi ønsker at formatere det til langt format, kan det gøres med melt()-funktionen fra data.table-pakken. Først angiver vi hvilket data, der skal transformeres. Datasættet skal angives som en data.table, hvilket fås ved at bruge as.data.table. Hernæst angives i id.vars den variabel, som angiver individ-niveauet, og i measure angiver vi hvilke variable, der skal samles til en enkelt variabel. I dette tilfælde er det variablene med navnene a, b, c, d, e og f. Med variable.name kan vi navngive den nye variabel, der beskriver hvilken søjle observationen står i i det oprindelige data, og med value.name kan vi angive navnet på den nye variabel, der indeholder værdien af observationen i det oprindelige data.

library(data.table)

lang <- melt(as.data.table(kanin), 
             id.vars = c("kanin"), # Id-variabel
             measure = c("a","b","c","d","e","f"), # Navne på de variable som indeholder målinger
             variable.name = "sted", # Navn på den variabel som indeholder sted
             value.name = "swelling") # Navn på vores outcome variabel
head(lang, 10)
    kanin sted swelling
 1:     1    a      7.9
 2:     2    a      8.7
 3:     3    a      7.4
 4:     4    a      7.4
 5:     5    a      7.1
 6:     6    a      8.2
 7:     1    b      6.1
 8:     2    b      8.2
 9:     3    b      7.7
10:     4    b      7.1

12.2 Plot af individuelle profiler

For at illustrere observationer over tid for hvert individ, kan man lave et såkaldt spagetti-plot. Dette kan laves ved brug af ggplot.

Vi betragter her eksemplet med calciumtilskud uge 10.

I ggplot() angiver vi datasættet (calcium) og inden for aes() angiver vi, hvad der skal være på x-aksen og y-aksen i figuren, samt hvordan vi ønsker data grupperet (her efter girl, da det angiver individniveauet) samt at vi ønsker forskellig farve, afhængigt af om pigen har alle målinger (dropout=0) eller ej (dropout=1). Her skal vi definere dropout som en faktor for at undgå at få farvekoden vist som en farvegradient. Da jeg gerne vil have completers i blå (dropout=0, dvs skal være farvekode nr 2 og ikke farvekode nr 1 som ellers er default fordi variabelkode 0 kommer før 1), dropouts i rød (dropout=1, skal være farvekode nr 1 og ikke nr 2), bytter jeg rundt på kode 0 og 1 med dropout-1. Hernæst tilføjer vi datapunkter med + geom_point() og forbinder observationerne med linjer med + geom_line().

library(ggplot2)

ggplot(calcium, aes( x = visit, y = bmd, group = girl, color=as.factor(1-dropout)) ) + 
  geom_line() + geom_point()

Ønsker vi at adskille observationerne for pigerne, som får placebo fra pigerne, som får calcium, kan vi nemt konstruere to separate plots for de to grupper ved at tilføje +facet_wrap(~grp), hvor grp angiver placebo/calcium. Samtidigt justeres placeringen af legend, så den optræder under figuren vha. theme(legend.position="bottom"):

ggplot(calcium, aes( x = visit, y = bmd, group = girl, color=as.factor(1-dropout))) + 
  geom_line() + geom_point() + facet_wrap(~grp) +
  theme(legend.position="bottom")

12.3 Varianskomponentmodel

12.3.1 Estimation

Vi kan estimere en varianskomponentmodel (mixed model) vha. lme()-funktionen fra nmle-pakken. I lme() kan vi angive både de systematiske (fixed) effekter og de tilfældige (random) effekter.

De systematiske effekter angives via en formel som vi plejer: \[ \text{respons} \sim \text{variabel1}+ \text{variabel2} + ...\] Interaktioner mellem variabel kan angives med *.

De tilfældige effekter, der beskriver korrelationen mellem observationerne, angives i et argument random. Efter =~ angiver vi hvilke variable der skal benyttes til at beskrive afhængighedsstrukturen, og efter | angives den variabel, som angiver hvilken person, målingerne er taget på (patientens idnummer). Den samlede form af argumentet er derfor random =~ variabel1 | idnummer. Sættes 1 som variabel1 svarer det til at have en tilfældig effekt på et intercept hørende til hver patient, er variabel1 kvantitativ svarer det til tilfældig effekt på både intercept og hældning (dvs. en random regression).

I eksemplet her har vi angivet en model, hvor middelværdien af bmd beskrives ved en interaktion mellem grp og visit (systematiske effekter), og med et tilfældigt intercept for hver pige (random intercept):

# install.packages("nlme")
library(nlme)

calcium$girl = as.factor(calcium$girl)
calcium$visit = as.factor(calcium$visit)

model.var = lme(bmd ~ grp*visit , data = calcium, na.action = na.exclude, 
                random =~ 1 | girl)
  • Med fixef() får vi koefficienterne hørende til de systematiske effekter
  • Med summary() får vi de systematiske estimater (dvs køfficienter, helt som vi plejer). Vi får også en korrelationsmatrix frem, som angiver korrelationen mellem de estimerede parametre. Den skal vi ikke bruge til noget
  • Med intervals() får vi parameterestimater med tilhørende konfidensintervaller for både de systematiske effekter og spredningen på de tilfældige effekter.
  • Variansen på de tilfældige effekter kan fås med VarCorr (her får vi både varians og spredning (=sqrt(varians)), med intervals() får spredningen incl et CI).
fixef(model.var)
summary(model.var)
intervals(model.var)
VarCorr(model.var)

Ønsker vi at teste for interaktion, kan vi benytte anova():

anova(model.var)
            numDF denDF   F-value p-value
(Intercept)     1   381 20974.226  <.0001
grp             1   110     2.217  0.1394
visit           4   381   616.848  <.0001
grp:visit       4   381     5.297  0.0004

HUSK at vi her kun må aflæse p-værdien i nederste linje. drop1() som vi ellers er så glade for, virker ikke her. Ønsker vi at sammenligne to modeller, hvor den ene er en grovere model end den anden mht de variable vi har med som systematiske effekter, kan vi gøre det ved at fitte begge modeller med ekstra argument method=ML i lme() og derefter benytte anova(model1,model2). Skal vi rapportere estimaterne fra en af disse to modeller, skal vi dog fitte den igen uden method=ML-argumentet.

12.3.2 Modeller med fælles baseline

Se script hørende til uge 10, hvordan man tvinger middelværdien i to grupper til at være ens ved baseline

12.4 Modelkontrol

Modelkontrol for en varianskomponentmodel udføres ved at tjekke fordelingerne af de

  • de betingede residualer (små residualer, observeret minus prædikteret værdi på individniveau, dvs både systematiske og tilfældige effekter trækkes fra).
  • de sædvanlige residualer (store residualer, observeret minus prædikteret middelværdi, dvs kun systematiske effekter trækkes fra).

12.4.1 Små residualer

Begge typer residualer fås ved hjælp af resid()-funktionen. Her fås pr. automatik de betingede / små residualer.

Vi kontrollerer normalfordeling af disse:

par(mfrow=c(1,2))
hist( resid(model.var) )
qqnorm( resid(model.var) )
qqline( resid(model.var) )

Plotter mod de prædikterede værdier at undersøge varianshomogenitet:

scatter.smooth( predict(model.var), resid(model.var),
                lpars=list(col="red",lwd=2))
abline(h=0, lwd=.5)

12.4.2 Store residualer

De store residualer får man frem med ekstra argument level=0 i resid(). Tilsvarende kan vi få prædiktionerne frem på populations/gennemsnits-niveau med predict(, level=0)

scatter.smooth( predict(model.var, level=0), resid(model.var, level=0),
                lpars=list(col="red",lwd=2))

Her ligger plottet i pølser - der er en mulig prædikteret gennemsnitsværdi for hver kombination af visit og grp (5*2=10)

Når vi skal checke normalfordelingsantagelsen, skal man sno sig lidt i R. Nu bliver det teknisk … De store (sædvanlige) residualer er ikke uafhængige. Hvis korrelationen mellem målingerne for hver pige er meget høj vil de sædvanlige residualer være meget korrelerede. Hvis der samtidigt er stor forskel på, hvor højt pigerne ligger i niveau betyder det, at de sædvanlige residualer ligner hinanden meget for hver enkelt pige og at disse residualer så kommer til at ligge i ‘klumper’. Vi kan i så fald ikke forvente at få en normalfordeling frem i qq-plot, når vi lægger disse ‘klumper’ sammen. I dette calcium-eksempel er det ikke et problem fordi inter- og intraindivid variationen ikke er så høj. Men i øvelsesopgaven i denne uge får I et QQ-plot frem, som ser ret sjovt ud, hvis I bare regner på de store residualer som her:

qqnorm( resid(model.var, level=0), main="FORKERT QQ-plot")
qqline( resid(model.var, level=0))

I stedet skal man skalere residualerne på en måde, så de bliver uafhængige og vi derved kan vurdere normalfordeling. Koden kan I ikke forstå.

#install.packages('mgcv')
library(mgcv) 
res <- residuals(model.var, level=0, na.action="na.exclude")
cov <- extract.lme.cov(model.var, calcium, ) 
L <- t(chol(cov))  
scaled.residuals <- solve(L) %*% na.omit(res)
qqnorm(scaled.residuals)
qqline(scaled.residuals)

Vi ser, at QQ-plottet faktisk bliver lidt pænere, når vi laver det rigtigt.

Hver gang du skal checke normalfordeling af de store residualer, skal du benytte koden ovenfor - lav copy and paste, erstat model.var med navnet på den model du vil kontrollere.

12.5 Korrelationsstrukturer

Når vi har en korrelationsstruktur, vi ikke modellere med random effects, skal vi have fat gls()-funktionen (Generalized Least Squares’). Her angiver vi igen de systematisk effekter med formelnotation, mens korrelationsstrukturen skal specificeres i et ekstra argument corr.

12.5.1 Ustruktureret (UN)

Korrelations strukturen angives i argumenterne corr= og weight. For at opnå ustruktureret korrelation sættes corr = corSymm(form=~1 | girl) og weight = varIdent(form =~ 1 | visit).

model.un = gls(bmd ~ grp*visit, data = calcium, na.action = na.exclude, 
               corr = corSymm(form=~1 | girl), weight = varIdent(form =~ 1 | visit), method = "REML")

Parameterestimater for modellen kan fås vha. summary().

12.5.2 Autoregressiv korrelation

En model med autoregressiv korrelationsstruktur kan i R estimeres med gls-funktionen, hvor vi angiver correlationsstrukturen corr = corAR1(). Er tidsmålingerne ækvidistante kan vi sætte corr=corAR1(form =~ 1 |girl). Er de ikke ækvidistante sætter vi corr=corCAR1(form=~visit|girl)).

model.ar = gls(bmd ~ grp*visit, data = calcium, na.action = na.exclude, 
               corr = corAR1(form =~ 1 |girl), method = "REML")

Igen kan vi udtrække parameter estimater med summary().

12.6 Random regression

Random regression kan i R laves ved enten lme() som ovenfor eller vha. lmer-funktionen fra lme4-pakken. Sidstnævnte pakke er nødvendig, hvis man har lidt mere komplicerede kovariansstrukturer som kræver flere tilfældige effekter (det har I ikke i dette kursus).

12.6.1 Via lme()

Her tilføjer vi den variabel, som skal angive random slope i random-linjen:

rr <- lme(bmd ~ grp*years, random=~years| girl, 
          data=calcium, na.action = na.exclude)

Derefter kan vi benytte funktioner som beskrevet ovenfor til at hive estimater, residualer mm ud.

12.6.2 Via `lmer()

Her kan vi angive de systematiske effekter i formelnotation som vi plejer, og de tilfældige effekter tilføjes i parentes i selve formellinjen. Det er her mulig at tilføje flere tilfældige effekter i samme model.

I eksemplet nedenfor har vi som faste effekter en interaktion mellem grp og years. Derudover tillader vi som en tilfældig effekt at hældingerne over tid (years) er forskellig for hvert individ (girl). Dette gøres ved at skrive (years | girl) (hvorved dett tilfældige intercept automatisk kommer med).

#install.packages("lme4")
library(lme4)
rr <- lmer(bmd ~ grp*years + (years | girl), data=calcium)

Ønsker vi kun at have en tilfældig effekt på interceptet for hvert individ, kan vi i stedet skrive (1|girl).

Vi kan få parameteresimater og konfidensintervaller med hhv. summary() og confint().

summary(rr)
confint(rr)

Prøv at køre koden - og bemærk dog at forfatteren til denne pakke har sørget for at vi ikke får p-værdierne med, fordi han er modstander af p-værdier… Så man enten bestemme dem selv (nemmest bare at lave Wald-test, han angiver trods alt teststørrelserne), eller køre modellen via lme() hvis muligt.