Gebruik de tutorial om deze opdrachten te maken.

Question 1: Coagulation

Beschrijving data

Blood coagulation times by diet

The example dataset we will use is a set of 24 blood coagulation times. 24 animals were randomly assigned to four different diets and the samples were taken in a random order. This data comes from Box, Hunter, and Hunter (1978). . This data frame contains the following columns

  • coag is the coagulation time in seconds
  • diet are the diet type, here A,B,C or D

Vragen

Q1A. Wat voor type data zijn coag en diet ?

Coag is numeriek, diet is categorisch.

Q1B. Wat voor type statistische model(len) kunnen we gebruiken voor de analyse van deze data?

ANOVA

Q1C. Welke variabele is de response en welke de verklarende variabele?

Coag is de response, diet is de verklarende variabele

Q1D. Wat is de null-hypothese

De verschillende diets hebben geen invloed op de blood coagulation.

Inlezen van de data.

In deze vraag gebruiken we data uit het faraway pakket. Dit pakket moet je eerst downloaden als dat nog niet eerder is gedaan. Je kan controleren of het faraway pakket al is geinstaleerd via de Package tab in het venster rechtonder en type faraway in het zoekveld. Als deze niet is geinstalleerd kan je de install knop aanklikken in hetzelfde venster.

library(faraway)
library(car)
## 
## Attaching package: 'car'
## The following objects are masked from 'package:faraway':
## 
##     logit, vif
data(coagulation)
attach(coagulation)
str(coagulation)
## 'data.frame':    24 obs. of  2 variables:
##  $ coag: num  62 60 63 59 63 67 71 64 65 66 ...
##  $ diet: Factor w/ 4 levels "A","B","C","D": 1 1 1 1 2 2 2 2 2 2 ...
summary(coagulation)
##       coag       diet 
##  Min.   :56.00   A:4  
##  1st Qu.:61.75   B:6  
##  Median :63.50   C:6  
##  Mean   :64.00   D:8  
##  3rd Qu.:67.00        
##  Max.   :71.00
plot(coag~diet)

pairwise.t.test(coag,diet)
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  coag and diet 
## 
##   A       B       C      
## B 0.01141 -       -      
## C 0.00090 0.31755 -      
## D 1.00000 0.00345 0.00014
## 
## P value adjustment method: holm
model<-aov(coag~diet)
str(model)
## List of 13
##  $ coefficients : Named num [1:4] 6.10e+01 5.00 7.00 2.99e-15
##   ..- attr(*, "names")= chr [1:4] "(Intercept)" "dietB" "dietC" "dietD"
##  $ residuals    : Named num [1:24] 1 -1 2 -2 -3 ...
##   ..- attr(*, "names")= chr [1:24] "1" "2" "3" "4" ...
##  $ effects      : Named num [1:24] -3.14e+02 5.66 1.40e+01 4.88e-15 -3.53 ...
##   ..- attr(*, "names")= chr [1:24] "(Intercept)" "dietB" "dietC" "dietD" ...
##  $ rank         : int 4
##  $ fitted.values: Named num [1:24] 61 61 61 61 66 66 66 66 66 66 ...
##   ..- attr(*, "names")= chr [1:24] "1" "2" "3" "4" ...
##  $ assign       : int [1:4] 0 1 1 1
##  $ qr           :List of 5
##   ..$ qr   : num [1:24, 1:4] -4.899 0.204 0.204 0.204 0.204 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:24] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:4] "(Intercept)" "dietB" "dietC" "dietD"
##   .. ..- attr(*, "assign")= int [1:4] 0 1 1 1
##   .. ..- attr(*, "contrasts")=List of 1
##   .. .. ..$ diet: chr "contr.treatment"
##   ..$ qraux: num [1:4] 1.2 1.1 1.13 1.27
##   ..$ pivot: int [1:4] 1 2 3 4
##   ..$ tol  : num 1e-07
##   ..$ rank : int 4
##   ..- attr(*, "class")= chr "qr"
##  $ df.residual  : int 20
##  $ contrasts    :List of 1
##   ..$ diet: chr "contr.treatment"
##  $ xlevels      :List of 1
##   ..$ diet: chr [1:4] "A" "B" "C" "D"
##  $ call         : language aov(formula = coag ~ diet)
##  $ terms        :Classes 'terms', 'formula'  language coag ~ diet
##   .. ..- attr(*, "variables")= language list(coag, diet)
##   .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. ..$ : chr [1:2] "coag" "diet"
##   .. .. .. ..$ : chr "diet"
##   .. ..- attr(*, "term.labels")= chr "diet"
##   .. ..- attr(*, "order")= int 1
##   .. ..- attr(*, "intercept")= int 1
##   .. ..- attr(*, "response")= int 1
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. ..- attr(*, "predvars")= language list(coag, diet)
##   .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "factor"
##   .. .. ..- attr(*, "names")= chr [1:2] "coag" "diet"
##  $ model        :'data.frame':   24 obs. of  2 variables:
##   ..$ coag: num [1:24] 62 60 63 59 63 67 71 64 65 66 ...
##   ..$ diet: Factor w/ 4 levels "A","B","C","D": 1 1 1 1 2 2 2 2 2 2 ...
##   ..- attr(*, "terms")=Classes 'terms', 'formula'  language coag ~ diet
##   .. .. ..- attr(*, "variables")= language list(coag, diet)
##   .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. ..$ : chr [1:2] "coag" "diet"
##   .. .. .. .. ..$ : chr "diet"
##   .. .. ..- attr(*, "term.labels")= chr "diet"
##   .. .. ..- attr(*, "order")= int 1
##   .. .. ..- attr(*, "intercept")= int 1
##   .. .. ..- attr(*, "response")= int 1
##   .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. .. ..- attr(*, "predvars")= language list(coag, diet)
##   .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "factor"
##   .. .. .. ..- attr(*, "names")= chr [1:2] "coag" "diet"
##  - attr(*, "class")= chr [1:2] "aov" "lm"
plot(model$residuals~coag)

