Thuis
Contacten

    Hoofdpagina


Docent: drs. Rob Flohr

Dovnload 2 Mb.

Docent: drs. Rob Flohr



Pagina18/25
Datum05.12.2018
Grootte2 Mb.

Dovnload 2 Mb.
1   ...   14   15   16   17   18   19   20   21   ...   25

> 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=c(rep("A",9),rep("B",9),rep("C",9))

> migraine=data.frame(pain,drug)

> migraine

pain drug

1 4 A

2 5 A


3 4 A

4 3 A


5 2 A

6 4 A


7 3 A

8 4 A


9 4 A

10 6 B


11 8 B

12 4 B


13 5 B

14 4 B


15 6 B

16 5 B


17 8 B

18 6 B


19 6 C

20 7 C


21 6 C

22 6 C


23 7 C

24 5 C


25 6 C

26 5 C


27 5 C

> plot(pain~drug,data=migraine)




> results=aov(pain~drug,data=migraine)

> summary(results)

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



> pairwise.t.test(pain,drug,p.adjust="bonferroni")
Pairwise comparisons using t tests with pooled SD
data: pain and drug
A B

B 0.00119 -

C 0.00068 1.00000
P value adjustment method: bonferroni

>


=============================================================================

Bayesian counterpart using MCMCpack

R version 2.15.2 (2012-10-26) -- "Trick or Treat"

Copyright (C) 2012 The R Foundation for Statistical Computing

ISBN 3-900051-07-0

Platform: x86_64-w64-mingw32/x64 (64-bit)


R is free software and comes with ABSOLUTELY NO WARRANTY.

You are welcome to redistribute it under certain conditions.

Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.

Type 'contributors()' for more information and

'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or

'help.start()' for an HTML browser interface to help.

Type 'q()' to quit R.
[Previously saved workspace restored]
> 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=c(rep("A",9),rep("B",9),rep("C",9))

> migraine=data.frame(pain,drug)

> migraine

pain drug

1 4 A

2 5 A


3 4 A

4 3 A


5 2 A

6 4 A


7 3 A

8 4 A


9 4 A

10 6 B


11 8 B

12 4 B


13 5 B

14 4 B


15 6 B

16 5 B


17 8 B

18 6 B


19 6 C

20 7 C


21 6 C

22 6 C


23 7 C

24 5 C


25 6 C

26 5 C


27 5 C

> local({pkg <- select.list(sort(.packages(all.available = TRUE)),graphics=TRUE)

+ if(nchar(pkg)) library(pkg, character.only=TRUE)})

Loading required package: coda

Loading required package: lattice

Loading required package: MASS

##

## Markov Chain Monte Carlo Package (MCMCpack)



## Copyright (C) 2003-2014 Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park

##

## Support provided by the U.S. National Science Foundation



## (Grants SES-0350646 and SES-0350613)

##

Warning messages:



1: package ‘MCMCpack’ was built under R version 2.15.3

2: package ‘coda’ was built under R version 2.15.3



> library(MCMCpack)

> m8=MCMCregress(pain~drug,data=migraine)

> summary(m8)
Iterations = 1001:11000

Thinning interval = 1

Number of chains = 1

Sample size per chain = 10000


1. Empirical mean and standard deviation for each variable,

plus standard error of the mean:


Mean SD Naive SE Time-series SE

(Intercept) 3.665 0.3793 0.003793 0.003834

drugB 2.106 0.5367 0.005367 0.005367

drugC 2.222 0.5338 0.005338 0.005615


2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%

(Intercept) 2.9358 3.419 3.663 3.914 4.415

drugB 1.0363 1.755 2.112 2.464 3.158

drugC 1.1673 1.878 2.224 2.572 3.252




Note that the intercept equals the group mean of the baseline drug (drug A), being 3.67 and that the regression coefficient for drug B equals the difference between the mean of drug B and drug A (5.78 - 3.67 = 2.11), and the same counts of course for drug C (5.89 - 3.67 = 2.22)!!

Keeping drug A as baseline, drug B as well as drug C differ significantly from drug A as the zero-value (of the regression coefficient) is not included in the 95% prediction interval.




> plot(m8,trace=FALSE)

>


-------------------------------------------------------------------------------------------------------------

Jacob Cohen (1923-1998), Amerikaans psycholoog ontwikkelde Cohen's kappa coëfficiënt : een voor toeval gecorrigeerde maat voor de overeenkomst tussen de scores van beoordelaars ('raters'). De mate van geobserveerde overeenstemming wordt gerelateerd aan en aldus gecorrigeerd voor de mate van overeenstemming die je op grond van toeval/kansberekening zou verwachten.

