Øvelser dag 4

Rat growth

1. Indlæs data.

Husk først at sætte working directory og læs derefter data ind med kommandoerne:

library(foreign)
d <- read.dbf('ratgrowth.dbf')

Gennemsnittet af responsen (growth) for hver kombination af antibiotic og vitamins:

mean( d$growth[ d$antibiotic=='no' & d$vitamins =='no' ])

[1] 1.19

mean( d$growth[ d$antibiotic=='no' & d$vitamins =='yes' ])

[1] 1.22

mean( d$growth[ d$antibiotic=='yes' & d$vitamins =='no' ])

[1] 1.033333

mean( d$growth[ d$antibiotic=='yes' & d$vitamins =='yes' ])

[1] 1.543333

I tabelform er gruppegennemsnittene:

	vit yes	vit no
antibiotika no	1.22	1.19
antibiotika yes	1.54	1.03

2. Kvantific?r effekten af antibiotika.

I formiddags fandt vi en interaktion mellem vitamin og antibiotika. Det betyder at vi er nødt til at kvantificere effekten af antibiotika for hvert niveau af vitamin (dvs. for hhv. yes og no).

Af tabellen ovenfor kan vi se:

for vit=yes er forskellen mellem yes og no for antibiotika 1.54-1.22=0.32, dvs. en øgning på 0.32.
for vit=no er forskellen mellem yes og no for antibiotika 1.03-1.19=-0.16, dvs. et fald på 0.16.

Vi vil nu prøve at få R til at bestemme tallene (inspireret af slide 27-28):

lm3 <- lm( growth ~ vitamins + antibiotic:vitamins, data=d )
summary( lm3 )


Call:
lm(formula = growth ~ vitamins + antibiotic:vitamins, data = d)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.110000 -0.025000  0.003333  0.016667  0.110000 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)                1.19000    0.03496  34.039 6.06e-10 ***
vitaminsyes                0.03000    0.04944   0.607  0.56082    
vitaminsno:antibioticyes  -0.15667    0.04944  -3.169  0.01322 *  
vitaminsyes:antibioticyes  0.32333    0.04944   6.540  0.00018 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.06055 on 8 degrees of freedom
Multiple R-squared:  0.9336,    Adjusted R-squared:  0.9087 
F-statistic: 37.48 on 3 and 8 DF,  p-value: 4.659e-05

Vi fylder tallene ind i en tabel for at være helt sikre på, at vi forstår output rigtigt / har fittet den rigtige model

	vit yes	vit no
antibiotika no	1.19+0.03	1.19
antibiotika yes	1.19+0.03+0.32	1.19-0.16

Vi bemærker at summen af tallene i denne tabel stemmer overens med gruppegennemsnittene angivet i tabellen i spørgsmål 1 ovenfor.

Vi ser, at for de rotter som har fået vitamin, er effekten af at give antibiotika en øgning på 0.32, p=0.0002.

For de rotter som ikke har fået vitamin, er der en negativ effekt af at give antibiotika, idet effekten er -0.16, p=0.01.

Vi konkluderer altså, at antibiotika kun er gavnligt sammen med vitamin, hvorimod der kan være en skadelig effekt af at give antibiotika uden vitamin.

3. Undersøg normalfordelingsantagelsen

Vi laver et QQ-plot

plot( lm3, which=2 )

og ser umiddelbart ikke tegn på, at normalfordelingsantagelsen skulle være problematisk. Med datasæt af denne størrelse, vil vi dog sjældent være i stand til at finde afvigelser mm. afvigelsen er ekstrem.

4. Undersøg antagelsen om varianshomogenitet

Vi laver dels et grafisk check

plot( lm3, which=1 )

Vi ser, at spredningen er størst i en af grupperne på midten. Men igen er datasættet så lille (3 rotter i hver gruppe) at vi ikke kan konkludere der er forskelle.