plot(model, which=1)

fligner.test(coag,diet)
## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  coag and diet
## Fligner-Killeen:med chi-squared = 1.5559, df = 3, p-value = 0.6694
plot(model, which=2)

qqPlot(model$residuals)

summary(model)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## diet         3    228    76.0   13.57 4.66e-05 ***
## Residuals   20    112     5.6                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(model)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = coag ~ diet)
## 
## $diet
##     diff         lwr       upr     p adj
## B-A    5   0.7245544  9.275446 0.0183283
## C-A    7   2.7245544 11.275446 0.0009577
## D-A    0  -4.0560438  4.056044 1.0000000
## C-B    2  -1.8240748  5.824075 0.4766005
## D-B   -5  -8.5770944 -1.422906 0.0044114
## D-C   -7 -10.5770944 -3.422906 0.0001268

Analyseer de data door gebruik te maken van de statistische werkvolgorde, zoals weergegevn in de ANOVA tutorial, en beantwoordt de volgende vragen.

Q1D. Zijn de aannamens voor het uitvoeren van een ANOVA correct? Ja

Q1E. Mogen we de null-hypothese verwerpen? Ja

Q1F. Welke waarde heeft de test-statistiek?

Q1G. Welke diet-afhankelijke coagulatie-tijden verschillen statistisch significant van elkaar? A-B, A-C, B-D, C-D

Question 2 Rabbit data

Beschrijving van de data

A nutritionist studied the effects of six diets, on weight gain of domestic rabbits. From past experience with sizes of litters, it was felt that only 3 uniform rabbits could be selected from each available litter. There were ten litters available forming blocks of size three. In this analysis we only investigate the effect of diets on weight gained.

  • treat is the diet
  • gain is the weight gain

Vragen.

Q2A. Wat voor type data zijn Treat en Weight?

Treat is categorisch, Weight is numeriek

Q2B. Wat voor een type model(len) kan je gebruiken om deze data te analyseren? Verklaar.

ANOVA, want hierin werk je ook met twee soorten variabelen, namelijk categorisch en numeriek

Q2C. Wat zijn de de response variabele en verklarende variabele in deze dataset ? Treat is verklarend, Weight is response.

De data komen weer uit het faraway pakket.

data(rabbit)
attach(rabbit)
summary(rabbit)
##  treat      gain           block   
##  a:5   Min.   :20.60   b1     : 3  
##  b:5   1st Qu.:35.73   b10    : 3  
##  c:5   Median :39.30   b2     : 3  
##  d:5   Mean   :38.47   b3     : 3  
##  e:5   3rd Qu.:42.08   b4     : 3  
##  f:5   Max.   :51.90   b5     : 3  
##                        (Other):12
plot(gain~treat)

pairwise.t.test(gain,treat)
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  gain and treat 
## 
##   a    b    c    d    e   
## b 0.95 -    -    -    -   
## c 1.00 1.00 -    -    -   
## d 1.00 1.00 1.00 -    -   
## e 0.51 1.00 1.00 1.00 -   
## f 1.00 0.88 1.00 1.00 0.46
## 
## P value adjustment method: holm
model<-aov(gain~treat)
plot(model)

summary(model)
##             Df Sum Sq Mean Sq F value Pr(>F)
## treat        5  293.4   58.68   1.886  0.134
## Residuals   24  746.5   31.10
TukeyHSD(model)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = gain ~ treat)
## 
## $treat
##      diff        lwr       upr     p adj
## b-a -6.46 -17.366149  4.446149 0.4657149
## c-a -3.20 -14.106149  7.706149 0.9408659
## d-a -5.80 -16.706149  5.106149 0.5790395
## e-a -7.82 -18.726149  3.086149 0.2667934
## f-a  0.28 -10.626149 11.186149 0.9999995
## c-b  3.26  -7.646149 14.166149 0.9363140
## d-b  0.66 -10.246149 11.566149 0.9999636
## e-b -1.36 -12.266149  9.546149 0.9987465
## f-b  6.74  -4.166149 17.646149 0.4200188
## d-c -2.60 -13.506149  8.306149 0.9751430
## e-c -4.62 -15.526149  6.286149 0.7769135
## f-c  3.48  -7.426149 14.386149 0.9177686
## e-d -2.02 -12.926149  8.886149 0.9919355
## f-d  6.08  -4.826149 16.986149 0.5303186
## f-e  8.10  -2.806149 19.006149 0.2340703

