8 Lineære modeller

I litteraturen skelner man i analyse af et kvantitativt outcome (f.eks. vitamin D) mellem

  • ANOVA - de forklarende variable er kategoriske (kvalitative/faktorer)
    • one-way én forklarende variabel
    • two-way to forklarende variable
    • multi-way 3+ forklarende variable
  • regression - de forklarende variable er kvantitative
    • simpel (emne uge 4)
    • multipel
  • ANCOVA - én kvantitativ og én kategorisk (emne uge 5)
  • generel lineær model / multipel regression / multivariabel lineær model - blanding af kvantitative og kvalitative (emne uge 6)

I bund og grund er de alle specialtilfælde af den samme klasse af modeller - den Generelle Lineære Model (GLM). I SAS kan vi benytte proc glm til alle disse analyser. I uge 6 introduceres I til en anden procedure, proc reg, som kun kan håndtere kvantitative forklarende variable, men til gengæld kan foretage automatisk model selektion og har lidt flere muligheder mht diagnostics.

Generel syntaks proc glm:

  • I proc glm-linjen tilføjes udover et data argument også plots=all for at få modelkontrolplots mm ud
  • I class-linjen fortæller man, hvilke variable SAS skal opfatte som kategoriske (class/faktorer)
  • I model-linjen
    • sættes responsen på venstre side af lighedstegnet
    • sættes en eller flere forklarende variable på højre side af lighedstegnet
    • efter en / skrives ALTID clparm - ellers får man ikke konfidensintervaller på parametrene frem
    • i ANOVA (dvs når man kun har kvalitative/class-variable i sin model) skriver man desuden solution efter /, ellers får man ikke engang parametrene med i output.

Syntaks med et tænkt datasæt med navn datasaet, to forklarende variable (forklarende1, forklarende2) hvor forklarende1 er kategorisk mens forklarende2 er kvantitativ:

proc glm data=datasaet;
   class forklarende1;
   model respons = forklarende1 forklarende2 / solution clparm;
run;

Generel syntaks proc reg:

Meget lig syntaksen for proc glm, men man kan dog ikke benytte en class-linje, man skriver clb i stedet for clparm for at få CI.

Syntaks med et tænkt datasæt med navn datasaet, to forklarende variable (forklarende1, forklarende2) som begge er kvantitative (svarende til multipel regression, uge 6):

proc reg data=datasaet;
   model respons = forklarende1 forklarende2 / clb;
run;

8.1 One-way ANOVA

Eksempel: Vitamin D (uge 2)

Vitamin D (vitd) som respons, land (country) som forklarende variabel:

proc glm data=vit;
    class country;
    model vitd=country / solution clparm;
run;

Her får vi et overordnet test for association (ANOVA, p-værdien aflæses under Type III test) og estimater for de individuelle forskelle med højeste gruppenummer som reference.

8.1.1 Test for ens varianser

Vi kan bede proc glm lave test for om varianserne i grupperne er ens ved at tilføje en means-linje (her skal vi specificere navnet på den variabel som angiver den inddeling, for hvilken vi ønsker at sammenligne varianserne, dvs her country):

proc glm data=vit;
    class country;
    model vitd=country / solution clparm;
    means country / hovtest welch;
run;

Her giver hovtest (Homogeneity Of Variance test) Levene’s test for ens varianser, mens welch giver os det overordnede (Welch-)test for association uden antagelse om at variansen i grupperne(=landene) er ens.

8.1.2 Korrektion for multiple sammenligninger

Her tilføjes en lsmeans-linje, hvor den variabel vi ønsker at justere for multiple sammenligninger for angives. I adjust angives den metode, man ønsker at benytte (f.eks tukey, sidak eller dunnett mfl, se evt SAS dokumentation). Ekstra options pdiff og cl beder om p-værdier og konfidensintervaller hørende til de multiple sammenligninger.