Vi supplerer med Bartlett’s test, men er her nødt til først at definere en ny variabel, som angiver de 4 mulige kombinationer af vores forklarende variable, idet Bartlett kun kan køres i den ensidede variansanalysemodel.

Vi starter med at definere den nye variabel:

d$av <- 'nn'
d$av[ d$antibiotic=='yes' & d$vitamins=='no']  <- 'yn'
d$av[ d$antibiotic=='no'  & d$vitamins=='yes'] <- 'ny'
d$av[ d$antibiotic=='yes' & d$vitamins=='yes'] <- 'yy'

Som altid, når der laves en ny variabel, laves en kontrol af at den er rigtigt defineret. Vi tæller først, hvor mange rotter der er for hver kombination af de to forklarende variable

table( d$antibiotic, d$vitamins )

     
      no yes
  no   3   3
  yes  3   3

og sammenligner med den variabel vi lige har defineret og ser at det passer:

table( d$av )


nn ny yn yy 
 3  3  3  3

Vi kan nu udføre Bartlett’s test:

bartlett.test( d$growth ~ d$av )


    Bartlett test of homogeneity of variances

data:  d$growth by d$av
Bartlett's K-squared = 5.7271, df = 3, p-value = 0.1257

Med en p-værdi på 0.13 kan vi således ikke afvise hypotesen om, at varianserne er ens.

Vi konkluderer derfor, at det ikke er problematisk at antage at variansen i de 4 grupper er ens.

5. Interaktionsplot

Først installeres pakken gplots med kommandoen install.packages('gplots'). Dernæst gør vi indholdet af denne pakke tilgængelig med library():

library(gplots)

Vi laver plottet som foreslået i teksten:

plotmeans( growth ~ interaction(vitamins,antibiotic), data=d )

For alle kombinationer af vitamin og antibiotika har vi fået plottet mean incl. 95% konfidensinterval. Vi har fået forbundet alle 4 punkter. Vi vil hellere forbinde punkterne, så vi ser effekten af vitaminer:

plotmeans( growth ~ interaction(vitamins,antibiotic), data=d ,
           connect=list(c(1,2),c(3,4)))

Her ser vi, at der ingen effekt er af at give vitaminer til rotter, der har ikke har fået antibiotika. Omvendt ser vi en øgning for de rotter, som fik antibiotika.

Vi kan også illustrere effekten af antibiotika:

plotmeans( growth ~ interaction(vitamins,antibiotic), data=d ,
           connect=list(c(1,3),c(2,4)))

Her ser vi, at for de rotter, som ikke har fået vitaminer er der en skadelig effekt af at give antibiotika. At dette rent faktisk er signifikant, fandt vi ud af ovenfor. Vi ser også, at effekten af at give antibiotika til rotter, der har fået vitaminer, er gavnlig.

6. En anden parameterisering

Vi prøver nu at køre modellen:

lm4 <- lm( growth ~ antibiotic:vitamins, data=d )
summary( lm4 )


Call:
lm(formula = growth ~ antibiotic:vitamins, data = d)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.110000 -0.025000  0.003333  0.016667  0.110000 

Coefficients: (1 not defined because of singularities)
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)                1.54333    0.03496  44.145 7.65e-11 ***
antibioticno:vitaminsno   -0.35333    0.04944  -7.147 9.74e-05 ***
antibioticyes:vitaminsno  -0.51000    0.04944 -10.315 6.73e-06 ***
antibioticno:vitaminsyes  -0.32333    0.04944  -6.540  0.00018 ***
antibioticyes:vitaminsyes       NA         NA      NA       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.06055 on 8 degrees of freedom
Multiple R-squared:  0.9336,    Adjusted R-squared:  0.9087 
F-statistic: 37.48 on 3 and 8 DF,  p-value: 4.659e-05

og fylder ind i tabellen

	vit yes	vit no
antibiotika no	1.54-0.32	1.54-0.35
antibiotika yes	1.54	1.54-0.51