Analyseer de data door gebruik te maken van de tutorial en beantwoord de volgende vragen.

Q1D. Zijn de aannamens voor het uitvoeren van een ANOVA correct? Verklaar. Ik neem aan van wel

Q1E. Mogen we de null-hypothese verwerpen? Verklaar. Nee, de p-value>0,05

Q1F. Welke niveaus of de diet factor verschillen statistisch significant van elkaar? Geen niveay verschilt significant

Q3

A drug company tested three formulations of a pain relief medicine for migraine headache sufferers. For the experiment 27 volunteers were selected and 9 were randomly assigned to one of three drug formulations. The subjects were instructed to take the drug during their next migraine headache episode and to report their pain on a scale of 1 to 10 (10 being most pain).

Inlezen data

we voeren de data nu met de hand in.

Pain = c(4, 5, 4, 3, 2, 4, 3, 4, 4, 6, 8, 4, 5, 4, 6, 5, 8, 6, 6, 7, 6, 6, 7, 5, 6, 5, 5)
Drug = gl(3,9,labels=c("A", "B", "C"))
DATA<-data.frame(Pain, Drug)
summary(DATA)
##       Pain       Drug 
##  Min.   :2.000   A:9  
##  1st Qu.:4.000   B:9  
##  Median :5.000   C:9  
##  Mean   :5.111        
##  3rd Qu.:6.000        
##  Max.   :8.000
plot(Pain~Drug)

pairwise.t.test(Pain,Drug)
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  Pain and Drug 
## 
##   A       B      
## B 0.00079 -      
## C 0.00068 0.83042
## 
## P value adjustment method: holm
model<-aov(Pain~Drug)
plot(model)

summary(model)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Drug         2  28.22  14.111   11.91 0.000256 ***
## Residuals   24  28.44   1.185                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(model)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Pain ~ Drug)
## 
## $Drug
##          diff        lwr      upr     p adj
## B-A 2.1111111  0.8295028 3.392719 0.0011107
## C-A 2.2222222  0.9406139 3.503831 0.0006453
## C-B 0.1111111 -1.1704972 1.392719 0.9745173

Gebuik de tutorial om de data te analyseren in de context van de volgende vraag: Is er een verschil in pijn bestrijding door de verschillende medicijnen bij migraine?

Q4: Typical exam questions

De tabel hieronder is gegeven maar er missen wat data (?).

             Df Sum Sq Mean Sq F value  Pr(>F)  
VerklarendeF  3    ?1    2.199     ?3   0.0945 .
Residuals    36  34.53     ?2

Q4.A Bereken de waarden die op de plaatsen van de vraagtekens moeten staan (?).

Q4.B Welke eigenschappen van de t-statistiek en F-statistiek komen overeen?

Simulatie ANOVA experiment.

Eerst maken we de factor. We hebben twee versies nodig van de factor, een numerieke en een categorische. De numerieke is alleen een lijst van de niveaus gemiddelde.

VerklarendeN<-c(rep(2,10),rep(2.2,10),rep(2.8,10),rep(3.4,10) )

dan de categorische versie wat een naam aan de gemiddelde koppeld.

VerklarendeF<-gl(4,10, labels=c("hokus", "Pokus", "Pilatus", "Pas"))

We maken nu een ANOVA model alleen met de gemiddelde en plotten het resultaat

Response<-VerklarendeN
plot(Response~VerklarendeF)

We zien alleen de verschillende factor gemiddelde liggen, de between-group variance.

Nu voegen we de rest term toe ofwel we geven de within-group varianec een waarde. We voegen daarom een stochatische variabele toe met een gemiddelde van 0 en een bepaalde standaard deviatie (sd) hier 1.2.

ResponseR<-VerklarendeN+rnorm(40,0,1.2)
Dat2<-data.frame(ResponseR,VerklarendeF)
summary(Dat2)
##    ResponseR       VerklarendeF
##  Min.   :-1.885   hokus  :10   
##  1st Qu.: 1.739   Pokus  :10   
##  Median : 2.555   Pilatus:10   
##  Mean   : 2.587   Pas    :10   
##  3rd Qu.: 3.805                
##  Max.   : 5.357
plot(ResponseR~ VerklarendeF)
points(ResponseR~ VerklarendeF)

Voer een analyse uit op deze data en pas de gemiddelde zo aan dat de test net significant is. Hiervoor moeten ze de VerklarendeF aanpassen. Letop ik gebruik geen set.seed.