proc glm data=vit;
    class country;
    model vitd=country / solution clparm;
    lsmeans country / adjust=tukey pdiff cl;
run;

8.2 Two-way ANOVA

Eksempel: Vitamin D (uge 2)

En lidt anden model end gennemgået i slides her. I en two-way ANOVA har vi to forklarende kvalitative variable. Vi kan her vælge land (country) og solvaner (sunexp). De sættes begge i class-linjen og på højre side af lighedstegnet i model-linjen.

proc glm data=vit;
   class country sunexp;
     model vitd=country sunexp / solution clparm;
run;

Det er muligt at tilføje en means såvel som en lsmeans-linje (se ovenfor, kap 8.1)

Interaktion:

  • til at teste om der er interaktion angives i model-linjen: country sunexp country*sunexp
  • til at estimere effekten af country i en model med interaktion angives i model-linjen: sunexp country*sunexp
  • til at estimere effekten af sunexp i en model med interaktion angives i model-linjen: country country*sunexp

8.3 Simpel lineær regression

Eksempel: Vietnam veteraner (IQ og lungegfunktion, uge 2)

FEV1 som outcome, iq som forklarende kvantitativ variabel:

proc glm data=viet;
    model FEV1 = iq / clparm;
run;

Analysen kan også laves vha proc reg (her skriver man dog clb i stedet for clparm for at få CI i output):

proc reg data=viet;
    model FEV1 = iq / clb;
run;

8.4 Multipel lineær regression

Eksempel: Lungefunktion (uge 6)

Hvis de forklarende variable alle er kvantitative kan vi benytte proc reg (men får det samme frem med proc glm):

proc reg data=pemax;
    model pemax=age sex height weight bmp fev1 rv frc tlc;
run;

proc glm data=pemax;
    model pemax=age sex height weight bmp fev1 rv frc tlc;
run;

8.4.1 Modelselektion

Vi fraråder i det hele taget at benytte variabelselektionsprocedurer

Når vi har mange mulige kovariater i spil og gerne vil have reduceret dette antal uden at have en holdning til, hvilke variable der bør være med i modellen, vælger nogle at lave automatisk variabelselektion. Dette kan gøres på flere måder:

  • Backwards elimination: Fjerner i hvert skridt den mindst signifikante
  • Forwards selection: Tilføjer i hvert skridt den mest signifikante
  • Hybrid mellem forwards og backwards

Automatiseret modelselektion (forlæns og baglæns) kan man udføre i proc reg (dette er ikke en mulighed i proc glm). Man tilføjer blot et selection=forward eller selection=backward i modellinjen efter en slash. Uanset om man søger forlæns eller baglæns, skal man have alle variable i spil med i modellinjen. For baglæns søgning:

proc reg data=pemax;
    model pemax=age sex height weight bmp fev1 rv frc tlc
    / selection=backward;
run;

8.5 Modelkontrol

Eksempel: Vietnam veteraner (IQ og lungefunktion, uge 2)

Med plots=all, plots=diagnostics eller plots=DiagnosticsPanel i proc glm-linjen får vi automatisk QQ-plot og histogram (til vurdering af normalfordeling) og residualplot (til vurdering af varianshomogenitet), f.eks.

proc glm data=viet plots=DiagnosticsPanel;
    model  FEV1 = iq / clparm;
run;

Vil vi også have plot af residualerne mod en kvantitativ kovariat tilføjer vi yderligere residuals(smooth) så syntaksen lyder (bemærk parantes til at samle plot-options):

proc glm data=viet plots=(diagnostics residuals(smooth));
    model FEV1 = iq / clparm;
run;

Vi får samme plots frem uanset om vi laver ANOVA (med class-variable) eller har flere forklarende variable med i modellen.

Ovenstående virker også for proc reg.

8.5.1 Residualer mm

Eksempel: Vitamin D (uge 1 og 5)

