Thuis
Contacten

    Hoofdpagina


Docent: drs. Rob Flohr

Dovnload 2 Mb.

Docent: drs. Rob Flohr



Pagina10/25
Datum05.12.2018
Grootte2 Mb.

Dovnload 2 Mb.
1   ...   6   7   8   9   10   11   12   13   ...   25


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

Huiswerkopgave voor les 4


I) Lees de oratie van Prof. dr. Jan-Willem Romeijn en formuleer vragen die je tijdens de les kunt stellen; de tekst van de oratie staat op Blackboard.


II) Los het volgende vraagstuk op:

Lineaire regressie

How well does the number of beers a student drinks predict his or her blood alcohol content?


Sixteen student volunteers at Ohio State University drank a randomly assigned number of 12-ounce cans of beer. The students were equally divided between men and women and differed in weight and usual drinking habits.
Thirty minutes later, a police officer measured their blood alcohol content (BAC). Here are the data:

Student

1

2

3

4

5

6

7

8

Beers

BAC


5

0.10


2

0.03


9

0.19


8

0.12


3

0.04


7

0.095


3

0.07


5

0.06


Student

9

10

11

12

13

14

15

16

Beers

BAC


3

0.032


5

0.05


4

0.07


6

0.10


5

0.085


7

0.09


1

0.01


4

0.05


Opdracht:

Bepaal de regressiecoëfficiënten en toets die op significantie volgens
a) de klassieke (frequentistische) statistiek
b) de Bayesiaanse statistiek m.b.v. MCMCpack (zonder priors, dus alleen het model en de data invullen)

Uitwerking:

a) klassiek m.b.v. R:

> beers=c(5,2,9,8,3,7,3,5,3,5,4,6,5,7,1,4)



> bac=c(0.10,0.03,0.19,0.12,0.04,0.095,0.07,0.06,0.02,0.05,0.07,0.10,0.085,0.09,0.01,0.05)

> plot(beers,bac)


b)
> lm(bac~beers)


Call:

lm(formula = bac ~ beers)


Coefficients:

(Intercept) beers



-0.01270 0.01796


> regres=lm(bac~beers)

> abline(regres)



c)
> cor(beers,bac)

[1] 0.8943381



> 0.8943381^2

[1] 0.7998406


d)
Not surprisingly, we find that BAC increases as BAC increases; the relationship is quite strong, with beer consumption explaining 80% of the variation in BAC.


e)
> model=lm(bac~beers)



> summary(model)
Call:

lm(formula = bac ~ beers)


Residuals:

Min 1Q Median 3Q Max

-0.027118 -0.017350 0.001773 0.008623 0.041027
Coefficients:

Estimate Std. Error t value Pr(>|t|)



(Intercept) -0.012701 0.012638 -1.005 0.332

beers 0.017964 0.002402 7.480 2.97e-06 ***

---


Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.02044 on 14 degrees of freedom

Multiple R-squared: 0.7998, Adjusted R-squared: 0.7855

F-statistic: 55.94 on 1 and 14 DF, p-value: 2.969e-06
> confint(model)

2.5 % 97.5 %

(Intercept) -0.03980535 0.01440414

beers 0.01281262 0.02311490



> summary.aov(model)

Df Sum Sq Mean Sq F value Pr(>F)

beers 1 0.02337 0.023375 55.94 2.97e-06 ***

Residuals 14 0.00585 0.000418

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1



N.B. 0.02337/(0.02337+0.00585)=0.7998

>


Conclusion: there is very strong evidence that drinking more beers increases BAC in the population)

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

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]
> utils:::menuInstallPkgs()

--- Please select a CRAN mirror for use in this session ---

trying URL 'http://cran-mirror.cs.uu.nl/bin/windows/contrib/2.15/MCMCpack_1.3-3.zip'

Content type 'application/zip' length 2950090 bytes (2.8 Mb)

opened URL

downloaded 2.8 Mb


package ‘MCMCpack’ successfully unpacked and MD5 sums checked
The downloaded binary packages are in

C:\Users\Eigenaar\AppData\Local\Temp\RtmpuchU0P\downloaded_packages

> 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)

> help(package=MCMCpack)

starting httpd help server ... done

> beers=c(5,2,9,8,3,7,3,5,3,5,4,6,5,7,1,4)

> bac=c(0.10,0.03,0.19,0.12,0.04,0.095,0.07,0.06,0.02,0.05,0.07,0.10,0.085,0.09,0.01,0.05)

> m=data.frame(beers,bac)

> m

beers bac

1 5 0.100

2 2 0.030

3 9 0.190

4 8 0.120

5 3 0.040

6 7 0.095

7 3 0.070

8 5 0.060

9 3 0.020

10 5 0.050

11 4 0.070

12 6 0.100

13 5 0.085

14 7 0.090