Vi ser at vi med denne parameterisering får forskelle mellem grupperne i forhold til antibiotika=yes og vitamin=yes. Dvs. alle sammenligninger mod denne reference gruppe. Det kunne vi også få frem ved at bruge vores interaktionsvariabel av ovenfor med referencegruppe yy. Til dette kan man benytte kommandoen relevel():

lm4a <- lm( growth ~ relevel( factor(av), 'yy'), data=d )
summary( lm4a )


Call:
lm(formula = growth ~ relevel(factor(av), "yy"), data = d)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.110000 -0.025000  0.003333  0.016667  0.110000 

Coefficients:
                            Estimate Std. Error t value Pr(>|t|)    
(Intercept)                  1.54333    0.03496  44.145 7.65e-11 ***
relevel(factor(av), "yy")nn -0.35333    0.04944  -7.147 9.74e-05 ***
relevel(factor(av), "yy")ny -0.32333    0.04944  -6.540  0.00018 ***
relevel(factor(av), "yy")yn -0.51000    0.04944 -10.315 6.73e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.06055 on 8 degrees of freedom
Multiple R-squared:  0.9336,    Adjusted R-squared:  0.9087 
F-statistic: 37.48 on 3 and 8 DF,  p-value: 4.659e-05

7. En tredie parameterisering

Vi kører den foreslåede model:

lm5 <- lm( growth ~ antibiotic:vitamins -1, data=d )
summary( lm5 )


Call:
lm(formula = growth ~ antibiotic:vitamins - 1, data = d)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.110000 -0.025000  0.003333  0.016667  0.110000 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
antibioticno:vitaminsno    1.19000    0.03496   34.04 6.06e-10 ***
antibioticyes:vitaminsno   1.03333    0.03496   29.56 1.86e-09 ***
antibioticno:vitaminsyes   1.22000    0.03496   34.90 4.97e-10 ***
antibioticyes:vitaminsyes  1.54333    0.03496   44.15 7.65e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.06055 on 8 degrees of freedom
Multiple R-squared:  0.9985,    Adjusted R-squared:  0.9977 
F-statistic:  1300 on 4 and 8 DF,  p-value: 2.783e-11

Vi ser her, at vi ikke længere har nogen referencegruppe (der er ikke noget intercept og alle 4 kombinationer af de to variable optræder i vores summary). Det vi får ud her er de 4 estimerede middelværdier i hver gruppe (dvs. gruppegennemsnittene). Det er disse 4 tal, som er plottet i interaktionsplottene ovenfor. Vi kan nu få CI for disse 4 gennemsnit ved kommandoen confint()

confint( lm5 )

                              2.5 %   97.5 %
antibioticno:vitaminsno   1.1093814 1.270619
antibioticyes:vitaminsno  0.9527147 1.113952
antibioticno:vitaminsyes  1.1393814 1.300619
antibioticyes:vitaminsyes 1.4627147 1.623952

Dette er faktisk den eneste måde vi kan få R til at beregne konfidensintervallerne i de 4 grupper på - ellers ville vi være nødt til at lave dem i hånden. Bemærk også at e.g konfidensintervallet for kombinationen antibiotic=='no' og vitamins=='no' er noget smallere end det, vi får når vi laver et interaktionsplot. Det skyldes, at interaktionsplottet beregner en SD (og dermed en SEM) i hver af de fire grupper. Dvs hver af disse SD’er er baseret på n=3 observationer. I vores model har vi antaget at variansen er ens i de 4 grupper (som absolut ikke var urimeligt ifølge vores undersøgelser ovenfor). Det betyder at vi får et mere præcist bud på den fælles varians (baseret på de 12 observationer). Omvendt bliver CI’et for yes-yes kombinationen en anelse større (denne gruppe har den mindste variation). Idet vi har præcis lige mange observationer i hver gruppe (3 stk) bliver CI’erne lige brede i alle 4 grupper, når vi bruger den fælles SD modellen estimerer.