Først laver vi en id-variabel så det er nemt at identificere individerne, når vi laver print af data. Derudover laver vi en ny variabel sol i to kategorier

proc format;
     value solf 0 = "Nej" 1 = "Ja";
run;

data vit; 
     set vit; 
     id=_N_;

     if sunexp=1 then sol=0;
     if sunexp in (2,3) then sol=1;
     format sol solf.;
run;

Vi benytter residualer til modelkontrol. Nogle gange ønsker vi at arbejde videre med disse, eller benytte størrelser, som ikke er benyttet i SAS’ automatiske plots. Man kan beregne et hav af størrelser for hver observation og gemme dem i et datasæt ved at benytte et output-statement, f.eks.

proc glm data=vit;
     model vitd = sol bmi / solution clparm;
     output out=ny p=yhat r=resid student=stresid
        rstudent=uresid press=press
        cookd=cook;
run;

Med out= navngiver vi datasættet (her er valgt ny) og med de øvrige argumenter specificerer vi, hvad vi vil have beregnet i datasættet. Hvad der kan stå på venstre side af hvert lighedstegn er bestemt af SAS, feks står p for Prediction. På højre side angiver vi, hvad vi ønsker variablen skal hedde i vores nye datasæt. Dvs når vi skriver p=yhat får vi en variabel med navn yhat som indeholder prædiktionen for hver observation. En komplet liste af muligheder finder du på SAS’ hjemmeside her.

I eksemplet ovenfor får vi et nyt datasæt med navn ny (som vi evt kan tage et kig på i ‘Output data’-fanen eller ved proc print data=ny;). Dette datasæt indeholder:

  • en kopi af de oprindelige variable (her vitd, sol og bmi)
  • yhat som er den prædikterede værdi (dvs a + b x iq + c x sol), a er intercept, b er hældning, c forskel i intercept)
  • resid som er residualet, dvs FEV1-yhat
  • stresid som er Studentized residual, dvs residualet divideret med dens spredning
  • rstudent er Studentized residual, men baseret på en beregning hvor den pågældende observation er udeladt
  • press er residualet FEV1-yhat men hvor prædiktionen yhat er baseret på en beregning hvor den pågældende observation er udeladt
  • cookd Cook’s distance

Hvordan vi bruger nogle af disse størrelser fremgår af nedenstående afsnit.

8.5.2 Ekstra check af varianshomogenitet

Eksempel: Vietnam veteraner (IQ og lungefunktion, uge 2)

Nogle gange kan det være nemmere at vurdere varianshomogenitet ud fra et plot af kvadratroden af den absolutte værdi af residualerne mod de prædikterede værdier. Dette plot skal vi selv lave, da det ikke er en del af de i SAS implementerede plots

Vi skal selv beregne kvadratroden af residualerne, inden vi plotter dem mod de prædikterede værdier. Vi benytter datasættet ny defineret i kap 8.5.1, hvor vi gemte residualerne i variablen resid og laver en ny variabel, sqrt_abs_res som er kvadratroden af den numeriske værdi af residualet:

data ny;
    set ny;
    sqrt_abs_res=sqrt(abs(resid));
run;

Vi plotter nu disse residualer mod prædikterede værdier og lægger en udglatning oveni:

proc sgplot data=ny;
    loess x=yhat y=sqrt_abs_res;
run;

8.5.3 Linearitet

Eksempel: Vietnam veteraner (IQ og lungefunktion, uge 2)

Kun relevant for lineær regression. Hertil plottes residualerne mod kovariaten. Dette kan tilføjes ved plots=residuals(smooth) (med (smooth) angiver vi, at vi ønsker en ‘smoother’ lagt ovenpå (en slags glidende gennemsnit, også kaldt en loess-kurve)). Ønsker man at angive flere argumenter til plots, fordi vi ønsker flere slags plot, kan vi angive disse samlet i ():

proc glm data=viet plots=(diagnostics residuals(smooth));
    model FEV1 = iq / clparm;