15 1 0.010

16 4 0.050



> m6=MCMCregress(bac~beers,data=m)

> summary(m6)
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) -0.0124923 0.014940 1.494e-04 1.470e-04

beers 0.0179221 0.002848 2.848e-05 2.773e-05
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%

(Intercept) -0.0417037 -0.0219154 -0.0124212 -0.0031053 0.017676



(de nul-waarde zit in het 95%-interval, dus niet significant)


beers 0.0122021 0.0161481 0.0179487 0.0197013 0.023553



(de nul-waarde zit niet in het 95%-interval, dus wel significant)






> plot(m6,trace=FALSE)

>

Nogmaals met MCMCpack (14 april 2015)

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]
> beers

[1] 5 2 9 8 3 7 3 5 3 5 4 6 5 7 1 4

> bac

[1] 0.100 0.030 0.190 0.120 0.040 0.095 0.070 0.060 0.032 0.050 0.070 0.100



[13] 0.085 0.090 0.010 0.050

> library(MCMCpack)

Loading required package: coda

Loading required package: lattice

Loading required package: MASS

##

## Markov Chain Monte Carlo Package (MCMCpack)



## Copyright (C) 2003-2015 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

3: package ‘lattice’ was built under R version 2.15.3

> posterior<-MCMCregress(bac~beers)

> plot(posterior)

> raftery.diag(posterior)


Quantile (q) = 0.025

Accuracy (r) = +/- 0.005

Probability (s) = 0.95

Burn-in Total Lower bound Dependence

(M) (N) (Nmin) factor (I)

(Intercept) 3 4028 3746 1.080

beers 2 3851 3746 1.030

sigma2 2 3680 3746 0.982


> summary(posterior)
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) -0.0103031 0.0145196 1.452e-04 1.428e-04

beers 0.0176230 0.0027678 2.768e-05 2.695e-05

sigma2 0.0005433 0.0002503 2.503e-06 2.892e-06


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

(Intercept) -0.03869 -0.0194611 -0.0102341 -0.001180 0.019016

beers 0.01206 0.0158990 0.0176489 0.019352 0.023095

sigma2 0.00025 0.0003779 0.0004868 0.000638 0.001157


>

Les 4

ONDERDEEL I)
-terugkomen op uitgeversvraagstuk (zie bijlage aan het eind van dit bestand)

a) i.v.m. essentie Bayesiaanse statistiek: tot nog toe overwegend voorbeelden besproken met non-informatieve prior waardoor getalsmatig vrijwel dezelfde uitkomsten werden verkregen als bij toepassing van de klassieke (frequentistische) statistiek.


Dit voorbeeld laat zien: zelfs bij gelijksoortige uitkomsten is er een belangrijk verschil: het gaat om de posterior waarmee kansuitspraken over een parameter (onbekende grootheid die van belang is, hier: gemiddeld aantal typografische fouten per 100 blz.) kunnen worden gedaan.

Bovendien: het uitgeversprobleem is een exemplarisch voorbeeld van Statistical Decision Making: er wordt gebruik gemaakt van een norm (hier: 10) , van een kansuitspraak over een parameter, en van een beslissingscriterium (wanneer wel/niet maatregelen nemen om aantal fouten te verminderen?)

b) ivm posterior predictive distribution (alle informatie over de parameter wordt gebruikt).

Voorbeeld van de munt (les 1): de posterior predictive distribution wordt in dit voorbeeld gevonden door:



Wat we hier in feite doen is een groot aantal parameterwaarden aselect uit de posterior verdeling trekken en dan elke waarde gebruiken om het aantal successen (k) te voorspellen.

Het concept van de posterior predictive distribution en de uitwerking m.b.v. WinBUGS zal tijdens les 5 aan de orde komen, onder meer in de context van de vraag hoe adequaat een statistisch model is in het beschrijven van de werkelijkheid.
De posterior predictive distribution kan gezien worden als een verzameling van voorspellingen ten aanzien van welke data op basis van het model verwacht kunnen worden.

ONDERDEEL II)


- oratie Romeijn:

Discussie in de groep



1   ...   6   7   8   9   10   11   12   13   ...   25

  • > lm(bacbeers) Call: lm(formula = bac beers) Coefficients: (Intercept) beers -0.01270 0.01796 > regres=lm(bacbeers)
  • > cor(beers,bac)
  • (Intercept) -0.012701 0.012638 -1.005 0.332 beers 0.017964 0.002402 7.480 2.97e-06 ***
  • 0.7998
  • > library(MCMCpack)
  • > m6=MCMCregress(bacbeers,data=m) > summary(m6)
  • -0.0124923
  • 0.0122021 0.0161481 0.0179487 0.0197013 0.023553 (de nul-waarde zit niet in het 95%-interval, dus wel significant)

  • Dovnload 2 Mb.