Voorbeeld:
Dertig scripties worden beoordeeld door twee beoordelaars met voldoende of onvoldoende.
De resultaten zijn:









B

B







voldoende

onvoldoende

A

voldoende

10

2

A

onvoldoende

5

13

'observed agreement' : (10+13)/30 = 0.767

Berekenen van kansen:

beoordelaar A geeft 12 van de 30 scripties voldoende -> 12/30


beoordelaar B geeft 15 van de scripties voldoende -> 15/30

Kans dat A én B voldoende geven = 12/30 x 15/30

Op dezelfde manier krijg je de kans dat beiden onvoldoende geven, 18/30 x 15/30

-> 'overall agreement' of

'expected agreement' : (12/30 x 15/30) + (18/30 x 15/30) = 0.50

uitgeschreven:




-> Cohen's kappa: (0.767 - 0.50) / ( 1 - 0.50) = 0.534


Gehanteerde grenzen:


< 0.40 : 'poor agreement'
0.4 < kappa < 0.75 : 'fair to good agreement'
> 0.75 : 'excellent agreement'

Conclusie: overeenstemming is matig


Tot zover ging het over beschrijvende statistiek.

Nu een voorbeeld in het kader van 'statistical inference' met behulp van WinBUGS:

( = de mate waarin de objectieve methode, kortweg O.M., een 1 geeft; = de mate waarin de surrogaatmethode, kortweg S.M., een 1 geeft wanneer de O.M. een 1 geeft, = de mate waarin de S.M. een 0 (nul) geeft wanneer de O.M. een 0 geeft.

->

#Cohen's kappa om de mate van overeenkomst tussen een surrogaat methode en een objectieve (gouden standaard methode) te bepalen. Voor 233 kinderen wordt een goedkope grieptest gebruikt. Van de 18 kinderen die daadwerkelijk griep hebben identificeert deze test er 14 als zodanig en 4 gevallen van griep worden gemist. Van de 215 kinderen die geen griep hebben geeft de test in 210 gevallen de correcte uitslag en in 5 gevallen niet. Welke uitspraak over deze goedkope grieptest (surrogaat methode) kan gedaan worden?

model{

#Priors:


alpha~dbeta(1,1)

beta~dbeta(1,1)

gamma~dbeta(1,1)
#De kans dat voor een bepaalde trekking ('count') zowel O.M. als S.M. een 1 geven (pi[1]), enz.

pi[1]<-alpha*beta

pi[2]<-alpha*(1-beta)

pi[3]<-(1-alpha)*(1-gamma)

pi[4]<-(1-alpha)*(gamma)
#Likelihood:de data - bestaande uit het aantal keren dat O.M. en S.M. beide een 1 geven (a), het aantal keren dat O.M. een 1 geeft en S.M. een 0 (b), etc. - worden gegenereerd door een multinomiale verdeling met n proefnemingen ('trials') waarbij voor elke 'trial' de kans op een a-uitkomst pi[1] is, op een b-uitkomst pi[2], enz.
y[1:4]~dmulti(pi[],233)
#Mate van overeenstemming : O.M. beslist 1 én S.M. beslist 1 (alpha*beta) of O.M. beslist 0 én S.M. beslist 0 (1-alpha)*gamma:
xi<-alpha*beta+(1-alpha)*gamma
#De 'overall' kans op overeenstemming (agreement): de kans dat O.M. tot 1 (a + b) beslist én S.M. tot 1 beslist (a + c) = (a + b)(a + c) of de kans dat S.M. tot 0 beslist (b + d) én dat O.M. tot 0 beslist (c + d), dus samen:
(a + b)(a + c) + (b + d)(c + d)

psi<-(pi[1]+pi[2])*(pi[1]+pi[3])+(pi[2]+pi[4])*(pi[3]+pi[4])


kappa<-(xi-psi)/(1-psi)
}
list(y=c(14,4,5,210))

OUTPUT:


model is syntactically correct

data loaded

model compiled

initial values generated, model initialized


Time series






Kernel density










Node statistics


1   ...   14   15   16   17   18   19   20   21   ...   25

  • > plot(paindrug,data=migraine)
  • > results=aov(paindrug,data=migraine) > summary(results)
  • > pairwise.t.test(pain,drug,p.adjust="bonferroni")
  • _migraine=data.frame(pain,drug)">> 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=c(rep("A",9),rep("B",9),rep("C",9)) > migraine=data.frame(pain,drug)
  • > library(MCMCpack) > m8=MCMCregress(paindrug,data=migraine) > summary(m8)
  • > plot(m8,trace=FALSE)
  • statistical inference

  • Dovnload 2 Mb.