run;

Vi kan også bare lave plottet selv ved at plotte residualerne mod kovariaten i vores datasæt ny lavet vha output-statement (kap 8.5.1):

proc sgplot data=ny;
    loess x=iq y=resid;
run;

8.6 Estimate / kvantificering

Eksempel: AMH og P-piller (uge 4)

I P-pille eksemplet så vi i forelæsningerne på en model med interaktion mellem alder og P-piller, logaritmetransformeret AMH som outcome. Vi kan regne videre på koefficienterne (=estimaterne), når vi har fittet en model. Her kan vi benytte estimate-linjer. I videoen nedenfor (20 min) gennemgår jeg estimaterne og viser hvordan man kan regne videre på disse i modellen med interaktion:

proc glm data=pp (plots=all);
   class ppiller;
   model logamh = ppiller alder ppiller*alder / solution clparm;
run;

Videoen kan du finde her, hvis den ikke vises automatisk herunder.

Videoens indhold:

Jeg finder frem til linjens ligning for hver af de to P-pille-grupper:

  • Ikke-brugere: 2.32 - 0.03*alder
  • Brugere: 2.50 - 0.04*alder

og viser hvordan man bestemmer CI for hældningen for ikke-P-pille brugere med estimate:

proc glm data=pp;
   class ppiller;
   model logamh = ppiller alder ppiller*alder / solution clparm;
   estimate "hældning ikke-brugere"" alder 1 ppiller*alder 1 0;
run;

Forskellen mellem brugere og ikke-brugere for fastholdt alder finder man ved brug af estimate i følgende form:

proc glm data=pp;
   class ppiller;
   model logamh = ppiller alder ppiller*alder / solution clparm;
   estimate "dif 25 år" ppiller -1 1 ppiller*alder -25 25;
   estimate "dif 40 år" ppiller -1 1 ppiller*alder -40 40;
run;

8.7 Diagnostics

8.7.1 Cook’s afstand

Eksempel: Vitamin D (uge 5)

Vi kan plotte Cook’s afstand mod f.eks. observationsnummeret for at se, hvem der evt har en stor afstand. Id-variabel og Cook beregnede vi i datasættet ny ovenfor (kap 8.5.1) og kan tage et print af observationerne med stor Cook:

proc print data=ny;
    var id cook;
    where cook > .0188; * =4/n=4/213;
run;

Vi kan let få et plot frem vha proc reg (bemærk at vi her ikke kan benytte kvalitative (class)-variable, så disse skal være kodet 0/1 (og har man flere niveauer af en kvalitativ variabel må man lave en 0/1 variabel for hvert niveau (dummy-variable)):

proc reg data=vit plots(label)=(CooksD);
    id id;
    model vitd=sol bmi;
run;

8.7.2 Influence

Eksempel: Vitamin D (uge 5)

Influence er Cook’s afstand spaltet ud på de enkelte estimerede parametre (dvs intercept og hældninger). Dette laves vha proc reg og et influence i modellinjen:

proc reg data=vit plots(label)=(dfbetas);
    model vitd = sol bmi / influence;
    ods output OutputStatistics=OutputStatistics;
run;
proc print data=OutputStatistics (obs=10);
run;

Influence kan kun fås frem ved brug af proc reg, dvs her dur det ikke med proc glm.

8.7.3 Kollinearitet (VIF)

Eksempel: Lungefunktion (uge 6)

I en multipel regressionsmodel er der en risiko for, at der er kollinearitet mellem nogle af kovariaterne.

For at undersøge kollinearitet benyttes Variace Inflation Factor (VIF), Tolerance (1/VIF) og standardiserede koefficienter som fås frem med ekstra argumenter i model-linjen:

proc reg data=pemax;
    model pemax=age sex height weight bmp fev1 rv frc tlc 
    / VIF tol stb;
run;

Disse størrelser kan kun fås frem ved brug af proc reg, dvs her dur det ikke med proc glm.

8.8 Prædiktion

Eksempel: Vietnam veteraner (IQ og lungefunktion, uge 2)

Som beskrevet i kap 8.5.1 kan vi vha en output-linje i proc glm bestemme prædikterede værdier for hvert individ i datasættet (den variabel vi kaldte yhat). Vi ønsker dog ofte at prædiktere for en værdi, som ikke nødvendigvis findes i datasættet - og vil gerne have usikkerheden (CI) såvel som prædiktionsinterval for denne værdi også.

8.8.1 Prædiktion incl CI

Vi kan få prædiktionen incl CI frem på to måder:

  1. Ved at flytte y-aksen hen i den værdi, vi ønsker at prædiktere for (e.g. iq=100). Det gør vi i praksis ved at definere en ny variabel iq100=iq-100; og benytte den i analysen i stedet for iq:
data viet;
    set viet;
    iq100=iq-100;
run;
proc glm data=viet;
    model FEV1 = iq100 / clparm;
run;

I denne model svarer interceptet til niveauet for de, som har den forklarende variabel lige 0 (iq100=0), dvs iq=100. Vi aflæser derfor direkte prædiktionen svarende til iq på 100 i intercept-linjen i output - incl CI.

  1. Vha en estimate linje. Her skal vi først angive en tekst vi ønsker SAS benytter i output, så vi ved hvor vi skal aflæse prædiktionen. Denne tekst skal angives i citationstegn.
    Vi ønsker at beregne a + b x10, svarende til 1x intercept + b x 10, hvor altså b er hældningen hørende til iq. Det angiver vi med intercept 1 iq 100:
proc glm data=viet;
    model FEV1 = iq / clparm;
    estimate "IQ=100" intercept 1 iq 100;
run;

8.8.2 Prædiktionsinterval

Når prædiktionen er bestemt som angivet ovenfor, f.eks. vha en estimate-linje, beregnes prædiktionsinterval som prædiktion +/- 2x Root MSE, idet Root MSE angiver spredningen omkring linjen.

8.9 Splines

Hvis en ikke-lineær sammenhæng ikke kan transformeres til linearitet (f.eks. ved log-transformation af kovariaten) kan man benytte spline-funktioner. Her findes forskellige varianter, den simpleste, som præsenteres her, er lineære splines. Her skæres en kurve op i stykkevis lineære bidder. Fordelene ved at bruge lineære splines i forhold til andre typer af splines er, at vi får en model som er let at rapportere - vi kan give en hældning i hvert af intervallerne.

I eksemplet med væksthormon som funktion af alder gennemgået i uge 7, lavede vi sammenhængen stykkevis lineær i intervaller < 10, 10-12, 12-13, 13-15, > 15 år, dvs kurven knækker ved 10, 12, 13 og 15 (svært at se knæk ved 10 og 12 år):

I praksis fittes en sådan kurve ved først at specificere knækpunkterne og derefter definere en spline-variabel for hvert tidspunkt:

data juul;
     set juul;
     age10=max(age-10,0);
   age12=max(age-12,0);
     age13=max(age-13,0);
     age15=max(age-15,0);
run;

Disse benyttes derefter som kovariater i proc glm eller proc reg sammen med den oprindelige variabel (her age):

proc glm data=juul;
    where age ge 5 and age le 20 and sex="male";
    model ssigf1 = age  age10 age12 age13 age15 / solution clparm;
run;

Fortolkningen af spline-variablene er her: age10 angiver forskellen i hældning i intervallet 10-12 ifht < 10, age12 angiver forskellen i hældning i intervallet 12-13 ifht 10-12 etc. Er alle disse forskelle i hældninger 0, dvs ingen effekt af splinevariablene, har vi kun age tilbage - og dermed linearitet. Dette kan vi teste vha et contrast-statement. Dette minder meget om estimate-linjerne illustreret i kapitel 8.8.1. Her ønsker vi at teste, at 1x koefficienten til age10 er 0, 1x koefficienten til age12 er 1 etc. Det giver et test med 4 frihedsgrader, fordi vi forsøger at teste 4 parametre lig 0. Dette dur ikke i proc reg.

proc glm data=juul;
    where age ge 5 and age le 20 and sex="male";
    model sigf1=age  age10 age12 age13 age15 / solution clparm;
  contrast 'test linearitet' age10 1, age12 1, age13 1, age15 1;
run;

Er der som i eksemplet her ikke linearitet (du finder en p-værdi for linearitet fra contrast-linjen på <.0001), vil vi gerne rapportere hældningen i hvert interval frem for forskellen i hældninger mellem intervaller. Lene viser i slides, hvordan man kan bestemme nye-spline variable til dette formål, men det er meget simplere at finde disse hældninger frem vha estimate-linjer. Hældningen i første interval er jo bestemt af koefficienten til age, andet interval age+age10 etc. Vi kan derfor få hældningerne incl CI frem med:

proc glm data=juul;
    where age ge 5 and age le 20 and sex="male";
      model ssigf1=age  age10 age12 age13 age15 / solution;
    estimate 'haeldn < 10'  age 1;
    estimate 'haeldn 10-12' age 1 age10 1;
    estimate 'haeldn 12-13' age 1 age10 1 age12 1;
    estimate 'haeldn 13-15' age 1 age10 1 age12 1 age13 1; 
    estimate 'haeldn > 15'  age 1 age10 1 age12 1 age13 1 age15 1;
run;

8.10 Test af flere parametre

Eksempel: Vitamin d uge 7, den generelle lineære model.

Vi kan teste vilkårlige kombinationer af parametrene lig 0 i et overordnet test vha contrast-linjer i proc glm (dur ikke i proc reg).

I eksemplet gennemgået i forelæsningerne ønsker vi at teste, om der samlet set er effekt af bmi, sunexp og lvitdintake på log-transformeret vitamin D (lvitd), justeret for country. Vi spørger altså om vi kan udelade disse tre variable og kun beholde country i modellen. Det bliver et test med 4 frihedsgrader (1 for bmi, 2 for sunexp og 1 for lvitdintake). I opbygningen af testet håndterer vi kvantitative og kvalitative (class) variable lidt forskelligt. Her specificerer vi at

  • 1x koefficienten til bmi skal være 0 (bmi 1)
  • koefficienten hørende til første niveau af sunexp skal være 0. sunexp har 3 niveauer, 3. niveau er referencegruppe og har derfor en koefficient på 0. Dermed skal vi angive, at koefficienten hørende til 1. niveau skal være lig koefficienten hørende til 3. niveau - det angiver vi som om at 1x koefficient til 1. niveau minus 1x koefficienten til 3. niveau skal være 0 med syntaksen sunexp 1 0 -1
  • koefficienten hørende til 2. niveau skal også være 0 / lig koefficient hørende til 3. niveau, dvs sunexp 0 1 -1
  • 1x koefficienten til lvitdintake skal være 0 (lvitdintake 1)

I contrast-linjen separerer vi de enkelte hypoteser med komma. For den kvalitative variabel er det ligegyldigt om vi skriver sunexp 1 0 -1 eller sunexp -1 0 1 - der skal placeres ialt to 1-taller ud for de parametre vi vil sammenligne, det ene 1-tal skal være negativt.

data vit;
   set vit;
   lvitd=log2(vitd);
   lvitdintake=log2(vitdintake);
run;

proc glm data=vit;
    class country sunexp;
    model lvitd=bmi sunexp lvitdintake country / solution clparm;

    contrast "fjern alle tre paa en gang"
        bmi 1,
        sunexp 1 0 -1,
        sunexp 0 1 -1,
        lvitdintake 1;
run;