Inleiding tot R bij Statistiek voor de Gedragswetenschappen

Inleiding tot R bij Statistiek voor de Gedragswetenschappen Joris Pieters

Peter Theuns

Alain Isaac

Januari 2014

Inhoudsopgave 1 Inleiding 1.1 Over R . . . . . . . . . 1.2 Command line interface 1.3 Enkele begrippen . . . . 1.4 R-functies in deze cursus

. . . .

2 2 2 3 4

2 Basistechnieken 2.1 Rekentechnieken, variabelen en vectoren . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Analyse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Matrixrekenen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5 5 9 15

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

3 Beschrijvende statistiek 19 3.1 Eéndimensionaal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.2 Relaties tussen 2 variabelen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4 Kans- en steekproevenverdelingen 40 4.1 Kansverdelingen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 4.2 Steekproevenverdelingen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 5 Statistische analyse van het gemiddelde 50 5.1 Betrouwbaarheidsinterval . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 5.2 Toetsen van een hypothese omtrent µ0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 5.3 Vergelijken van twee gemiddelden: onafhankelijke steekproeven (between-subjects vergelijking) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 5.4 Vergelijken van twee gemiddelden: gekoppelde data (paired samples, within-subject) . . 54 Oplossingen

56

Appendix A: Overzicht functies

67

Appendix B: Packages

76

Appendix C: Databestanden

77

1

1

Inleiding

1.1

Over R

R is een computertaal en bijhorend softwarepakket dat voornamelijk gebruikt wordt voor statistische berekeningen. In tegenstelling tot de meeste andere statistische paketten is R gratis, en bovendien kan iedere gebruiker het programma aanvullen met (eventueel zelfgemaakte) “packages”1 . Dit is echter niet de enige reden waarom we met R leren werken in plaats van met commerciële software. Waar andere software voornamelijk de nadruk legt op gebruikersgemak wordt van een R-gebruiker verwacht dat hij of zij ook weet waarom een bepaalde procedure gekozen wordt en niet een andere. We zijn van mening dat R niet enkel een gratis programma is dat de berekeningen voor ons doet, maar vooral ook dat het een leeromgeving is die studenten helpt om statistiek beter te begrijpen. R vereist echter wel enige oefening, en dat is het doel van deze handleiding. Eerst en vooral zal je natuurlijk R moeten installeren: Windows: http://cran.freestatistics.org/bin/windows/ Linux: http://cran.freestatistics.org/bin/linux/ OS X: http://cran.freestatistics.org/bin/macosx/

Vanaf hier heeft het enkel zin om verder te lezen als je achter een computer zit waarop R is ge¨ınstalleerd. Zoals je zal zien neemt deze tekst een “tutorial-stijl” aan: Er wordt gezegd wat je moet doen (welke input je R moet geven), en wat de output is die R geeft. Je moet dus telkens proberen zelf te doen wat we in de tekst voordoen. Gebruik liever geen copy/paste, maar typ de input over. Dit betekent echter niet dat je veel gaat bijleren door hersenloos alles over te typen. We moedigen je namelijk aan om zoveel mogelijk te experimenteren met je input door kleine wijzigingen aan te brengen aan de input die hier wordt voorgesteld en na te gaan wat het effect daarvan is. Laat je vooral niet ontmoedigen wanneer dit een “error” oplevert: R begrijpt niet wat je bedoelt als je iets verkeerd hebt ingegeven. Dat is niet erg, je kan veel leren uit fouten. Gebeurt dit, adem dan diep in en uit en probeer opnieuw! Verder zijn er ook oefeningen op het einde van ieder deel. Die zou je steeds moeten kunnen maken aan de hand van wat je daarvoor hebt geleerd. De oplossingen van de meeste oefeningen staan achteraan in deze bundel.

1.2

Command line interface

Je bent wellicht het meest gewend om een computer te besturen met een grafische gebruikersomgeving (Graphical User Interface, GUI): de menu’s, knoppen, tekstvakken, et cetera waarmee je een programma aanstuurt. Dit is echter slechts één methode, een andere manier om een computer te besturen is met een “Command Line Interface” (CLI), waarbij de opdrachten die je aan een programma geeft moeten worden ingetypt. R werkt grotendeels via de CLI. Het grote voordeel van een CLI is de veel grotere flexibiliteit. Stel bijvoorbeeld dat je een programma een bepaalde procedure meerdere keren wilt laten uitvoeren, telkens op een iets andere manier. Zonder CLI kan dat een hele hoop muisgeklik opleveren, terwijl zoiets met gebruik van een CLI niet echt veel meer werk is dan de procedure slechts één keer uit te voeren. Een CLI heeft als nadeel dat je de commando’s (best uit het hoofd) moet kennen. Bovendien is het belangrijk dat je ze juist intypt. In tijden van SMS-taal en spellingcorrectie is de vaardigheid om woorden correct te typen echter niet meer iedereen gegeven. De CLI van R is bovendien hoofdlettergevoelig. Enkel voor het gebruik van spaties geldt enige tolerantie in de R-CLI. Het gebruik van R is dus een stuk eenvoudiger voor mensen die blind kunnen typen. Als je dat niet kunt: oefenen! 1 Zo

zijn er letterlijk duizenden beschikbaar.

2

1.3

Enkele begrippen

Een variabele is een naam die we geven aan een verzameling gegevens. Die verzameling kan leeg zijn, één enkel element bevatten, uit meerdere verschillende soorten elementen bestaan, tot en met een zeer complexe structuur omschrijven.2 Een functie is een deeltaak van een (computer)programma. We noteren een functie steeds met de naam van de functie, gevolgd door ronde haken. In de context van autorijden hebben we bijvoorbeeld de functie starten() die op zichzelf zal bestaan uit functies als draai sleutel() en handrem afzetten(). Een functie kan één of meerdere parameters hebben. Dit zijn variabelen waarop de functie wordt toegepast. Ze worden tussen de ronde haakjes van de functie genoteerd, gescheiden door een komma. Bijvoorbeeld, de functie rijden(startpunt, bestemming) heeft twee parameters die allebei een plaats voorstellen. Parameters kunnen verplicht of optioneel zijn. In tegenstelling tot verplichte parameters moet je optionele parameters niet altijd ingeven. Een optionele parameter wordt in deze tekst, maar ook in de help-functie van R, genoteerd door de naam van de parameter te laten volgen door een gelijkheidsteken en de standaardwaarde van die optionele parameter, ook wel default genoemd. In de functie schakel omhoog(versnellingen = 1) geeft een optionele parameter weer hoeveel versnellingen we omhoog gaan schakelen indien deze parameter niet ingevuld wordt. Of we deze functie nu aanroepen door middel van schakel omhoog() of schakel omhoog(versnellingen = 1), in beide gevallen zal de wagen één versnelling hoger schakelen. Willen we in één keer twee versnellingen hoger schakelen dan roepen we schakel omhoog(versnellingen = 2) aan, enzoverder. Wanneer een parameter optioneel is, terwijl er geen zinvolle standaardwaarde is, dan stellen we de standaardwaarde NULL in3 . Indien een functie zowel verplichte als optionele parameters bevat worden de verplichte parameters meestal eerst geplaatst. Bijvoorbeeld twee verplichte en drie optionele parameters: doe een testrit( auto, stad, maximum snelheid = 120, tijd minuten = 30, radiostation = NULL). Wanneer een functie meerdere optionele parameters heeft zou verwarring kunnen ontstaan: Stel we hebben een functie met meerdere optionele parameters waarvan we er slechts enkele zelf willen invullen en voor de overige de standaardwaarde behouden. Het programma weet nu niet voor welke parameters de ingevulde waarden dienen. Om dit op te lossen geven we bij optionele parameters steeds de naam van deze parameter mee. Zo is het duidelijk dat doe een testrit(Volvo, Brussel, maximum snelheid = 70) mij een testrit laat maken aan maximum 70km/u, en dit voor 30 minuten (standaardwaarde), terwijl doe een testrit(Volvo, Brussel, tijd minuten = 50) mij 50 minuten laat rijden met een maximum van 120km/u.4 Of een parameter verplicht dan wel optioneel is kan soms voor de hand liggen, maar in andere gevallen is dit, net als de standaardwaarde, een arbitraire keuze van de programmeur. Bij verplichte parameters mag de naam van de parameter meegegeven worden maar dit is niet noodzakelijk aangezien de volgorde voldoende duidelijkheid biedt omtrent welke waarden voor welke parameter bedoeld zijn. Functies kunnen overloaded zijn, wat betekent dat één functie met verschillende soorten parameters kan worden aangeroepen. Zo zal lengte(auto) de auto meten van vooraan tot achteraan, terwijl lengte(mens) van onder tot boven wordt gemeten. Overloading zorgt er dus voor dat we niet voor ieder type variabele andere functies moeten kennen.

2 Een

variabele is in de context van een computerprogramma dus een ruimer begrip dan in de statistiek. 6= 0. Terwijl 0 een cijfer is met een waarde (zij het dan de waarde die de afwezigheid van een kenmerk voorstelt), is NULL niets, nada, nothing. 4 Merk op dat verplichte parameters dit probleem niet hebben. Het is echter wel toegelaten om de naam van verplichte parameters toch mee te geven. 3 NULL

3

1.4

R-functies in deze cursus

Omschrijving van de functie. functienaam ( verplichte parameter 1, verplichte parameter 2, . . ., optionele parameter 1 = standaardwaarde 1, optionele parameter 2 = standaardwaarde 2, . . . = . . .)

JorisPieters,PeterTheuns,AlainIsaac

Alle R-functies zullen steeds op dezelfde manier voorgesteld worden, in een kader zoals hieronder:

Omschrijving van de eerste verplichte parameter. Omschrijving van de tweede verplichte parameter. Omschrijving van de eerste optionele parameter. Omschrijving van de tweede optionele parameter.

Voor verplichte parameters geldt dat enkel de waarde van die parameter verplicht moet opgegeven worden en dit in de volgorde waarin die verplichte parameters worden vermeld in het kader (of in de help-functie van R), voor optionele parameters geldt dat wanneer deze niet expliciet worden opgegeven (met de naam van die parameter, een gelijkheidsteken en de gewenste waarde) hun standaardwaarde wordt toegepast. Sommige functies in R hebben tientallen optionele parameters. We zullen hier echter enkel de parameters vermelden die we nodig hebben. In de help-functie van R vind je ongetwijfeld meer details. In R zijn bovendien heel wat functies overloaded; het kan dus voorkomen dat we een bepaalde functie meerdere keren bespreken waarbij de waarde van een parameter de ene keer een getal is, de andere keer een vector, en nog een andere keer een geheel databestand. Het is niet de bedoeling dat je alle functies die in deze tekst worden besproken uit het hoofd zou leren. In de bijlagen vind je een overzicht van de functies die we gebruiken. Anderzijds is er geen reden om je te beperken tot de functies die wij hier voorstellen: “Google” je vragen en de kans is groot dat iemand je voor is geweest met een gelijkaardige vraag die beantwoord is geweest op een forum of mailinglist.

4

2

Basistechnieken

In dit deel gaan we R toepassen op enkele eenvoudige wiskundige voorbeelden zodat je gewend geraakt aan het programma. Pas in het volgende deel zullen we ingaan op het gebruik van R voor statistiek. Verwacht niet dat we hier het hele middelbaar gaan herhalen, wel geven we voldoende basis mee opdat je zelf aan de slag kunt om andere problemen op te lossen. Als je R opstart krijg je het scherm van de R-console te zien en knippert de cursor achter het > teken. Daar kan je instructies typen. Een + vooraan op een lijn geeft aan dat die lijn deel uitmaakt van de input op voorgaande regel(s). Alle overige regels die je te zien krijgt zijn dus output. Je kan de toets <↑> op je toetsenbord gebruiken om eerder getypte regels terug op te vragen. Dit is zeer handig als je een fout hebt getypt. In plaats van alles opnieuw te typen druk je op <↑>, en je krijgt de laatste lijn opnieuw, maar nu kan je deze aanpassen. Zeker wanneer je meer dan n lijn terug nodig zal hebben kan het handig zijn deze vanuit R te kopiëren en in een teksteditor te plakken. Je kan meerdere opdrachtlijnen tegelijk plakken in R; deze worden dan van boven naar onder uitgevoerd.

2.1

Rekentechnieken, variabelen en vectoren

Om te beginnen, typ achter het > teken volgende input (spaties zijn optioneel): > 1 + 1 Gevolgd door <ENTER>. Dit resulteert in onderstaande output: [1] 2 [1] verwijst hier naar het feit dat dit het eerste element is van de output. De eigenlijke output (het resultaat) staat hierachter en is uiteraard 2. Uiteraard kent R ook andere bewerkingen (Merk hier op dat R het punt gebruikt als decimaalteken en niet met een komma.): > 17 / 5 [1] 3.4 R hanteert de gebruikelijke voorrangsregels van bewerkingen:5 > (3 + 4) * 7 + 2 [1] 51 Voor machten, gebruik een caret-teken (ˆ) : > (17-5)^2 [1] 144 Voor wortels gebruiken we machtnotatie, dus

√

1

25 = 25 2 :

> 25^(1/2) [1] 5 We hadden ook 25^.5 kunnen schrijven, maar echter niet 25^1/2 aangezien dat volgens de volgorde van bewerkingen ge¨ınterpreteerd zou worden als 5 Haakjes

→ Macht/Wortel → Product/Deling → Som/Verschil.

5

251 2 .

Breuken geven we in als deling, dus 1 2

+

1 3

1 4

wordt: > ((1/2) + (1/3)) / (1/4) [1] 3.333333 We kunnen een waarde toewijzen aan een variabele (in dit geval de letter x), dit doen we als volgt: > x <- 7 Merk op dat je geen output krijgt. Echter, wel als je de variabele x opvraagt: > x [1] 7 We kunnen nu met x werken alsof het een getal was: > (x + 3) * 2 [1] 20 Dit kan handig zijn wanneer we grote bewerkingen moeten uitvoeren, zoals bijvoorbeeld: 7 4 + − (4 + 8)2 (7 − 4) ∗ 43

17 18

9 ∗

14 3

Als we dit allemaal in één instructieregel ingeven is de kans op fouten groot. Eén ontbrekend haakje en we krijgen geen oplossing. Een mogelijke oplossing: > a > b > c > a [1]

<- 4 / ((4 + 8)^2) <- 7 / ((7 - 4) * 4^3) <- 9 / ((17 / 18) * (14 / 3)) + b - c -1.977781

Naast deze wiskundige operatoren zijn uiteraard ook de logische en booleaanse operatoren beschikbaar: R-code

< <

≤ <=

> >

≥ >=

= ==

6= !=

EN &

OF |

NIET !

Bijvoorbeeld: > (7 > 3) & (5 < 7) [1] TRUE Later zullen we zien dat een variabele in R veel meer kan voorstellen dan alleen getallen, zoals een vector getallen, de uitkomst van een statistische test, een heel databestand,... Het abstraheren zorgt ervoor dat we geen getallen moeten onthouden omdat we die kunnen vervangen door een zelf gekozen variabelenaam. Stel dat we de lengte (in centimeter) van enkele kinderen willen bijhouden en vervolgens het gemiddelde willen berekenen: > jantje_cm <- 107 > marieke_cm <- 103 > fien_cm <- 112 > gemiddelde <- (jantje_cm + marieke_cm + fien_cm) / 3 > gemiddelde [1] 107.3333 Als we echter veel kinderen zouden hebben wordt dit nogal veel typwerk. In plaats van voor ieder kind een variabele aan te maken kunnen we ook één vector aanmaken met als lengte (= aantal elementen) het aantal kinderen: 6


Vector aanmaken. c ( x)

Elementen van de vector gescheiden door een komma.

Stel dat er acht kinderen zijn: > lengte_kinderen <- c(107, 103, 112, 115, 98, 120, 97, 98) Tussen vierkante haken plaats je de index waarvoor je een waarde wil opvragen. Bijvoorbeeld de lengte van het derde kind (Fien): > lengte_kinderen[3] [1] 112

Vector.

Lengte van een vector weergeven. length ( x)


sum ( x)


Sommeer de elementen van een vector.

Vector.

Gebruikmakend van bovenstaande functies kunnen we de gemiddelde lengte berekenen als volgt: > sum(lengte_kinderen) / length(lengte_kinderen) [1] 106.25

Rekenkundig gemiddelde berekenen. mean ( x)


Aangezien R vooral voor statistiek bedoeld is kent de taal uiteraard onderstaande fuctie.

Vector.

> mean(lengte_kinderen) [1] 106.25 Een rekenkundige bewerking uitgevoerd op een vector wordt toegepast op elk element van deze vector: > lengte_in_meter <- lengte_kinderen / 100 > lengte_in_meter [1] 1.07 1.03 1.12 1.15 0.98 1.20 0.97 0.98 We kunnen ook bewerkingen uitvoeren met meerdere vectoren: > gewicht_begin_dieet <- c(110, 83, 98, 88, 79, 94, 80) > gewicht_einde_dieet <- c(102, 84, 95, 86, 79, 93, 76) > gewicht_verloren <- gewicht_begin_dieet - gewicht_einde_dieet > gewicht_verloren [1] 8 -1 3 2 0 1 4

Reeks getallen met een vast interval aanmaken. seq ( from = 1, to = 10, by = 1)

Eerste getal. Laatste getal. Verschil.

7


Stel, we willen een vector creëren die bestaat uit een reeks getallen met een vast interval:

> seq(from = 1, to = 10, by = 1) [1] 1 2 3 4 5 6 7 8 9 10 Aangezien twee van de parameters in bovenstaande functie toch hun standaardwaarde aannemen kan dit ook afgekort worden tot: > seq(to = 10) [1] 1 2 3 4

5

6

7

8

9 10

Een ander voorbeeld: > seq(from = 20, to = 100, by = 10) [1] 20 30 40 50 60 70 80 90 100 Eerder zagen we dat je tussen vierkante haakjes de index kan opgeven van een element uit een vector dat je wilt opvragen. Je kan ook meerdere elementen tegelijk opvragen door in plaats van één index een vector met indices mee te geven. Stel bijvoorbeeld dat we het gewicht aan het begin van het dieet willen opvragen van de even elementen in de vector, dan kan dit zo: > gewicht_begin_dieet[c(2, 4, 6)] [1] 83 88 94 Of, met wat net gezien werd betreffende de functie seq(): > gewicht_begin_dieet[seq(from = 2, to = 6, by = 2)] [1] 83 88 94 Waarbij de tweede methode natuurlijk de voorkeur geniet als je met grotere vectoren werkt. Een waarde kan soms ontbreken, dan spreken we van een “missing”. Stel, we hebben enkele mensen een vragenlijst gegeven waarin een item peilde naar het aantal sekspartners dat met het afgelopen jaar had. Dit is typisch een item dat nogal wat mensen open laten. Voor ontbrekende waarden gebruiken we NA (“not available”). Merk op dat niet iedere functie hier even goed mee overweg kan: > sekspartners <- c(3, 1, NA, NA, 2, 7, NA, 0, 1) > length(sekspartners) [1] 9 > mean(sekspartners) [1] NA Willen we toch het gemiddeld aantal sekspartners op basis van zij die wel geantwoord hebben6 , dan geven we die functie de parameter na.rm = TRUE mee. 7 > mean(sekspartners, na.rm = TRUE) [1] 2.333333

Ga na of een object “missing” is. is.na ( x) > is.na(sekspartners) [1] FALSE FALSE TRUE


We kunnen ook opvragen of een bepaalde waarde wel of niet bestaat:

Object.

TRUE FALSE FALSE

TRUE FALSE FALSE

Variabelen, maar ook elementen in een vector hoeven niet numeriek te zijn (let op de aanhalingstekens, zonder aanhalingstekens zou R op zoek gaan naar een variabele met die naam.): > naam <- "Jantje" > beoordelingen <- c("Zeer slecht", "Slecht", "Neutraal", "Goed", "Zeer goed") 6 Met

de bedenking dat dit misschien een vertekend beeld geeft aangezien het ontbreken hier zeker niet random is. deze parameter op bijna iedere functie toepasbaar is wordt deze niet vermeld in de functiebeschrijvingen.

7 Aangezien

8

Teksteditor openen om een object aan te passen. fix ( x)


Om een variabele aan te passen gebruiken we onderstaande functie:

Naam van het object dat je wil aanpassen.

> fix(beoordelingen) In het kader dat nu verschijnt kan je de optie “Niet van toepassing” toevoegen achter “Zeer goed”. Sluit het kader hierna af, en bevestig.

Help opvragen. help ( topic)


Indien je hulp nodig hebt bij een bepaalde R-functie, dan kan je die als volgt opvragen:

Naam van een functie. Plaats tussen aanhalingstekens indien de parameter speciale tekens bevat.

Oefeningen 2.1.1 a) 1 + 2 ∗ 3 + 4 b) (1 + 2) ∗ (3 + 4) c) 1 + (2 ∗ 3) + 4 √ 2.1.2 4 2401 2.1.3 1+

1 1+

1 2

2.1.4 We vroegen enkele personen naar hun geboortejaar en kregen als antwoorden: 1980, 1973, 1992, 1990, 1987 en 1985. a) Maak een vector voor deze gegevens. b) Maak een nieuwe vector aan op basis van de voorgaande die de leeftijd van deze personen bevat. c) Bereken de gemiddelde leeftijd. 2.1.5 Onderstaande tabel toont de score van 10 kinderen op wiskunde voor- en na extra begeleiding. Maak een vector aan met de scores voor de begeleiding, en één met de scores na de begeleiding. Kind 1 2 3 4 5 6 7 8 9 10 Score voor 74 60 79 72 69 60 78 69 56 63 Score na 75 65 78 74 69 57 80 73 60 72 a) Bereken de gemiddelden van beide vectoren en bereken op basis hiervan de gemiddelde toename. b) Laat R een vector maken die de toename in scores weergeeft. Bereken het gemiddelde van deze vector en vergelijk met a). Verklaar.

2.2

Analyse

Stel, we hebben deze wiskundige functie: f (x) = x2 − 3x

9

Door middel van onderstaande functie geven we aan dat we een wiskundige functie willen aanmaken.8 Hierachter schrijven we het functievoorschrift (dit wijkt dus enigszins af van de andere R-functies waarbij alles tussen de ronde haakjes valt. Wiskundige functie aanmaken. function ( arglist) expr

Variabelen van de functie. Functievoorschrift.

> f <- function(x) x^2 - 3 * x Onze wiskundige functie f kan nu in R opgeroepen worden om bijvoorbeeld functiewaarden te berekenen. Bijvoorbeeld: > f(7) [1] 28 Wat dus hetzelfde is als: > 7^2 - 3 * 7 [1] 28 Een efficiënte manier om de eerste tien veelvouden van 3 te berekenen en deze in te vullen in de functie: > f(seq(from = 3, to = 30, by = 3)) [1] 0 18 54 108 180 270 378 504 648 810

Teken een curve van een wiskundige functie of een puntenwolk. plot ( x, y, from = 0, to = 0, col = "black", xlab = NULL, ylab = NULL, axes = TRUE, add = FALSE)


We kunnen R ook vragen om de functie te tekenen:

X-coördinaten of uitdrukking. Y -coördinaten indien x geen uitdrukking is, anders optioneel. Beginpunt. Eindpunt. Kleur. Label X-as. Label Y -as. Teken assen. Toevoegen aan bestaande grafische voorstelling.

FALSE is een variabele van het type Boolean. De andere mogelijke waarde is uiteraard TRUE. "black" is meestal de standaardwaarde voor kleur, voor andere mogelijkheden, zie de volgende functie: JorisPieters,PeterTheuns,AlainIsaac

Overzicht van beschikbare kleurnamen. colors () > plot(f, from = -10, to = 10)

Dit levert volgende tekening van bijhorende dalparabool op: 8 function()

is dus een functie (in de ruime zin van het woord) die een wiskundige functie cre¨ eert.

10

120 100 80 60 0

20

40

f(x)

−10

−5

0

5

10

x

Sinusfunctie.


De sinus- en de cosinusfunctie kunnen we als volgt bekomen:

sin ( x)

cos ( x)


Cosinusfunctie.

Hoek in radialen.

Hoek in radialen.

Opgelet: pi is geen functie maar een constante in R die ons uiteraard het getal π geeft. > sinusFunctie <- function(x) sin(x) > cosinusFunctie <- function(x) cos(x) > plot(sinusFunctie, from = - 2 * pi, to = 2 * pi, col = "red", xlab = "", ylab = "sinus (rood) en cosinus (blauw)") > plot(cosinusFunctie, from = - 2 * pi, to = 2 * pi, add = TRUE, col="blue")

11

1.0 0.5 0.0 −0.5 −1.0

sinus (rood) en cosinus (blauw)

−6

−4

−2

0

2

4

6

R laat ook toe om de afgeleide van een functie te bepalen. Daarvoor hebben we 2 (functie)onderdelen nodig:

Vergelijking van een functie. body ( fun)

Uitdrukking. Naam van de variabele waarnaar we afleiden (tussen aanhalingstekens). JorisPieters,PeterTheuns,AlainIsaac

deriv ( expr, name)


Afgeleide van een uitdrukking bepalen.

Functie.

Die kunnen dan worden gecombineerd tot: > bodyf <- body(f) > deriv(bodyf, "x") Wat volgende output geeft: expression({ .value <- x^2 - 3 * x .grad <- array(0, c(length(.value), 1L), list(NULL, c("x"))) .grad[, "x"] <- 2 * x - 3 attr(.value, "gradient") <- .grad .value }) De afgeleide van de functie vinden we in deze output achter .grad[, "x"] <- , en inderdaad, de afgeleide van y = x2 − 3x is gelijk aan y 0 = 2x − 3. De bepaalde integraal van een functie vinden we via via onderstaande functie:

12

integrate ( f, lower, upper)


Integraal van een functie berekenen. Functie. Ondergrens. Bovengrens.

Stel dat we onderstaande wiskundige functie wensen te integreren met ondergrens = 1 en bovengrens = 1.5: g(x) =

10 5 + ex2

Exponenti¨ ele functie ex berekenen. exp ( x)


Waarbij ex verkregen kan worden door:

Getal of vector.

> g <- function(x) (10 / (5 + exp(x^2))) > integrate(g, 1, 1.5) 0.5064493 with absolute error < 5.6e-15 Met andere woorden, de oppervlakte tussen de X-as en deze functie tussen 1 en 1.5 op de X-as bedraagt ongeveer 0.5. Grafisch: > plot(g, from = 0, to = 3, xlab = "", ylab = "")

Lijnen tekenen op een grafische voorstelling. lines ( x, y, col = "black") > > > >


Bakenen we het stuk af waarin we ge¨ınteresseerd zijn:

Vector met x-coördinaten. Vector met y-coördinaten Kleur.

lines(c(1, 1), c(0, g(1)), col = "green") lines(c(1.5, 1.5), c(0, g(1.5)), col = "green") lines(c(1, 1.5), c(0, 0), col = "green") lines(c(1, 1.5), c(g(1), g(1.5)), col = "green")

Aangezien g(x) quasi lineair verloopt tussen 1 en 1.5 kunnen we met eenvoudige meetkunde een goede schatting geven van deze oppervlakte. Deel de oppervlakte op in een driehoek en een vierhoek: > lines(c(1, 1.5), c(g(1.5), g(1.5)), col="green") Berekenen van de oppervlakten (respectievelijk basis · hoogte en (basis · hoogte) / 2): > vierhoek <- 0.5 * g(1.5) > driehoek <- (0.5 * (g(1) - g(1.5))) / 2 > vierhoek + driehoek [1] 0.496466

13

1.5 1.0 0.5 0.0 0.0

0.5

1.0

1.5

2.0

2.5

3.0

Dit resultaat ligt zeer dicht bij de berekening door middel van de integraal. Aan de hand van de tekening kan je zien waarom de schatting iets lager uitvalt dan de berekening. Uiteraard kunnen we dergelijke trukjes eigenlijk alleen toepassen op (bij benadering) lineaire stukken in een functie, net daarom zijn er bepaalde integralen.

Grafische voorstelling van een veelterm weergeven. polygon ( x, y, col = "black")


Om een oppervlakte volledig op te vullen hebben we een veelhoek nodig:

Vector met x-coördinaten. Vector met y-coördinaten. Kleur.

Willen we bijvoorbeeld bovenstaande oppervlakte exact gevuld zien (en niet benaderd door een driehoek en een rechthoek), dan hebben we een polygoon nodig met volgende punten: x 1.0 1.0 1.0 . . . 1 1.0 . . . 2 .. .

y 0 g(1.0) g(1.0 . . . 1) g(1.0 . . . 2) .. .

1.5 1.5

g(1.5) 0

Met andere woorden: het punt links onder, een hele reeks punten op de curve, en het punt rechts onder. Als we die reeks punten om de 0.01 zetten zal het resultaat nauwkeurig genoeg zijn voor een mooie figuur: > > > > >

plot(g, from = 0, to = 3, xlab = "", ylab = "") r <- seq(from = 1, to = 1.5, by = 0.01) x <- c(1, r, 1.5) y <- c(0, g(r), 0) polygon(x,y, col = "blue")

14

1.5 1.0 0.5 0.0 0.0

0.5

1.0

1.5

2.0

2.5

3.0

Oefeningen 2.2.1 Gegeven: f (x) = 2 − (x − 2)2 Integreer deze functie met ondergrens 1 en bovengrens 3. Stel grafisch voor tussen 0 en 4, en probeer de oppervlakte ook te schatten.

2.3

Matrixrekenen

Gegeven zijn twee vectoren:     7 1 9 2   b= a= −1 3 3 4 > a <- c(1, 2, 3, 4) > b <- c(7, 9, -1, 3) Eerder hebben we al gezien dat we de som of het verschil tussen twee vectoren van gelijke lengte kunnen maken: > a + b [1] 8 11

2

7

Het scalaire veelvoud van een vector is de vector die je bekomt door elk element van die vecor te vermenigvuldigen met hetzelfde getal: > 3 * a [1] 3 6

9 12

Het (matrix)product van een rijvector met een kolomvector geeft als uitkomst een getal. In R wordt het matricieel product voorgesteld door %*%, en om een matrix te transponeren (aT ) gebruiken we onderstaande:

15


Transponeren van een vector of matrix. t

Vector of matrix.

( x) > t(a) %*% b [,1] [1,] 34

Wanneer we echter een kolomvector vermenigvuldigen met een rijvector krijgen we een matrix: > a %*% t(b) [,1] [,2] [,3] [,4] [1,] 7 9 -1 3 [2,] 14 18 -2 6 [3,] 21 27 -3 9 [4,] 28 36 -4 12 Tenslotte zijn een kolom- en een rijvector niets anders dan matrices met respectievelijk één kolom of één rij. Het product van een matrix A van formaat (m × n) (m rijen en n kolommen) met een matrix B van formaat (n × p) levert een matrix C van formaat (m × p) 9 : A(m×n) × B(n×p) = C(m×p) Laten we dat eens proberen met: 

A=

2 7

5 1

4 2

3 4

1 8

1 2  B= 3 4 3

 7 5  3  8 6

Geef de waarden in als één vector met de kolommen achter elkaar, en maak hiervan een matrix door middel van:

matrix ( data, nrow = 1, ncol = 1)


Matrix cre¨ eren op basis van een vector.

Vector. Aantal rijen. Aantal kolommen.

> A <- matrix(c(2, 7, 5, 1, 4, 2, 3, 4, 1, 8), nrow = 2, ncol = 5) > A [,1] [,2] [,3] [,4] [,5] [1,] 2 5 4 3 1 [2,] 7 1 2 4 8 > B <- matrix(c(1, 2, 3, 4, 3, 7, 5, 3, 8, 6), nrow = 5, ncol = 2) > B [,1] [,2] [1,] 1 7 [2,] 2 5 [3,] 3 3 [4,] 4 8 [5,] 3 6 > A %*% B 9 Om

een product van 2 matrices mogelijk te maken moet het aantal kolommen van de eerste matrix gelijk zijn aan het aantal rijen van de tweede matrix.

16

[,1] [,2] [1,] 39 81 [2,] 55 140 > B %*% A [,1] [,2] [,3] [,4] [,5] [1,] 51 12 18 31 57 [2,] 39 15 18 26 42 [3,] 27 18 18 21 27 [4,] 64 28 32 44 68 [5,] 48 21 24 33 51 Indien we enkel het element op rij 2, kolom 5 willen opvragen uit A: > A[2, 5] [1] 8 Stel, we hebben  4 M = 8 6

onderstaande vierkante matrix M, en we willen daar de inverse van bekomen:  7 1 5 4 9 2

> M <- matrix(c(4, 8, 6, 7, 5, 9, 1, 4, 2), nrow = 3, ncol = 3) > M [,1] [,2] [,3] [1,] 4 7 1 [2,] 8 5 4 [3,] 6 9 2

Inverse van een matrix berekenen. solve ( a)


Voor het berekenen van de inverse M−1 gebruiken we:

Vierkante matrix.

En niet ^(-1), dat neemt de inverse van de matrixelementen: > invM <- solve(M) > invM [,1] [,2] [,3] [1,] 4.333333 0.8333333 -3.833333 [2,] -1.333333 -0.3333333 1.333333 [3,] -7.000000 -1.0000000 6.000000 Hoe weten we nu zeker dat dit de inverse is van M? Vermenigvuldiging met het origineel moet uiteraard de eenheidsmatrix opleveren: > M %*% invM [,1] [,2] [,3] [1,] 1.000000e+00 -5.551115e-16 0 [2,] -3.552714e-15 1.000000e+00 0 [3,] -3.552714e-15 -8.881784e-16 1

Afronden.


Fout? Nee, gewoon een kwestie van afronding:

round ( x, digits = 0)

Getal of vector. Aantal decimalen.

17

> round(M %*% invM) [,1] [,2] [,3] [1,] 1 0 0 [2,] 0 1 0 [3,] 0 0 1 Oefeningen 2.3.1 

18 17  X= 18 19 18

1 1 0 0 1

Bereken XT X

−1

 16 8  4  12 12

18

3

Beschrijvende statistiek

3.1

E´ endimensionaal

Stel dat we onderstaande gegevens hebben van twintig eerstejaarsstudenten: “ID”: Het nummer van de student, om ze gemakkelijk te kunnen identificiëren. “Geslacht” “Hoe goed in wiskunde”: Hierbij hebben we de student zelf gevraagd een inschatting te geven van haar of zijn wiskundige capaciteiten, op een schaal van “Heel slecht” tot “Heel goed”. “Leeftijd” “Lengte”

ID 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

Geslacht Man Man Vrouw Vrouw Vrouw Vrouw Man Vrouw Vrouw Man Man Vrouw Man Vrouw Vrouw Vrouw Vrouw Vrouw Vrouw Man

Hoe goed in wiskunde Heel slecht Neutraal Slecht Goed Heel slecht Heel slecht Neutraal Goed Heel slecht Neutraal Goed Slecht Goed Heel goed Slecht Slecht Goed Neutraal Neutraal Heel goed

Leeftijd 18 19 18 18 17 18 18 19 18 21 18 46 18 19 20 18 18 22 18 21

Lengte 171.3 169.8 158.7 164.9 173.5 169.9 168.3 169.4 162.7 180.5 172.5 174.9 164.7 156.0 159.3 169.0 171.1 175.4 168.0 170.2

Merk op dat er één variabele op nominale schaal is, één op ordinale schaal, en twee die op minstens intervalschaal gemeten zijn. Omdat voor het merendeel van de statistische analyses het onderscheid tussen intervalschaal en ratioschaal niet cruciaal is zullen we deze hier niet apart bespreken. De reden dat we wel twee variabelen op minstens intervalschaal zullen bespreken betreft een ander onderscheid. De ene, namelijk leeftijd, is gemeten op een discrete schaal. We hebben iedere student gevraagd hoe oud hij of zij is, en hoewel leeftijd onderliggend uiteraard continu is (iemand is enkel precies 18 jaar op de dag en op het exacte moment van zijn of haar geboorte, 18jaar later ...) zal een bevraging toch een discrete leeftijd opleveren (niemand zal antwoorden dat hij of zij 18.394851314... jaar oud is). Lengte daarentegen zullen we hier behandelen alsof die werd gemeten op een continue schaal10 . Het zou kort voor de bocht zijn om te zeggen dat discreet steeds overeenkomt met afwezigheid van decimalen, terwijl continu steeds overeen zou komen met waarden waarbij decimalen betrokken zijn. Echter, in de praktijk komt het vaak hier op neer.11 10 Strikt

genomen zou je kunnen zeggen dat ook lengte op een discrete schaal werd gemeten, weliswaar met kleine intervallen (1mm), maar voor echte continue gegevens zouden we cijfers tot op een oneindig aantal decimalen moeten beschouwen. Hieruit volgt dan ook dat strikt genomen data wellicht nooit “echt”continu kunnen worden geregistreerd. 11 Een mooi voorbeeld is geld: Het aantal euro dat je op zak hebt heeft decimalen (je hebt bijvoorbeeld 12.37 euro op zak), maar toch is dit discreet aangezien je een eurocent fysiek niet verder kan opdelen. In financi¨ ele instellingen daarentegen zal men geld toch (in sterkere mate) continu beschouwen door met meer decimalen te werken. Dit gebeurt omdat anders kleine afwijkingen door (foutieve) afrondingen kunnen resulteren in zeer grote totale afwijkingen na duizenden transacties (De Vancouver Stock Exchange maakte ooit deze fout en stond na twee jaar op slechts de helft van de eigenlijke waarde).

19

Databestand openen.


We gaan deze data niet intypen aangezien er reeds een databestand beschikbaar is. Om zogenaamde “platte tekst” zoals het databestand “studenten.txt” in te lezen gebruiken we deze functie:

read.table ( file, header = FALSE)

Bestandsnaam (tussen aanhalingstekens). Eerste lijn bevat naam variabele

Huidige werkmap tonen.


Om niet iedere keer het volledige pad naar het bestand te moeten intypen stellen we een werkmap in. Op die manier weet R in welke map het bestand gevonden kan worden. Via File → Change dir... stel je de werkmap in op de map waar je “studenten.txt” hebt opgeslagen (Op Mac: R → Preferences → Startup → Directory → Change). Wil je weten wat de huidige werkmap is, typ dan:

getwd () Nu kan je het databestand openen: > studenten <- read.table("studenten.txt", header = TRUE) > studenten ID Geslacht Wiskunde Leeftijd Lengte 1 1 Man 1 18 171.3 2 2 Man 3 19 169.8 3 3 Vrouw 2 18 158.7 . . . . . . . . . . . . . . . . . . 20 20 Man 5 21 170.2 Er zijn twee opvallende verschillen tussen de tabel in de tekst en de dataset. Ten eerste hebben we “Hoe goed in wiskunde” veranderd in “Wiskunde” omdat een naam van een variabele voor R geen spaties mag bevatten. Als alternatief zou je ook de spaties kunnen weglaten, of ze bijvoorbeeld vervangen door een laag streepje, zoals we eerder in de tekst hebben gedaan. Echter, iedere keer Hoe goed in wiskunde intypen is natuurlijk tijdverspilling, vandaar onze voorkeur om de naam van de variabele te vereenvoudigen. Ten tweede werden de waarden “Heel slecht” tot “Heel goed” vervangen door de cijfers 1 tot en met 5. Het betreft hier een ordinale schaal, maar R kan (uiteraard) de volgorde van de waarden niet kennen. Bij gebrek aan beter zou R alfabetische volgorde gebruiken, maar dat zou hier niet correct zijn. Uiteindelijk is een datatabel eigenlijk niets anders dan een vector van variabelen (die op zichzelf natuurlijk ook vectoren kunnen zijn). R biedt een grote hoeveelheid opties om verschillende soorten data-bestanden te openen. Dit is handig aangezien er heel wat verschillende manieren zijn waarop een databestand kan zijn aangemaakt afhankelijk van het gebruikte programma, wie het databestand heeft gemaakt, enzoverder. Opdat we ons hier kunnen focussen op statistiek in R zijn alle databestanden bij deze cursus op dezelfde manier opgemaakt. In een latere cursus zal je ook leren met andere soorten bestanden om te gaan. Als je de databestanden die horen bij deze cursus opent vanuit je besturingssysteem dan zullen die normaal gezien in een teksteditor openen omdat ze de extentie .txt hebben. Je kan zien dat de eerste rij telkens de namen van de variabele(n) bevat, daarom dat we header = TRUE typten. Verder zie je dat iedere volgende rij overeenkomt met één case. De waarden zijn van elkaar gescheiden door een (dat is niet hetzelfde als meerdere spaties!), en dat het decimaalteken het punt is en niet de komma. In andere situaties zou je read.table() kunnen aanvullen met de parameters sep of dec die respectievelijk aangeven hoe de verschillende waarden scheiden zijn, en wat het decimaalteken is. De standaardwaarden hiervan komen echter al overeen met hoe onze databestanden zijn opgesteld, zodat we gebruik kunnen maken van de standaardwaarden (en ze dus niet moeten specifiëren). Om zelf databestanden aan te maken kan je voorlopig best gebruik maken van een teksteditor, al dan niet

20

in combinatie met een spreadsheet-programma. Tonen hoe je dit in R doet zou ons momenteel te ver afleiden van de doelstellingen van deze cursus (maar ook dit zal in een latere cursus goedgemaakt worden). Variabele op nominale schaal Door middel van het dollar-teken kan je één vector uit de data kiezen: > geslacht <- studenten$Geslacht > geslacht [1] Man Man Vrouw Vrouw Vrouw Vrouw Man Vrouw Vrouw Man [11] Man Vrouw Man Vrouw Vrouw Vrouw Vrouw Vrouw Vrouw Man Levels: Man Vrouw Een frequentietabel kunnen we bekomen door middel van deze functie: JorisPieters,PeterTheuns,AlainIsaac

Frequentietabel weergeven. table ( ...)

Eén of meer objecten.

> table(geslacht) geslacht Man Vrouw 7 13 Voor wie liever relatieve frequenties ziet: > table(geslacht) / length(geslacht) geslacht Man Vrouw 0.35 0.65

Barplot tonen.


Een kolomdiagram kunnen we bekomen door middel van:

barplot ( height)

Vector met waarden die de hoogte van de balkjes bevat.

Merk op dat we de frequentietabel van de variabele moeten opgeven, niet de variabele zelf:

0

2

4

6

8

10

12

> barplot(table(geslacht))

Man

Vrouw

21

Variabele op ordinale schaal > wiskunde <- studenten$Wiskunde > wiskunde [1] 1 3 2 4 1 1 3 4 1 3 4 2 4 5 2 2 4 3 3 5

Cumulatieve frequenties berekenen. cumsum ( x)


Vanaf ordinale schaal kunnen we cumulatieve frequenties opvragen:

Vector.

> cumsum(table(wiskunde)) 1 2 3 4 5 4 8 13 18 20 Relatieve cumulatieve frequenties opvragen: > cumsum(table(wiskunde)) / length(wiskunde) 1 2 3 4 5 0.20 0.40 0.65 0.90 1.00 Wat uiteraard ook kan is (let op het verschil in haakjes in vergelijking met voorgaande en ga na waarom beide exact hetzelfde resultaat opleveren): > cumsum(table(wiskunde) / length(wiskunde)) 1 2 3 4 5 0.20 0.40 0.65 0.90 1.00

Vijf-getallen-samenvatting weergeven (frequentietabel voor schaal zonder cijfers). summary ( object)


Stel, we willen de vijf-getallen-samenvatting:

Een object waarvoor een samenvatting gemaakt moet worden.

> summary(wiskunde) Min. 1st Qu. Median 1.00 2.00 3.00

Mean 3rd Qu. 2.85 4.00

Max. 5.00

Merk op dat R ook het gemiddelde geeft terwijl dit hier een ordinale schaal is. Dat is uiteraard klinkklare onzin! Het is steeds aan jou om na te denken over welke resultaten je kan of mag interpreteren en hoe. De software doet enkel de berekeningen.

Boxplot tonen. boxplot ( x, names = NULL)


Probeer de vijf getallen terug te vinden in de boxplot die je genereert door middel van:

Data waarop de boxplot is gebaseerd. Vector met namen voor verschillende boxplots in één afbeelding.

> boxplot(wiskunde)

22

5 4 3 2 1

Variabele op interval- of ratio schaal: Discreet > leeftijd <- studenten$Leeftijd > leeftijd [1] 18 19 18 18 17 18 18 19 18 21 18 46 18 19 20 18 18 22 18 21 De frequentietabel: > table(leeftijd) leeftijd 17 18 19 20 21 22 46 1 11 3 1 2 1 1 > summary(leeftijd) Min. 1st Qu. Median 17.00 18.00 18.00

Mean 3rd Qu. 20.10 19.25

Max. 46.00

Hoewel de boxplot vooral voor variabelen op oridinale schaal populair is mag het zeker ook op een hogere meetschaal:

20

25

30

35

40

45

> boxplot(leeftijd)

Merk op dat we twee outliers hebben. R maakt geen onderscheid tussen verschillende soorten outliers.12 12 Sommige boeken en software (waaronder SPSS) geven enkel outliers tussen 1.5 en 3 interkwartielafstanden van de box aan met een bolletje en zullen extremere outliers aanduiden met een sterretje.

23

Eerder hadden we al gezien hoe het gemiddelde op te vragen:

n

x ¯=

1X xi n i=1

> n <- length(leeftijd) > sum(leeftijd) / n [1] 20.1 Ofwel: > mean(leeftijd) [1] 20.1 De formule voor de variantie (als beschrijving van de steekproef) is:

n

s2 =

1X 2 (xi − x ¯) n i=1

Met de hand hebben we dit berekend door eerst de verschilscores ten opzichte van het gemiddelde te berekenen: > verschilscores <- leeftijd - mean(leeftijd) > verschilscores [1] -2.1 -1.1 -2.1 -2.1 -3.1 -2.1 -2.1 -1.1 -2.1 [17] -2.1 1.9 -2.1 0.9

0.9 -2.1 25.9 -2.1 -1.1 -0.1 -2.1

Ter herinnering, de som van deze getallen zal steeds gelijk zijn aan nul. Echter, het kwadraat van deze verschilscores is wel bruikbaar: > kwadrverschil <- verschilscores^2 > kwadrverschil [1] 4.41 1.21 4.41 4.41 9.61 [12] 670.81 4.41 1.21 0.01 4.41

4.41 4.41

4.41 3.61

1.21 4.41

4.41 0.81

0.81

4.41

De variantie is net het gemiddelde van gekwadrateerde afwijkingen ten opzichte van het gemiddelde, dus: > variantie <- mean(kwadrverschil) > variantie [1] 36.89 Uiteraard kan dit ook in één keer: > varileeftijd <- mean((leeftijd - mean(leeftijd))^2) > varileeftijd [1] 36.89

Bereken de vierkantswortel. sqrt ( x)


En de standaardafwijking is uiteraard de vierkantswortel van de variantie:

Getal of vector.

> stdafwleeftijd <- sqrt(varileeftijd) > stdafwleeftijd [1] 6.073714

24

De variantie kan ook op een andere manier berekend worden: n

s2

=

1X 2 (xi − x ¯) n i=1

=

1X 2 xi − 2 xi x ¯+x ¯2 n i=1

=

1X 1X 2 1X xi − 2 xi x ¯+ x ¯ n i=1 n i=1 n i=1

=

1 X 2 2x ¯X xi − xi + x ¯ n i=1 n i=1

=

1X 2 x − 2x ¯2 + x ¯ n i=1 i

=

1X 2 x −x ¯2 n i=1 i

n

n

n

n

n

n

n

n

Laten we dat eens uittesten (let goed op hoe de haakjes staan!): > mean(leeftijd^2) - mean(leeftijd)^2 [1] 36.89 Uiteraard kent R een functie om de variantie en de standaardafwijking te berekenen:

var ( x)


Populatievariantie berekenen (schatter!).

Populatiestandaardafwijking berekenen (schatter!). sd ( x)


Vector.

Vector.

> var(leeftijd) [1] 38.83158 > sd(leeftijd) [1] 6.231499 Dit resultaat is niet hetzelfde als wat we hierboven vonden, verklaar dit! R berekent de schatter van de variantie in de populatie (delen door (n − 1) in plaats van door n): > var(leeftijd) * (n - 1) / n [1] 36.89 We kunnen ook de andere centrale momenten berekenen:13 n

mk =

1X k (xi − x ¯) n i=1

> m3 <- mean((leeftijd - mean(leeftijd))^3) > m3 [1] 862.332 > m4 <- mean((leeftijd - mean(leeftijd))^4) 13 Deze formule is weer voor de beschrijving van de steekproef. Voor schatting van een populatieparameter moet je delen door (n − 1) in plaats van door n

25

> m4 [1] 22515.55 Deze kunnen we vervolgens gebruiken om de parameters van symmetrie (g1 ) en gepiektheid (g2 ) te berekenen: > g1 <- m3/stdafwleeftijd^3 > g1 [1] 3.848677 > g2 <- m4/stdafwleeftijd^4 - 3 > g2 [1] 13.54494

g1 =

m3 = 3.848677 s3

g2 =

m4 − 3 = 13.54494 s4

Aangezien g1 groter is dan nul is de verdeling wellicht rechts scheef, terwijl g2 groter dan nul aangeeft dat de verdeling leptokurtisch is. Een grafische voorstelling zou hier uiteraard van pas komen. Het standaardhistogram in R blijkt echter niet geschikt te zijn wanneer we slechts een beperkt aantal mogelijke waarden hebben (we hebben “leeftijd” hier namelijk niet continu gemeten). We zullen daarom terug de barplot() gebruiken, maar aangepast met een bijkomende parameter. Gebruik best deze methode telkens wanneer je variabelen op interval- of ratioschaal hebt die discreet gemeten zijn: JorisPieters,PeterTheuns,AlainIsaac

Barplot weergeven. barplot ( height, space = NULL)

Vector met waarden die de hoogte van de balkjes bevat. Spatie tussen de balkjes (standaard wordt een spatie gebruik, wat correct is op nominale en ordinale schaal).

Frequentietabel met frequentie nul voor ontbrekende meetwaarden. factor ( x, levels)


Eerst vullen we de variabele “leeftijd” aan met frequenties gelijk aan nul voor alle mogelijke discrete waarden tussen de laagste en hoogste waargenomen leeftijd. Dit doen we door middel van:

Vector. Vector met alle mogelijke waarden.

> longleeftijd <- factor(leeftijd, seq(from = min(leeftijd), to = max(leeftijd)))

Laagste waarde.


Met vanzelfsprekend:

min ( ...) JorisPieters,PeterTheuns,AlainIsaac

Hoogste waarde.

Vector.

max ( ...)

Vector.

De frequentietabel ziet er nu zo uit: > table(longleeftijd) longleeftijd

26

17 18 19 20 21 22 23 24 ... 45 46 1 11 3 1 2 1 0 0 ... 0 1

0

2

4

6

8

10

> barplot(table(longleeftijd), space = 0)

17

19

21

23

25

27

29

31

33

35

37

39

41

43

45

Maak een stamdiagram. stem ( x)


Het stamdiagram:

Vector.

> stem(leeftijd) The decimal point is 1 digit(s) to the right of the | 1 2 3 4

| 788888888888999 | 0112 | | 6

Variabele op interval- of ratioschaal: Continu > lengte <- studenten$Lengte > lengte [1] 171.3 169.8 158.7 164.9 173.5 169.9 168.3 169.4 162.7 180.5 [11] 172.5 174.9 164.7 156.0 159.3 169.0 171.1 175.4 168.0 170.2 Al wat je kon toepassen op een discrete interval- of ratioschaal kan je ook hier toepassen. We bespreken data van een continue meting echter apart omdat het histogram hier anders wordt geproduceerd dan hierboven. Een histogram waarbij we R toelaten om de gegevens in klassen te delen (wat eigenlijk altijd een vereiste is voor histogrammen bij continuë data) bekomen we met:

27


Histogram tonen. hist ( x, freq = TRUE, nclass = NULL)

Vector. Absolute frequenties? (anders relatieve) Aantal klassen; indien weggelaten wordt automatisch gekozen.

> hist(lengte)

3 0

1

2

Frequency

4

5

6

Histogram of lengte

155

160

165

170

175

180

185

lengte

Merk op dat R zelf klassen kiest. Meestal is deze keuze vrij goed. Er is een parameter om dit bij te sturen, maar als je niet tevreden bent met het resultaat kan je wellicht beter de methode voor de gemodificeerde barplot die eerder werd besproken (één balkje per waarde) toepassen. Oefeningen 3.1.1 Het bestand “neerslag.txt” bevat de jaargemiddelden van de hoeveelheid neerslag (mm per m2 ) van 1833 tot en met 2011 in België (Bron: Vlaamse Milieumaatschappij). Bereken gemiddelde, variantie, symmetrie (g1 ) en gepiektheid (g2 ). Gebruik steeds de formules die de steekproef beschrijven. Maak een histogram en ga na of je de berekende parameters kan relateren aan wat je ziet in het histogram. 3.1.2 Het bestand “les.txt” bevat de beoordeling van een groep studenten betreffende de leskwaliteit enerzijds, en hoe interessant men het vak inhoudelijk vond anderzijds. Beide werden gescoord op een schaal van 0 (“Heel erg slecht”) tot en met 10 (“Heel erg goed”). a) Vergelijk de boxplots voor de twee beoordelingen. Bespreek de opvallende verschillen. b) Vergelijk de beoordeling van de leskwaliteit tussen de mannelijke en vrouwelijke studenten. Tips: boxplot() laat je ook toe meerdere variabelen tegelijk op te geven, en zet deze dan naast elkaar: boxplot(vector 1, vector 2, ...). Voor de duidelijkheid kan je de optionele parameter names meegeven wat een vector is met de respectievelijke namen van de variabelen zodat duidelijk is welke boxplot waarbij hoort. Om slechts een deel van een vector op te vragen (in dit geval de leskwaliteit voor enkel de mannen, en voor enkel de vrouwen) kan je tussen vierkante haken de gevraagde index of een vector van indices opgeven, zoals besproken werd in het deel “Basistechnieken”. Je kan echter ook een voorwaarde plaatsen tussen de vierkante haakjes, zoals vector 1[vector 2 == waarde].14

3.1.3 “kleur.txt” bevat de lievelingskleur van kinderen uit een kleuterschool. 14 Let op het dubbel gelijkheidsteken, dit is (zoals in wel meer computertalen) het symbool voor gelijkheid, terwijl een enkel gelijkheidsteken staat voor toewijzing (net als het in R populairdere <-).

28

a) Hoeveel kinderen zijn er in totaal? b) Wat is de meestvoorkomende lievelingskleur? Hoeveel kinderen verkiezen rood? c) Kan je de resultaten uit vraag b) ook terugvinden indien je een kolomdiagram maakt? 3.1.4 “voltijdsdeeltijds.txt” bevat twee variabelen uit een personeelsbestand, namelijk of iemand voltijds of deeltijds werkt voor het bedrijf, en hoeveel uren hij of zij per week werkt. a) Geef het interkwartiel van het aantal uren dat men wekelijks werkt. b) Maak een histogram van het aantal uren dat per week gewerkt wordt. We hebben twee methoden gezien, welke is hier het meest van toepassing en waarom? c) Maak een histogram van enkel die mensen die voltijds werken. Wat kan je concluderen betreffende de symmetrie? Kan je dit staven met behulp van de parameter g1 ? Hoe kan je dit verklaren? 3.1.5 We scoren alfabetisme als volgt: 0: Kan helemaal niet lezen of schrijven (met uitzondering van de eigen naam) 1: Kan losse woorden lezen en schrijven. Moeite met volzinnen. 2: Kan zinnen lezen en schrijven. Moeite met teksten. 3: Kan lezen en schrijven, maar beduidend minder snel dan personen met score 4. 4: Kan perfect lezen en schrijven.

Onderstaande tabellen geven het aantal 18-jarigen aan per categorie voor een steekproef genomen in 1900 en één in 2000. Aantal Aantal Analfabetisme in 1900 in 2000 0 136 5 1 214 12 2 183 43 3 92 213 4 36 873

Repliceer een object een welbepaald aantal keer. rep ( x, times = 1)


Bespreek eventuele verschillen in ligging en spreiding tussen 1900 en 2000. Probeer te verklaren waar eventuele verschillen vandaan kunnen komen. Tip:

Object. Aantal keren herhalen.

Dit is handig bij het ingeven van data. Heb je bijvoorbeeld een vector nodig met eerst tien keer de waarde 0, en vervolgens vijf keer de waarde 1, dan kan dit zo: waarden <- c(rep(0, 10), rep(1, 5)) In plaats van zo: waarden <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1)

Random een bepaald aantal getallen uit een vector kiezen. sample ( x, size, replace = FALSE)


3.1.6

Vector met elementen om uit te kiezen. Aantal trekkingen. Met teruglegging.

29

Zo kunnen we bijvoorbeeld 10 keer gooien met een dobbelsteen simuleren door middel van: dobbel <- sample(seq(from = 1, to = 6), size = 10, replace = TRUE) > dobbel [1] 4 1 1 2 6 6 5 4 4 3 Waarbij jouw resultaat natuurlijk niet hetzelfde zal zijn als wat hier staat. Herhaal dit meerdere keren, wat gebeurt er met het resultaat? Verhoog het aantal keren dat je gooit (100, 1000, ...) en maak telkens een histogram. Wat verandert er aan het histogram wanneer je meer keren gooit? 3.1.7 Je zou nu moeten weten wat R doet als je onderstaande opdracht geeft: > data <- rep(seq(1, 10), 50) a) Ga na of wat je veronderstelde correct is door middel van > table(data) b) Bereken de parameters van symmetrie (g1 ) en gepiektheid (g2 ). c) Kan je de conclusies van g1 en g2 ook terugvinden in een histogram? 3.1.8 Maak voor ieder van onderstaande een vector van tien getallen zodat aan de voorwaarde voldaan is. Er zijn uiteraard meerdere mogelijke antwoorden, de oplossingen zijn dus slecht één mogelijkheid. a) Mediaan is 7. b) Gemiddelde is 38. c) Derde kwartiel is 101. d) Variantie is 60.

3.2

Relaties tussen 2 variabelen

Associatiematen voor contingentietabellen Bij de Newfoundland hond wordt de neus soms gecategoriseerd als smal, gemiddeld of breed. We willen weten of we het geslacht van de hond beter kunnen voorspellen als we weten of de neus smal, gemiddeld of breed is. Reu Teef

Smal 50 72 122

Gemiddeld 68 39 107

Breed 111 18 129

229 129 358

Lambda (λ) geeft aan met welke proportie het aantal (#) foute voorspellingen van de waarde van de ene variabele afneemt, indien de waarde van de andere variabele bekend zou zijn. Met andere woorden, als je de contingentietabel voor 2 variabelen A en B kent, en je moet de waarde van A voorpellen, dan geeft λ aan met welke proportie het aantal foute voorspellingen over de waarde van A daalt indien de waarde van B wordt gegeven. λ=

# foute voorspellingen van A indien B onbekend − # foute voorspellingen van A indien B bekend # foute voorspellingen van A indien B onbekend

Zonder de neus te hebben gezien kunnen we bij het voorspellen van het geslacht best altijd “reu” voorspellen (omdat er in deze tabel meer reuen zijn dan teven). In 229 van de 358 gevallen is onze voorspelling van het geslacht dan correct, we maken 129 foute voorspellingen. Wanneer we echter de neus hebben gezien, dan zullen we op basis van deze gegevens bij een smalle neus best steeds “teef” voorspellen (dat is 72 keer correct, 50 keer fout), bij een middelmatige neus voorspellen we “reu” (68 keer correct, 39 fouten) en bij een brede neus voorspellen we ook best “reu” (111 keer correct, 18

30

fouten). λ=

129 − (50 + 39 + 18) = 0.1705 129

Maak een vector voor elke rij, met daarin de frequenties:

Vectoren als rijen aan elkaar binden tot een datatabel. rbind ( ...)


> reu <- c(50, 68, 111) > teef <- c(72, 39, 18)

Vectoren gescheiden door een komma.

> tabel <- rbind(reu, teef) > tabel [,1] [,2] [,3] reu 50 68 111 teef 72 39 18 Indien we willen dat in R ook de kolomnamen getoond worden:

colnames ( x)


Voeg kolomnamen toe aan een tabel. Tabel.

> colnames(tabel) <- c("smal", "gemiddeld", "breed") > tabel smal gemiddeld breed reu 50 68 111 teef 72 39 18

Vectoren als colommen aan elkaar binden tot een datatabel. Vectoren gescheiden door een komma.

Voeg rijnamen toe aan een tabel. rownames ( x)


cbind ( ...)


Merk op dat we deze tabel even goed als volgt hadden kunnen aanmaken met behulp van onderstaande functies:

Tabel.

> > > > > >

smal <- c(50, 72) gemiddeld <- c(68, 39) breed <- c(111, 18) tabel <- cbind(smal, gemiddeld, breed) rownames(tabel) <- c("reu", "teef") tabel smal gemiddeld breed reu 50 68 111 teef 72 39 18 Om R een λ-coëfficiënt te laten berekenen heb je het package rapport nodig, dit moet je eerst even ophalen met de Packages Installer en vervolgens laden (Zie appendix “Packages” voor meer uitleg.):

31


Package laden. library ( package)

Naam van het package

Je kan nu de package laden door: > library(rapport) Je zal merken dat automatisch ook enkele andere packages worden geladen. Vervolgens kan je λ opvragen met volgende functie:

lambda.test ( table)


Lambda (λ) berekenen (package: rapport)

Een tabel met twee variabelen.

> lambda.test(tabel) $row [1] 0.2358079 $col [1] 0.1705426 We krijgen twee resulaten: onder row vinden we λ waarbij de rij-variabele (geslacht van de hond) de verklarende variabele is, dat is niet de λ die we zoeken. Onder col vinden we λ waarbij de kolomvariabele (smalle, gemiddelde of brede neus) de verklarende variabele is, en dat is wel de λ waarin we ge¨ınteresseerd zijn. Wanneer we het geslacht van een hond uit de tabel moeten voorspellen maken we zo’n 17% minder fouten indien we weten of die hond een smalle, gemiddelde of brede neus heeft dan wanneer we dit niet weten. De Q-coëfficiënt is een associatiemaat die enkel van toepassing is op 2 × 2-tabellen. Voor een contingentietabel met celfrequenties gelijk aan A, B, C en D zoals deze: A C

B D

is Q gedefiniëerd als: Q=

(A × D) − (B × C) (A × D) + (B × C)

Q=

product diagonaal 1 − product diagonaal 2 product diagonaal 1 + product diagonaal 2

of

Verwijder de namen of labels van een object. unname ( obj) > > > >

a b c d

<<<<-

tabel[1, tabel[1, tabel[2, tabel[2,

Object. 1] 3] 1] 3] 32


We laten nu even de honden met een gemiddeld brede neus buiten beschouwing. Onderstaande functie gebruiken we om de labels - zoals “reu” of “teef” - weg te nemen van de celfrequentie, zodat we nog enkel de getallen overhouden, zonder de bijhorende rij-labels.):

> Q <- unname((a * d - b * c) / (a * d + b * c)) > Q [1] -0.7975709 De absolute waarde van Q zit dicht bij één, wat indicatief is voor een sterk verband tussen de breedte van de neus en het geslacht bij Newfoundlanders (gemiddeld brede neuzen buiten beschouwing gelaten). Een andere associatiemaat is de χ2 : χ2 =

r X k 2 X (Fij (o) − Fij (e)) i=1 j=1

Fij (e)

Bereken de χ2 contingentie. chisq.test ( x)


Waarbij Fij (o) de geobserveerde frequentie van een cel is, en Fij (e) de verwachte frequentie op basis van de marginale frequenties (indien er geen verband zou zijn), terwijl r en k zijn respectievelijk het aantal rijen en kolommen van de tabel. In R kunnen we de test uitvoeren door middel van:

Tweedimensionale tabel (matrix).

> chisq.test(tabel) Pearson’s Chi-squared test data: tabel X-squared = 55.2516, df = 2, p-value = 1.005e-12 De χ2 -waarde bedraagt 55.2516. De φ-coëfficiënt: r χ2 φ= n Opdat we de χ2 -waarde niet zouden moeten overtypen (wat afrondingsfouten kan opleveren) kunnen we deze uit de output extraheren. We wijzen eerst het volledige resultaat toe aan een nieuwe variabele, en vervolgens vinden we de χ2 -waarde door aan de naam van de nieuwe variabele $statistic toe te voegen, en dit binnen de eerder vermelde functie unname(variabele): > chitest <- chisq.test(tabel) > chi2 <- unname(chitest$statistic) > chi2 [1] 55.25158 > sqrt(chi2 / sum(tabel)) [1] 0.3928537 De contingentiecoëfficiënt C: s χ2 C= 2 χ +n > sqrt(chi2 / (chi2 + sum(tabel))) [1] 0.3656496

33

Rangcorrelatieco¨ effici¨ enten Twee experten in het domein van de Nihilistische Statistiek, namelijk I. Analasiac en Sir J. Poetiers bevinden zich op een congres. Tijdens de receptie discussiëren beide heren over de kwaliteit van de hapjes die worden rondgebracht. (De ratings van laag naar hoog zijn: 1: “Afschuwelijk”, 2: “Niet lekker”, 3: “Normaal”, 4: “Lekker”, 5: “Heel lekker”, 6: “Heerlijk”). Hapje Toostje met belugakaviaar Foie gras op zuurdesembrood Kikkerbilletjes Gefrituurde kreeftenpootjes Oesters ` a la meunière Slaatje met crème van witte truffel

Mening I. Analasiac Afschuwelijk Lekker Heerlijk Heel lekker Niet lekker Normaal

Mening Sir J. Poetiers Niet lekker Heerlijk Lekker Normaal Afschuwelijk Heel lekker

R kan uiteraard de volgorde van de waarderingen niet kennen, dus geven we meteen de bijhorende rangscores in: > analasiac <- c(1, 4, 6, 5, 2, 3) > poetiers <- c(2, 6, 4, 3, 1, 5) Correlatiecoëfficiënten kunnen we opvragen met deze functie:

cor ( x, y, method = "pearson")


Bereken een (rang)correlatieco¨ effici¨ ent.

Vector. Vector. Type correlatiecoëfficiënt: man”.

”pearson”, ”kendall” of ”spear-

De standaardwaarde voor de parameter method (Pearson) is hier uiteraard niet van toepassing omdat het ordinale schalen zijn. Wat wel kan: Kendalls Tau (τ ): > cor(analasiac, poetiers, method = "kendall") [1] 0.3333333 Spearmans rs : > cor(analasiac, poetiers, method = "spearman") [1] 0.4857143 Beide wijzen op een matig positief verband tussen de beoordelingen van onze twee experten. De twee zijn het dus in zekere mate met elkaar eens. Correlatie en regressie Het bestand “lengtegewicht.txt” bevat lengte (in centimeter) en gewicht (in kilogram) van een groep personen. > data <- read.table("lengtegewicht.txt", header = TRUE) Eerder hadden we summary() al toegepast op een vector, je kan dit ook gebruiken voor een volledig databestand: > summary(data) lengte Min. :150.9 1st Qu.:161.3

gewicht Min. :50.35 1st Qu.:60.29 34

Median :165.0 Mean :167.2 3rd Qu.:172.8 Max. :190.3

Median :64.13 Mean :65.06 3rd Qu.:67.47 Max. :92.10

Wat ons meteen de respectievelijke namen van de variabelen toont in het bestand, zodat we deze kunnen gebruiken: > lengte <- data$lengte > gewicht <- data$gewicht We zetten de variabele “lengte” uit op de X-as van het spreidingsdiagram, en “gewicht” op de Y -as (plot() werd in het vorige hoofdstuk reeds beschreven): > plot(lengte, gewicht) Zoals kon worden verwacht is er een duidelijk positief verband te zien. We hebben twee formules gezien om de correlatiecëfficiënt r te berekenen: door middel van het product van de gestandaardiseerde waarden, en via de covariantie (en uiteraard zullen we zien hoe R het in één keer kan). Voor beiden hebben we de standaardafwijkingen nodig: > s_lengte <- sqrt(mean((lengte - mean(lengte))^2)) > s_gewicht <- sqrt(mean((gewicht - mean(gewicht))^2)) Berekening van r door middel van de gestandaardiseerde waarden: Pn (Zx · Zy ) r = i=1 n Met: Zx =

xi − x ¯ sx

> Zlengte <- (lengte - mean(lengte)) / s_lengte > Zgewicht <- (gewicht - mean(gewicht)) / s_gewicht De essentie van standaardiseren is dat de gestandaardiseerde waarde Zx van x een gemiddelde heeft dat gelijk is aan nul en een standaardafwijking gelijk aan één, wat betekent dat Zx weergeeft hoeveel standaardafwijkingen x groter is dan x ¯ . Ga dit na voor beide standaard(meet)waarden. Vervolgens moeten we het product nemen van al de gestandaardiseerde waarden, en deze optellen en delen door het aantal, of nog: het gemiddelde van deze producten nemen: > mean(Zlengte * Zgewicht) [1] 0.7367061 Berekening van r door middel van de covariantie: r=

sxy sxy =q sx · sy s2x · s2y

Met covariantie: n

sxy =

1X (xi − x ¯) (yi − y¯) n i=1

> cov_lengte_gewicht <- mean((lengte - mean(lengte)) * (gewicht - mean(gewicht))) > cov_lengte_gewicht [1] 59.18579 35

> correl <- cov_lengte_gewicht / (s_lengte * s_gewicht) > correl [1] 0.7367061

Populatiecovariantie berekenen (schatter!). cov ( x, y)


We kunnen de covariantie ook opvragen door middel van:

Vector. Vector.

Hou er rekening mee dat R net als bij de berekening van de variantie een schatting maakt van de parameter in de populatie (deelt door n − 1), en dus niet de beschrijving van de steekproef geeft (delen door n). Als je echter de correlatie berekent op basis van varianties en covariantie berekent als schatters van de karakteristieke waarden voor de populatie zal je toch hetzelfde resultaat bekomen als wanneer je diezelfde correlatie berekent op basis van beschrijvende waarden (ga na waarom): > correl2 <- cov(lengte, gewicht) / sqrt(var(lengte) * var(gewicht)) > correl2 [1] 0.7367061 Uiteraard kan R ook in één keer de correlatie (= correlatiecoëfficiënt van Pearson) berekenen, dat gaat zoals eerder vermeld. Aangezien “pearson” de standaardwaarde is kunnen we dat hier achterwege laten. > cor(lengte, gewicht) [1] 0.7367061 De regressierechte waarmee we een variabele Y trachten te voorspellen aan de hand van X R.R.y/x heeft de vorm Y = b0 + b1 X, met: b0 = y¯ − (b1 · x ¯)

b1 =

sxy s2x

> b1 <- cov_lengte_gewicht / s_lengte^2 > b1 [1] 0.7473029 > b0 <- mean(gewicht) - b1 * mean(lengte) > b0 [1] -59.89801

Lineaire regressie.


Ofwel, in één keer laten berekenen door R met de functie:

lm ( formula)

Y -variabele ˜ X-variabele.

15

> lm(gewicht ~ lengte) Call: lm(formula = gewicht ~ lengte) Coefficients: 15 Opgelet, tussen de te verklaren variabele en de verklarende variabele staat een tilde en geen komma. Let ook op de volgorde van de twee variabelen. In tegenstelling tot de correlatie maakt dit voor de regressierecht wel degelijk verschil, er zijn immers 2 verschillende regressierechten.

36

(Intercept) -59.8980

lengte 0.7473

Om het resultaat grafisch voor te stellen vragen we opnieuw het spreidingsdiagram op, maar nu kunnen we de regressierechte toevoegen door middel van:

abline ( a, b, lwd = 1, col = "black")


Teken een rechte op een grafische voorstelling. Constante. Richtingscoëfficiënt. Lijndikte. Kleur.

> plot(lengte, gewicht) > abline(b0, b1) Tenslotte kunnen we ook voorspellingen doen. Bijvoorbeeld, stel dat iemand 170 cm groot is, hoeveel kilogram verwachten we dan dat die zou wegen? > b0 + b1 * 170 [1] 67.14347 Merk op dat alle voorspellingen punten op de regressierechte zijn. We kunnen de bekomen voorspelling grafisch voorstellen met een punt:

points ( x, y, lwd = 1, col = "black")


Teken ´ e´ en of meer punten op een grafische voorstelling. Getal of vector van x-coördinaten. Getal of vector van y-coördinaten. Lijndikte. Kleur.

70 50

60

gewicht

80

90

> points(170, b0 + b1 * 170, lwd = 10, col = "red")

150

160

170

180

190

lengte

Het kwadraat van de correlatiecoëfficiënt geeft aan welke proportie van de variantie in de Y -variabele wordt verklaard door de X-variabele: > r^2 [1] 0.5427358 We kunnen dit als volgt interpreteren: 54% van de variantie in gewicht wordt verklaard door de lengte van een persoon. De overige 46% van de variantie wordt dus door andere factoren bepaald (eetgewoonten, ...).

37

Oefeningen 3.2.1 Een assistent merkt op dat enkele studenten grote hoeveelheden chocoladetaart zitten te eten tijdens de oefeningenles. De studenten beweren dat taart eten hun vermogen om oefeningen te maken doet toenemen. Om dit te testen stelt de assistent een experiment op: aan enkele studenten geeft hij een (verschillende) hoeveelheid chocoladetaart en hij noteert hoeveel oefeningen iedere student maakt gedurende één les. Student A B C D E F P

Gram chocoladetaart (xi ) 52 100 80 96 120 128 576

Aantal oefeningen (yi ) 10 7 13 13 21 14 78

a) Bereken en interpreteer de correlatiecoëfficiënt. b) Laat R een spreidingsdiagram maken. Bereken bijhorende regressierechte en voeg deze toe aan het spreidingsdiagram. c) Stel dat een student 150 gram chocoladetaart zou eten, hoeveel oefeningen verwachten we dan dat die zou maken? d) Geef en interpreteer r2 . 3.2.2 De residuën van de regressie R.R.y/x zijn het verschil tussen iedere waarde yi en de bijhorende voorspelde waarde yî . De residuën geven weer welk deel van de variantie in Y we niet kunnen verklaren door X, en kunnen gebruikt worden om r2 te berekenen: Pn Pn 2 2 (yi − y¯) − i=1 (yi − yˆ) r2 = i=1 Pn 2 ¯) i=1 (yi − y Gebruik de gegevens uit voorgaande oefening. Bereken r2 aan de hand van deze formule. Om de voorspellingen yî te vinden neem je de waarden xi en doe een voorspelling aan de hand van de richtingscoëfficiënt en de constante (het intercept). 3.2.3 Twee psychiaters bespreken de ernst van het depressieve beeld bij zes patiënten. De mogelijke meetwaarden zijn: 1: “Zwaar depressief”; 2: “Matig depressief”; 3: “Licht depressief”; 4: “Dystymie”; 5: “Gelukkig”; 6: “Zeer gelukkig”. Bereken en interpreteer de correlatiecoëfficiënten die hier van toepassing zijn. Patiënt Psychiater 1 Psychiater 2

1 Licht depr. Matig depr.

2 Zwaar depr. Zwaar depr.

3 Dystymie Licht depr.

4 Zeer gelukkig Zeer gelukkig

5 Gelukkig Gelukkig

6 Matig depr. Dystymie

3.2.4 Onderstaande tabel geeft de aantallen deelnemers aan een congres ingedeeld naar geslacht en of zij al-dan-niet roken. Man Vrouw

Roker 45 42

Niet-roker 343 357

a) Bereken λ om na te gaan of we geslacht beter kunnen voorspellen indien we het rookgedrag kennen. 38

b) Bereken de Q-coëfficiënt. c) Bereken χ2 . d) Bereken de φ-coëfficiënt. e) Bereken de contingentiecoëfficiënt C. 3.2.5 Groep X Y

A 1 1265

A 2 1234

A 3 1813

A 4 1125

A 5 930

A 6 1050

A 7 893

A 8 753

A 9 822

A 10 836

Groep X Y

B 11 1638

B 12 1915

B 13 1539

B 14 1318

B 15 1246

B 16 1238

B 17 1393

B 18 1136

B 19 1298

B 20 1045

a) Laat R de regressierechte tekenen waarmee je Y voorspelt aan de hand van X. b) Doe het zelfde, maar enkel voor groep A, en vervolgens enkel voor groep B. c) Hoe noem je dit verschijnsel? Tip: gebruik verschillende kleuren voor de verschillende regressierechten. 3.2.6 Stel twee vectoren op met elk zes elementen opdat de correlatie tussen beiden tussen -0.60 en -0.80 ligt. 3.2.7 Onderstaande tabel toont hoe vaak een voetganger ongedeerd, gewond of dood is na aanrijding door een kleine gezinswagen, een grote gezinswagen, of een terreinwagen. Kleine gezinswegen Grote gezinswagen Terreinwagen

Ongedeerd 138 38 27

Gewond 92 39 53

Dood 17 12 82

a) Bereken λ om na te gaan of we de effecten van de aanrijding (ongedeerd, gewond, of dood) beter kunnen voorspellen als we het type wagen kennen. b) Bereken de Q-coëfficiënt (dit kan enkel 2 × 2, neem daarom enkel gezinswagens in beschouwing en neem de slachtoffers die gewond of dood zijn samen). c) Bereken χ2 . d) Bereken de φ-coëfficiënt. e) Bereken de contingentiecoëfficiënt C.

39

4

Kans- en steekproevenverdelingen

4.1

Kansverdelingen

De verschillende functies die we nodig hebben om in R met kansverdelingen te werken hebben steeds dezelfde vorm, namelijk één letter gevolgd door de benaming van de functie. We bespreken ze niet allemaal voor iedere verdeling apart. dverdeling() pverdeling() qverdeling() rverdeling()

Kansdichtheidsfunctie Cumulatieve kansdichtheidsfunctie Proportie Random getal

P [X = x] P [X ≤ x]

Kansdichtheid binomiaalverdeling.

Waarde. Aantal. Kans op succes.

Cumulatieve kansdichtheid binomiaalverdeling. Waarde. Aantal. Kans op succes.

Waarde waaronder een gegeven proportie ligt bij een binomiaalverdeling. Proportie. Aantal. Kans op succes.

Genereer random getallen uit een binomiaalverdeling. rbinom ( n, size, prob)


qbinom ( p, size, prob)


pbinom ( q, size, prob)


dbinom ( x, size, prob)


Binomiaalverdeling

Aantal random waarden. Aantal. Kans op succes.

Met 20 trekkingen en een kans op succes gelijk aan 0.7 per trekking is het aantal successen X binomiaal verdeeld X : B(20, 0.7). De kans om een bepaald aantal successen te hebben vinden we met dbinom(). Bijvoorbeeld de kans op 10 successen bij 20 trekkingen met kans op succes van 0.7 bedraagt: > dbinom(10, 20, 0.7) [1] 0.03081708 We kunnen natuurlijk in één keer de kansen op alle mogelijke uitkomsten opvragen: > kansen <- dbinom(seq(0, 20), 20, 0.7) > kansen [1] 3.486784e-11 1.627166e-09 3.606885e-08 5.049639e-07 5.007558e-06 3.738977e-05 [7] 2.181070e-04 1.017833e-03 3.859282e-03 1.200665e-02 3.081708e-02 6.536957e-02

40

[13] 1.143967e-01 1.642620e-01 1.916390e-01 1.788631e-01 1.304210e-01 7.160367e-02 [19] 2.784587e-02 6.839337e-03 7.979227e-04 Of, iets leesbaarder: > round(kansen, 6) [1] 0.000000 0.000000 0.000000 0.000001 0.000005 0.000037 0.000218 0.001018 [9] 0.003859 0.012007 0.030817 0.065370 0.114397 0.164262 0.191639 0.178863 [17] 0.130421 0.071604 0.027846 0.006839 0.000798

Maak een barplot.


Overzichtelijker in een histogram. Indien je barplot(kansen, space = 0) oproept zal je merken dat er geen waarden staan op de X-as. Dat is logisch, kansen is hier geen frequentietabel (wat eigenlijk vereist is voor barplot(), maar slechts een vector met kansen waarbij de eerste de kans op waarde nul voorstelt en zo verder. We voegen daarom een nieuwe optionele parameter toe:

barplot ( height, space = NULL,

Vector met waarden die de hoogte van de balkjes bevat. Spatie tussen de balkjes (standaard wordt een spatie gebruik, wat correct is op nominale en ordinale schaal). Een vector met benamingen bijhorende de balkjes. Toevoegen aan bestaande grafische voorstelling.

names.arg = NULL, add = FALSE)

0.00

0.05

0.10

0.15

> barplot(kansen, space = 0, names.arg = seq(0, 20))

0

2

4

6

8

10

12

14

16

18

20

En uiteraard: > sum(kansen) [1] 1 De kans om bij deze verdeling tussen 10 en 14 successen te hebben (grenzen inbegrepen): P [10 ≤ X ≤ 14] = P [X = 10] + P [X = 11] + P [X = 12] + P [X = 13] + P [X = 14] > sum(dbinom(seq(10, 14), 20, 0.7)) [1] 0.5664844 Of, via de cumulatieve kansdichtheidsfunctie (pbinom(waarde, aantal, kans op succes)): P [10 ≤ X ≤ 14] = P [X ≤ 14] − P [X ≤ 9] > pbinom(14, 20, 0.7) - pbinom(9, 20, 0.7) [1] 0.5664844

41

Stel dat we de mediaan van deze verdeling willen kennen (voor de mediaan geldt dat de proportie waarnemingen die kleiner zijn = 0.5) bijgevolg wordt de functie : qbinom(proportie, aantal, kans op succes)): > qbinom(0.5, 20, 0.7) [1] 14

Kansdichtheid normaalverdeling.

Waarde. Gemiddelde. Standaardafwijking.

Cumulatieve kansdichtheid normaalverdeling. Waarde. Gemiddelde. Standaardafwijking.

Waarde waaronder een gegeven proportie ligt bij een normaalverdeling. Proportie. Gemiddelde. Standaardafwijking.

Genereer random getallen uit een normaalverdeling. rnorm ( n, mean = 0, sd = 1)


qnorm ( p, mean = 0, sd = 1)


pnorm ( q, mean = 0, sd = 1)


dnorm ( x, mean = 0, sd = 1)


Normaalverdeling

Aantal random waarden. Gemiddelde. Standaardafwijking.

We kunnen bijvoorbeeld de (kans)dichtheid voor een bepaalde waarde opvragen, of nog, hoe hoog de curve van de verdeling op dat punt boven de X-as ligt. Merk op dat de standaardwaarden van de optionele parameters deze van de standaardnormaalverdeling zijn. Bijvoorbeeld, de hoogte van de curve van de standaardnormaalverdeling voor x = 0: > dnorm(0) [1] 0.3989423 We kunnen dus ook een kansdichtheidsfuntie tekenen (de x in dnorm() is hier een verwijzing naar de X-as en dus niet naar een object met de naam x). Bijvoorbeeld voor X : N (10, 3):16 > normfun <- function(x) dnorm(x, mean = 10, sd = 3) > plot(normfun, from = 0, to = 20) Het is aan te raden om van minstens drie standaardafwijkingen onder, tot drie standaardafwijkingen boven het gemiddelde te tekenen om een mooi resultaat te bekomen. Stel, we willen de oppervlakte tussen x = 5 en x = 8 inkleuren, dan kan dit door middel van polygon()17 . 16 Zie 17 Zie

“Basistechnieken: Analyse” voor uitleg over het tekenen van een functie. “Basistechnieken: Analyse” voor bespreking van de vectoren met co¨ ordinaten.

42

r <- seq(from = 5, to = 8, by = 0.01) x <- c(5, r, 8) y <- c(0, dnorm(r, mean = 10, sd = 3), 0) polygon(x, y, col = "orange")

0.08 0.00

0.04

normfun

0.12

> > > >

0

5

10

15

20

x

Deze oppervlakte komt overeen met de kans dat deze normaalverdeling een waarde aanneemt tussen x = 5 en x = 8. De bepaalde integraal van de normaalverdeling met ondergrens 5 en bovengrens 8 geeft ons die oppervlakte (zie “Basistechnieken”). P [5 ≤ X ≤ 8] > f <- function(x) dnorm(x, mean = 10, sd = 3) > integrate(f, lower = 5, upper = 8) 0.2047022 with absolute error < 2.3e-15 Of, via de cumulatieve kansdichtheidsfunctie: P [X ≤ 8] − P [X ≤ 5] > pnorm(8, mean = 10, sd = 3) - pnorm(5, mean = 10, sd = 3) [1] 0.2047022 Aangezien de normaalverdeling van −∞ tot ∞ loopt moet je dat ook kunnen ingeven. Gebruik hiervoor respectievelijk -Inf en Inf, bijvoorbeeld P [12 > X > ∞]: > integrate(f, lower = 12, upper = Inf) 0.2524925 with absolute error < 2.2e-05 Wat natuurlijk ook als volgt kan: P [12 > X > ∞] = 1 − P [X ≤ 12] > 1 - pnorm(12, mean = 10, sd = 3) [1] 0.2524925 Stel, we willen het derde kwartiel van N (10, 3) opvragen (dat is ook wel het 75e percentiel): > qnorm(.75, mean = 10, sd = 3) [1] 12.02347 Je kan random waarden genereren die worden getrokken uit een bepaalde normaalverdeling. Bijvoorbeeld, geef 10 random getallen uit de standaardnormaalverdeling (uw resultaat zal uiteraard anders zijn): > rnorm(10) [1] 0.80192789 -0.89995005 [7] -0.50133079 0.05944035

1.26771738 0.63118137 -2.17769884 -1.35954189 2.57231769 -1.37150059

43

Waarde. Aantal vrijheidsgraden.

Cumulatieve kansdichtheid t-verdeling. pt ( q, df)


Waarde waaronder een gegeven proportie ligt bij een t-verdeling. qt ( p, df)

Proportie. Aantal vrijheidsgraden. JorisPieters,PeterTheuns,AlainIsaac

Genereer random getallen uit een t-verdeling. rt ( n, df)


dt ( x, df)


Kansdichtheid t-verdeling.


t-verdeling

Aantal random waarden. Aantal vrijheidsgraden.

Bijvoorbeeld, de kans om bij een t-verdeling met 15 vrijheidsgraden een waarde van maximum 0.3 te bekomen: > pt(0.3, 15) [1] 0.6158527 De kans om bij een t-verdeling met 20 vrijheidsgraden een waarde groter dan 2.2 te bekomen: > 1 - pt(2.2, 20) [1] 0.01986429 De kans om bij een t-verdeling met 100 vrijheidsgraden een waarde tussen -1 en 1 te bekomen: > pt(1, 100) - pt(-1, 100) [1] 0.6802758

Kansdichtheid χ2 -verdeling. dchisq ( x, df)


χ2 -verdeling


pchisq ( q, df)


Cumulatieve kansdichtheid χ2 -verdeling.


44

qchisq ( p, df)

Proportie. Aantal vrijheidsgraden.


Genereer random getallen uit een χ2 -verdeling. rchisq ( n, df)


Waarde waaronder een gegeven proportie ligt bij een χ2 -verdeling.


Bijvoorbeeld, het 95e percentiel van een χ2 -verdeling met 10 vrijheidsgraden: > qchisq(.95, 10) [1] 18.30704

df ( x, df1, df2)

Waarde. Aantal 1e vrijheidsgraden. Aantal 2e vrijheidsgraden.

pf ( q, df1, df2)


Cumulatieve kansdichtheid F -verdeling.


Waarde waaronder een gegeven proportie ligt bij een F -verdeling. Proportie. Aantal 1e vrijheidsgraden. Aantal 2e vrijheidsgraden.

Genereer random getallen uit een F -verdeling. rf ( n, df1, df2)


qf ( p, df1, df2)


Kansdichtheid F -verdeling.


F -verdeling

Aantal random waarden. Aantal 1e vrijheidsgraden. Aantal 2e vrijheidsgraden.

Bijvoorbeeld een verhouding van 2 of lager bij 8 vrijheidsgraden in de teller en 10 in de noemer heeft onderstaande kans: > pf(2, 8, 10) [1] 0.8492264 Convergenties Een binomiaalverdeling met een redelijk groot aantal pogingen (n), en een kans op succes (p) die niet te dicht bij 0 of 1 ligt kan goed benaderd worden door een normaalverdeling met verwachting np en 45

p standaardafwijking np(1 − p). Stel, een multiple choice examen telt 40 vragen met elk vier mogelijke antwoorden waarvan er steeds slechts één juist is, zonder giscorrectie (de kans op succes bedraagt dus 1/4). We stellen dit grafisch voor, en bepalen de kans dat iemand minstens 10 juiste antwoorden geeft indien die persoon elke vraag zou gokken, en dit zowel exact als via de normaalbenadering. Merk op dat we bij het tekenen van de barplot een extra parameter (add) introduceren die maakt dat de balkjes op de curve getekend worden. De reden dat nog 0.5 wordt toegevoegd is omdat R het begin van de balkjes op de eenheden zet (en niet het midden)(zie cursus statistiek i.v.m. de “continu¨ıteitscorrectie).

0.10 0.05 0.00

normfun

0.15

> 1 - pbinom(9, 40, 0.25) [1] 0.5604603 > kansen <- dbinom(seq(0, 40), 40, 0.25) > normfun <- function(x) dnorm(x, mean = (40 * 0.25) + 0.5, sd = sqrt(40 * 0.25 * 0.75)) > plot(normfun, from = 0, to = 40, axes = FALSE) > barplot(kansen, space = 0, names.arg = seq(0, 40), add = TRUE)

0

2

4

6

8 10

13

16

19

22

25

28

31

34

37

40

x

Het is duidelijk dat de normaalverdeling goed past. Indien we nu de kans op slagen berekenen aan de hand van de normaalverdeling (de 0.5 is continu¨ıteitscorrectie), dan krijgen we nagenoeg hetzelfde resultaat als voorheen: > 1 - pnorm(9.5, mean = 40 * 0.25, sd = sqrt(40 * 0.25 * 0.75)) [1] 0.5724339 Oefeningen 4.1.1 Laten we er vanuit gaan dat er evenveel jongens als meisjes geboren worden. Op een dag worden aan het gemeentehuis 24 geboorten gemeld. Beantwoord onderstaande vragen: a) Wat is de kans dat er 10 jongens zijn? b) Wat is de kans dat er meer meisjes dan jongens zijn? c) Wat is het mediaan aantal jongens? 4.1.2 De scores op een bepaalde test zijn normaal verdeeld met gemiddelde = 500 en standaardafwijking = 50. Beantwoord onderstaande vragen met behulp van R: a) De kans om meer dan 500 te scoren. b) De kans om tussen 400 en 500 te scoren. c) De kans om tussen 400 en 600 te scoren. d) De kans om exact 500 te scoren.

46

e) Geef en interpreteer het 95e percentiel. 4.1.3 Een bepaald type ledlampen heeft een normaal verdeelde levensduur met een verwachting van 50000 branduren en een standaardafwijking van 2850 branduren. Als een lamp meer dan 55000 uren volbrengt wordt deze beschouwd als zeer goed. Ga na hoeveel procent van de lampen aan dit criterium voldoen en geef grafische weer. 4.1.4 Op de leeftijd van 10 jaar is de lichaamslengte van meisjes normaal verdeeld met gemiddelde = 143.4 cm en standaardafwijking = 6.9 cm. a) Indien de lichaamslengte van een meisje zich verder dan 2.5 standaardafwijkingen van het gemiddelde bevindt vereist dit bijkomende medische aandacht. Bepaal onder hoeveel centimeter en boven hoeveel centimeter dit het geval is. b) Wat is de kans om buiten die grenzen te vallen? c) Maak een bijhorende grafische voorstelling.

4.2

Steekproevenverdelingen

Normaal verdeelde populatie Stel, een bepaalde variabele is normaal verdeeld met verwachting = 1200 en standaardafwijking = 80. Eerder hebben we gezien hoe een reeks van 16 random waarden te bekomen uit een bepaalde verdeling (omdat het random waarden zijn zullen jouw resultaten uiteraard iets anders zijn dan wat hier wordt getoond): > x <- rnorm(16, mean = 1200, sd = 80) En daar kunnen we dan uiteraard het gemiddelde van nemen: > mean(x) [1] 1193.926 Iedere keer als je het gemiddelde van een reeks getallen x opvraagt zal je een iets andere uitkomst bekomen. Dit gemiddelde is dus zelf een toevalsvariabele die we de steekproevenverdeling X zullen noemen. Hoe is X verdeeld? We kunnen at random 1000 “steekproeven” genereren, en hiervan gemiddelde en standaardafwijking berekenen. De manier waarop we X hier aanmaken moet je niet kennen.18 > X <- 0 > for(i in 1:1000){X[i] <- mean(rnorm(16, mean = 1200, sd = 80))} > mean(X) [1] 1198.886 19

> sd(X) [1] 20.69369 Dit is (bijna) exact wat we verwachten. De som van onafhankelijke normaalverdelingen N : (µ, σ) is σ √ immers ook nomaal verdeeld, maar met parameters N : µ, n . Het gemiddelde van de steekproevenverdeling verschilt hier dan ook nauwelijks van 1200, en de standaardafwijking hiervan zit zeer zicht bij de veronderstelde √8016 = 20. 18 We kunnen hier de functie rep() niet gebruiken omdat we dan 1000 keer hetzelfde getal zouden verkrijgen. Daarom gebruiken we hier een zogenaamde loop. 19 Omdat we nu bezig zijn met statistische analyse en niet meer met beschrijvende statistiek zullen we voor de variantie en standaardafwijking steeds de schatter van de populatieparameter gebruiken. Voor de variantie betekent dit dus de kwadratensom gedeeld door n − 1 in plaats van gedeeld door n. Voor de variantie wordt de formule dan s2 = 2 1 Pn (x − x ¯ ) , zodat we de functies var() en sd() kunnen gebruiken, die zoals eerder vermeld de schatter van i i=1 n−1 respectievelijk populatie variantie en standaardafwijking geven.

47

Niet-normaal verdeelde populatie Stel, we hebben onderstaande functie: f (x) =

e−x/5 5

Een dergelijke exponentiële functie is kenmerkend voor bijvoorbeeld de tijd die het duurt voor een bepaalde gebeurtenis zich voltrekt. Bovenstaande verdeling zou de tijd kunnen voorstellen tot een volgende klant aan een winkel verschijnt, met een gemiddelde en standaardafwijking van elk vijf minuten. Het spreekt voor zich dat een normaalverdeling uiteraard niet mogelijk zou zijn: ten eerste is er een probleem qua begrenzing aangezien negatieve tijd uiteraard niet mogelijk is, ten tweede zijn dergelijke wachttijden sterk asymmetrisch zoals duidelijk te zien is als je een grafiek maakt:

0.00

0.05

0.10

0.15

0.20

> f <- function(x) exp(-x / 5) / 5 > plot(f, from = 0, to = 50, xlab = "", ylab = "")

0

10

20

30

40

50

Stel dat we nu verschillende steekproeven zouden nemen uit een dergelijke niet-normaal verdeelde variabele, wat zou dan de verdeling worden van de bekomen gemiddelden? Volgens de centrale limietstelling zouden we ook hier een normaal verdeelde steekproevenverdeling moeten bekomen, zij het wel pas wanneer het voldoende grote steekproeven betreft. We kiezen hier voor steekproeven van grootte 100. Je kan zelf nagaan wat het effect zou zijn mocht je met kleinere (of grotere) steekproefgrootten werken. > X <- 0 > for(i in 1:1000){X[i] <- mean(rexp(100, rate = 1/5))} > mean(X) [1] 4.995602 > sqrt(mean((X - mean(X))^2)) [1] 0.5110099 Ook voor een niet normaal verdeelde populatie lijkt dus te gelden dat de steekproevenverdeling hetzelfde gemiddelde heeft en een standaardafwijking gelijk aan de populatiestandaardafwijking gedeeld door de √ vierkantswortel van de steekproefgrootte: 5/ 100 = 0.5 ≈ 0.5110099. Maar is deze steekproevenverdeling ook normaal verdeeld? Een histogram kan verduidelijking brengen. We introduceren twee nieuwe optionele parameter bij het histogram. De eerste laat ons kiezen tussen absolute en relatieve frequenties op de Y -as, wat nodig is als we in dezelfde figuur ook een normaalverdeling willen tekenen, de tweede laat ons kiezen hoeveel klassen we willen zien aangezien R er standaard nogal weinig neemt: > plot(function(x) dnorm(x, mean = 4.995602, sd = 0.5110099), from = 3, to = 7) > hist(X, freq = FALSE, nclas = 20, add = TRUE)

48

0.8 0.6 0.4 0.2 0.0

function(x) dnorm(x, mean = 4.995602, sd = 0.5110099)

3

4

5

6

7

x

Het is duidelijk dat de steekproevenverdeling, van gemiddelden van steekproeven die uit een nietnormaal verdeelde populatie worden getrokken, dus eveneens bij benadering normaal verdeeld is. Je kan zelf nagaan dat dit echter niet geldt als je kleine steekproeven trekt.

49

5

Statistische analyse van het gemiddelde

5.1

Betrouwbaarheidsinterval

Betrouwbaarheidsinterval voor het gemiddelde. t.test ( x, conf.level = 0.95)


Het betrouwbaarheidsinterval voor het gemiddelde bekomen we door middel van:

Vector. Betrouwbaarheidsniveau.

Bijvoorbeeld, de file “erytrocyten.txt” bevat metingen van de hoeveelheid erytrocyten (rode bloedlichaampjes) per picoliter bloed bij 36 gezonde volwassenen. We bepalen op basis van deze data het 95% en het 90% betrouwbaarheidsinterval voor het gemiddeld aantal erytrocyten per picoliter bij gezonde volwassenen. Het resultaat vinden we telkens vlak onder ...percent confidence interval in de R output voor de t-test. Voor het 95% betrouwbaarheidsinterval kunnen we gebruik maken van de standaardwaarde voor conf.level, als je niets verder specificeert wordt namelijk het 95% betrouwbaarheidsinterval bepaald: > erytrocyten <- read.table("erytrocyten.txt", header = TRUE) > t.test(erytrocyten) One Sample t-test data: erytrocyten t = 110.1983, df = 35, p-value < 2.2e-16 alternative hypothesis: true mean is not equal to 0 95 percent confidence interval: 4.744837 4.922940 sample estimates: mean of x 4.833889 En bijgevolg is het 95% betrouwbaarheidsinterval : [4.744837 ; 4.922940]. Voor het 90% BI moeten we iets meer parameters ingeven: > t.test(erytrocyten, conf.level = 0.90) One Sample t-test data: erytrocyten t = 110.1983, df = 35, p-value < 2.2e-16 alternative hypothesis: true mean is not equal to 0 90 percent confidence interval: 4.759775 4.908003 sample estimates: mean of x 4.833889 Het 90% betrouwbaarheidsinterval voor het gemiddeld aantal erytrocyten per picoliter bij gezonde volwassenen is dus [4.759775 ; 4.908003].

50

Oefeningen 5.1.1 Bepaal de 95% betrouwbaarheidsintervallen voor lengte en gewicht op basis van de data over een steekproef deelnemers in het bestand “lengtegewicht.txt”. Welk van de twee intervallen is het breedste? Hoe komt dit?

5.2

Toetsen van een hypothese omtrent µ0

Toets een hypothese omtrent µ0 . t.test ( x, alternative = "two.sided", mu = 0)


Dezelfde functie, maar nu met andere parameters, wordt gebruikt om een hypothese (H0 : µ = µ0 ) te toetsen omtrent een vooropgestelde waarde µ0 :

Vector. Alternatieve hypothese, andere mogelijke waarden zijn “less” en “greater”. Waarde die we wensen te toetsen.

Bij de parameter voor de alternatieve hypothese gebruiken we "less" of "greater" als de alternatieve hypothese respectievelijk < of > bevat. We kunnen bijvoorbeeld toetsen of het gemiddeld aantal erytrocyten per picoliter bloed verschilt van 5 (H0 : µ = 5): > t.test(erytrocyten, mu = 5) One Sample t-test data: erytrocyten t = -3.7868, df = 35, p-value = 0.0005757 alternative hypothesis: true mean is not equal to 5 95 percent confidence interval: 4.744837 4.922940 sample estimates: mean of x 4.833889 Opgelet: de output ziet er exact hetzelfde uit als voorheen, maar nu kijken we niet naar het betrouwbaarheidsinterval, maar naar de lijnen die daarboven staan. Daar vinden de t-waarde, het aantal vrijdheidsgraden (df), de overschrijdingskans (p-value), en de alternatieve hypothese. We vinden t(35) = −3.79, p < 0.001. Gezien de overschrijdingskans kleiner is dan 0.0520 verwerpen we de nulhypothese. We rapporteren: “Het gemiddelde aantal erytrocyten per picoliter bloed verschilt significant van 5 (gemiddelde = 4.83, SD = 0.263, t(35) = −3.79, p < .00121 )”. De standaardafwijking wordt in de rapportering vermeld maar die is hier niet gegeven. Je weet ondertussen wel zelf hoe je die kan opvragen. Een wetenschappelijk artikel stelt dat het gemiddeld aantal erytrocyten per picoliter bij gezonde volwassenen kleiner is dan 5 (Ha : µ < 5). Toets nu of je deze stelling kunt bijtreden met de gegevens in de onderzochte steekproef: > t.test(erytrocyten, mu = 5, alternative = "less") One Sample t-test 20 Voor de hele verdere tekst geldt dat wanneer het betrouwbaarheidsniveau niet expliciet vermeld wordt 0.95 genomen wordt. 21 Bij waarden die nooit groter kunnen zijn, zoals p-waarden wordt de nul voor het decimaalteken weggelaten.

51

data: erytrocyten t = -3.7868, df = 35, p-value = 0.0002878 alternative hypothesis: true mean is less than 5 95 percent confidence interval: -Inf 4.908003 sample estimates: mean of x 4.833889 We aanvaarden de alternatieve hypothese, het gemiddeld aantal erytrocyten per picoliter bloed is kleiner dan 5 (gemiddelde = 4.83, SD = 0.263, t(35) = −3.79, p < .001). Stel dat je ook de steekproefgrootte (n) wenst te rapporteren, hoe vraag je die dan op? Probeer dit zelf! De andere éénzijdige toets: ga na of µ beduidend groter is dan 5 (Ha : µ > 5): > t.test(erytrocyten, mu = 5, alternative = "greater") One Sample t-test data: erytrocyten t = -3.7868, df = 35, p-value = 0.9997 alternative hypothesis: true mean is greater than 5 95 percent confidence interval: 4.759775 Inf sample estimates: mean of x 4.833889 Uiteraard kan deze toets niet significant zijn aangezien we reeds hadden aangetoond dat het gemiddelde kleiner is dan 5. Als het steekproefgemiddelde (4.8339) kleiner is dan de vooropgestelde waarde 5 kan je natuurlijk nooit statistische evidentie hebben dat het gemiddelde in de populatie groter zou zijn: “Het gemiddeld aantal erytrocyten per picoliter bloed is niet groter dan 5 (gemiddelde = 4.83, SD = 0.263, t(35) = −3.79, p = .999)” Merk op dat tussen de drie toetsen enkel de overschrijdingskans verschilt. Aangezien je deze toetsen eerder ook met de hand hebt uitgevoerd begrijp je natuurlijk waarom. In de praktijk ga je nooit deze drie toetsen allemaal uitvoeren, maar steeds enkel diegene die overeenkomt met je onderzoekshypothese. Oefeningen 5.2.1 Open “lengtegewicht.txt” en toets: a) Verschilt de gemiddelde lengte in de populatie van 170cm? b) Bedraagt het gemiddeld gewicht in de populatie meer dan 65kg? 5.2.2 Een grootschalig onderzoek werd uitgevoerd naar cannabisgebruik bij 18- tot 25-jarigen. “cannabis.txt” bevat voor iedere deelnemer het gerapporteerd aantal gram cannabis dat hij of zij wekelijks gebruikt. Toets of personen in deze leeftijdsgroep wekelijks gemiddeld minder dan een halve gram cannabis gebruiken.

52

5.3

Vergelijken van twee gemiddelden: onafhankelijke steekproeven (betweensubjects vergelijking)

Vergelijken van twee gemiddelden. t.test ( x, y, alternative = "two.sided", paired = FALSE)


De t-test voor twee onafhankelijke groepen:

Vector. Vector. Alternatieve hypothese, andere mogelijke waarden zijn “less” en “greater”. Afhankelijke steekproeven?

Volgens het afvalverwerkingsbedrijf is tussen 2011 en 2011 de totale hoeveelheid opgehaald afval per inwoner gedaald. We namen zowel in 2011 als in 2012 een steekproef (niet noodzakelijk bij dezelfde mensen) en noteerden de data in “afval.txt”. > data <- read.table("afval.txt", header = TRUE) > summary(data) afval jaar Min. :326.3 Min. :2011 1st Qu.:484.3 1st Qu.:2011 Median :522.8 Median :2011 Mean :521.0 Mean :2011 3rd Qu.:558.4 3rd Qu.:2012 Max. :688.2 Max. :2012

Data in groepen splitsen. split ( x, f)


We zien dat we een variabele hebben die het aantal kilogram afval bevat, en een variabele met het jaar (2011 of 2012). Iedere lijn in dit databestand komt dus overeen met één meting. Voor de t-test moeten we echter één vector geven die het aantal kilogram in 2011 weergeeft, en één voor 2012. Onderstaande functie kan hierbij van pas komen:

Te splitsen vector. Vector die de groep weergeeft.

> gesplitst <- split(data$afval, data$jaar) > gesplitst $‘2011‘ [1] 360.0 402.6 412.3 415.7 418.8 424.5 427.8 429.1 430.3 431.2 434.9 . . . [320] 625.3 625.5 627.5 628.4 629.6 631.4 633.6 639.3 642.3 $‘2012‘ [1] 326.3 331.7 372.0 386.5 387.9 396.9 397.7 405.3 405.7 407.4 409.4 . . . [287] 628.5 629.6 637.6 640.8 644.0 644.6 645.8 647.4 655.9 657.2 688.2 Merk op dat er in 2011 niet evenveel gegevens werden verzameld als in 2012, dat is ook niet nodig. We kunnen nu hieruit de twee vectoren extraheren die we nodig zullen hebben op basis van de namen 53

‘2011‘ en ‘2012‘ die we zien staan (de accentjes horen bij de naam van de variabele): > afval2011 <- gesplitst$‘2011‘ > afval2012 <- gesplitst$‘2012‘ De alternatieve hypothese stelt dat er een daling is dus Ha : µ2011 > µ2012 22 : > t.test(afval2011, afval2012, alternative = "greater") Welch Two Sample t-test data: afval2011 and afval2012 t = 2.0154, df = 559.61, p-value = 0.02217 alternative hypothesis: true difference in means is greater than 0 95 percent confidence interval: 1.646402 Inf sample estimates: mean of x mean of y 525.3192 516.2990 We verwerpen de nulhypothese, er werd minder afval geproduceerd in 2012 in vergelijking met 2011 (t(559) = 2.02, p = .022)23 . Oefeningen 5.3.1 Volgens een onderzoeker hebben mannen na een avondje uitgaan gemiddeld meer alcohol in hun bloed dan vrouwen. Een alcoholcontrole leverde volgende resultaten op (uitgedrukt in promille alcohol): 0.39 0.88 0.57 1.13

0.94 0.68 1.05 0.93

Mannen 0.24 0.69 0.91 0.87 0.59 0.90 0.46

0.73 0.35 1.61

0.61 0.11 0.87 0.17

0.54 0.23 0.00

Vrouwen 0.52 1.03 0.41 0.70 1.35 1.13

1.08 1.87 0.75

Toets de hypothese van de onderzoeker. 5.3.2 Is er een verschil in gemiddeld alcoholgehalte tussen Belgische en Duitse pils op basis van de gegevens in “pils.txt”?

5.4

Vergelijken van twee gemiddelden: gekoppelde data (paired samples, within-subject)

De t-test voor gekoppelde data verloopt in R gelijkaardig aan die voor onafhankelijke steekproeven. De twee vectoren zijn nu twee echter metingen bij dezelfde subjecten, en de parameter paired mag nu niet op zijn standaardwaarde staan. Een kleuterleidster noteert hoe vaak een kind in haar klas storend gedrag vertoont. Ze vindt dat dit veel te frequent gebeurt en besluit daarom operante conditionering toe te passen: bij iedere activiteit kan een kind een punt verdienen als het zich goed gedraagt. Wie 10 punten heeft krijgt een snoepje. Twee weken na invoering van dit systeem besluit de kleuterleidster het systeem te evalueren. Het databestand 22 Let

dus op met termen als “hoger”/“lager”, “meer”/“minder”, ... Schrijf je alternatieve hypothese altijd in de vorm <, >, of 6= en beslis op basis daarvan of het respectievelijk "less", "greater" of "two.sided" is. 23 Bij de t-test voor onafhankelijke populaties is het aantal vrijheidsgraden (meestal) geen geheel getal. Voor de rapportage nemen we steeds het eerste kleinere gehele getal (met andere woorden, laat de decimalen vallen).

54

is “kleuters.txt”. De kleuterleidster veronderstelt uiteraard dat het aantal storende gedragingen is gedaald, met andere woorden Ha : µvoor > µna . Bijgevolg wordt dit: > t.test(data$voor, data$na, paired = TRUE, alternative = "greater") Paired t-test data: data$voor and data$na t = 3.3514, df = 17, p-value = 0.001893 alternative hypothesis: true difference in means is greater than 0 95 percent confidence interval: 3.606944 Inf sample estimates: mean of the differences 7.5 We verwerpen bijgevolg de nulhypothese en besluiten dat het storend gedrag inderdaad is afgenomen (t(17) = 3.35, p = .002). Merk op dat we voor en na even goed hadden kunnen omdraaien: Ha : µna < µvoor . Dit zou resulteren in : > t.test(data$na, data$voor, paired = TRUE, alternative = "less") Wat uiteraard wel tot dezelfde conclusie zal leiden. Oefeningen 5.4.1 Dokter Hautakivi is psychiater op een afdeling waar patiënten met psychotische ziektebeelden verblijven. Van een farmaceutische firma krijgt hij een reis naar de Cara¨ıben aangeboden als hij kan aantonen dat hun nieuwe medicament het aantal psychotische symptomen dat een patiënt vertoont (hallucinaties, wanen, ...) kan doen dalen. Hij laat zijn patiënten observeren door verpleegkundigen die een hele dag lang noteren hoeveel psychotische symptomen de patiënten vertonen, en dit zowel op een dag wanneer de patiënt een placebo heeft gekregen als op een dag dat de patiënt het nieuwe medicijn heeft gekregen. Kan Dr. Hautakivi op vakantie? Het databestand heet “hautakivi.txt”. 5.4.2 Het bestand “IQpartners.txt” bevat IQ-scores van 138 koppels. a) Ga na of er een verschil is in gemiddeld intelligentiequotiënt. b) Het zijn hier geen twee metingen binnen één subject en toch spreken we hier van een within-subject opzet. Waarom?

55

Oplossingen 2.1.1 a) > 1 + 2 * 3 + 4 [1] 11 b) > (1 + 2) * (3 + 4) [1] 21 c) 1 + (2 * 3) + 4 [1] 11 2.1.2 > 2401^(1/4) [1] 7 2.1.3 > 1 + (1 / (1 + (1 / 2))) [1] 1.666667 2.1.4 a) > geboortejaar <- c(1980, 1973, 1992, 1990, 1987, 1985) b) Deze cursus werd geschreven in 2014, dus: > leeftijd <- 2014 - geboortejaar c) > mean(leeftijd) [1] 29.5 Indien je een ander jaartal hebt gebruikt (bijvoorbeeld 2015) zal de gemiddelde leeftijd natuurlijk een jaar hoger liggen. 2.1.5 > voor <- c(74, 60, 79, 72, 69, 60, 78, 69, 56, 63) > na <- c(75, 65, 78, 74, 69, 57, 80, 73, 60, 72) a) > mean(na) - mean(voor) [1] 2.3 b) > toename <- na - voor > mean(toename) [1] 2.3 2.2.1 > f <- function(x) 2 - (x - 2)^2 > integrate(f, 1, 3) 3.333333 with absolute error < 3.7e-14 > plot(f, 0, 4, xlab = , ylab = ) > lines(c(1, 3), c(0, 0), col = "purple") > lines(c(1, 1), c(0, f(1)), col = "purple") > lines(c(3, 3), c(0, f(3)), col = "purple") Niet de enige mogelijke manier: > > > >

lines(c(1, lines(c(1, lines(c(2, lines(c(1,

2), 2), 3), 3),

c(f(1), c(f(1), c(f(2), c(f(1),

f(2)), f(2)), f(3)), f(3)),

col col col col

= = = =

"purple") "purple") "purple") "purple")

Gezien f (1) = f (3): > (2 * f(1)) + ((2 * (f(2) - f(1))) / 2) [1] 3

56

2.3.1 > X <- matrix(c(18, 17, 18, 19, 18, 1, 1, 0, 0, 1, 16, 8, 4, 12, 12), 5, 3) > solve(t(X) %*% X) [,1] [,2] [,3] [1,] 0.004893964 0.006525285 -0.007748777 [2,] 0.006525285 1.092033714 -0.072831702 [3,] -0.007748777 -0.072831702 0.017477229 3.1.1 Vergeet niet je werkmap in te stellen indien nodig. Vervolgens: > neerslag <- read.table("neerslag.txt", header = TRUE) Als we neerslag opvragen zien we dat de variabele waarin we ge¨ınteresseerd zijn de naam “mm” heeft in de dataset: > mm <- neerslag$mm > gemiddelde <- mean(mm) > s2 <- mean((mm - mean(mm))^2) Of: > n <- length(mm) > s2 <- var(mm) * (n - 1) / n De rest van de berekeningen: > > > > >

s <- sqrt(s2) m3 <- mean((mm - mean(mm))^3) m4 <- mean((mm - mean(mm))^4) g1 <- m3 / s^3 g2 <- m4 / s^4 - 3

Resultaten en interpretatie: > gemiddelde [1] 792.5358 > s2 [1] 14343.83 > g1 [1] -0.1023693 > g2 [1] 0.3003735 g1 < 0, de verdeling is dus links scheef. g2 > 0, de verdeling is dus leptokurtisch. Beide liggen echter zeer dicht bij nul, dus echt duidelijk gaan we dit niet zien op het histogram. hist(mm) bevestigt dit. 3.1.2 > les <- read.table("les.txt", header = TRUE) a) > kwaliteit <- les$Leskwaliteit > inhoud <- les$Inhoud > boxplot(kwaliteit, inhoud, names = c("Leskwaliteit", "Inhoud")) De leskwaliteit wordt hoger gewaardeerd dan de inhoud van de les (mediaan). De spreiding is even groot (variatiebreedte en interkwartiel), maar de leskwaliteit werd duidelijk minder symmetrisch gescoord dan de inhoud. b) Nu moeten we de variabele kwaliteit opsplitsen naar gelang het geslacht, dus wijzen we eerst geslacht toe: > geslacht <- les$Geslacht We gebruiken table(vector) om een frequentietabel te bekomen van deze variabele: > table(geslacht) geslacht Man Vrouw

57

34

37

Dit leert ons dat de mogelijke waarden van deze variable "Man" en "Vrouw" zijn (dit had even goed 0 en 1, "man" en "vrouw", "V" en "M", ... kunnen zijn). > kwali m <- kwaliteit[geslacht == "Man"] > kwali v <- kwaliteit[geslacht == "Vrouw"] > boxplot(kwali m, kwali v, names = c("Leskwaliteit (M)", "Leskwaliteit (V)")) De vrouwen zijn duidelijk veel meer tevreden over de leskwaliteit dan de mannen (mediaan). Bovendien zijn zij het over het algemeen meer met elkaar eens dan de mannen (veel kleiner interkwartiel bij de vrouwen). Er is één uitzondering: een vrouw die duidelijk ontevreden is over de leskwaliteit. 3.1.3 > kleur <- read.table("kleur.txt", header = TRUE) > lievelings <- kleur$Lievelingskleur a) > length(lievelings) [1] 199 b) > table(lievelings) lievelings Blauw Geel Groen Oranje 28 14 33 64

Paars 39

Rood 21

Paars 39

Rood 21

Of: > summary(lievelings) Blauw Geel Groen Oranje 28 14 33 64 c) > barplot(table(lievelings)) 3.1.4 > werknemers <- read.table("voltijdsdeeltijds.txt", header = TRUE) > uren <- werknemers$Uur > statuut <- werknemers$Statuut a) > summary(uren) Min. 1st Qu. 18.00 30.00 > 39 - 30 [1] 9

Median 37.00

Mean 3rd Qu. 34.94 39.00

Max. 60.00

Aangezien summary niet meer is dan een (gelabelde) vector, kan je ook dit doen: > summary <- summary(uren) > Q1 <- summary[2] > Q3 <- summary[5] > Q3 - Q1 3rd Qu. 9 Dit kan handig zijn als je een bepaalde parameter later nog nodig hebt. b) Als je een frequentietabel hebt opgevraagd heb je gezien dat het aantal uren dat iemand werkt in dit voorbeeld discreet is. De aangepaste barplot() geniet hier dan ook de voorkeur boven hist: barplot(table(factor(uren, seq(from = min(uren), to = max(uren)))), space = 0) c) > > > > >

voltijds <- uren[statuut == "Voltijds"] barplot(table(factor(voltijds, min(voltijds):max(voltijds))), space = 0) s2 <- mean((voltijds - mean(voltijds))^2) s <- sqrt(s2) m3 <- mean((voltijds - mean(voltijds))^3)

58

> g1 <- m3 / s^3 > g1 [1] 3.550948 g1 is groter dan nul, de verdeling is dus rechts scheef. Het histogram toont dit bovendien duidelijk. Dit resultaat is hier niet eigenaardig doordat we enkel personen in beschouwing hebben genomen die voltijds werken. 3.1.5 > a1900 <- c(rep(0, 136), rep(1, 214), rep(2, 183), rep(3, 92), rep(4, 36)) > a2000 <- c(rep(0, 5), rep(1, 12), rep(2, 43), rep(3, 213), rep(4, 873)) > summary(a1900) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.000 1.000 1.000 1.513 2.000 4.000 > summary(a2000) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.00 4.00 4.00 3.69 4.00 4.00 Je kan ook nog een boxplot maken indien je dat duidelijker vindt. In het jaar 1900 kon de mediane 18-jarige enkel losse woorden lezen of schrijven, terwijl in 2000 de mediane 18-jarige een perfecte lezer/schrijver is. Minstens even opvallend is het verschil in spreiding: waar er in 1900 nog zeker wat spreiding is (interkwartiel, of zie boxplot), valt in 2000 het eerste kwartiel al samen met de maximumscore. In 2000 zit dus minsens 75% van de 18-jarigen aan de hoogste alfabetisme-score. Rond 1900 gingen de meeste kinderen slechts tot ongeveer 12 jaar naar school, om dan te beginnen werken op het land, in de mijnen, of in een fabriek. Enkel de gegoede burgerij kon zich permitteren om hun kinderen naar de middelbare school te laten gaan. Het alfabetisme ligt in 2000 echter niet enkel hoger dan in 1900 (meer kinderen gaan nu langer naar school), er is ook veel minder spreiding bij de 18-jarigen, wat te verklaren valt door de invoering van de leerplicht tot 18-jaar (vanaf 1914 was er in België leerplicht tot 12 jaar, in 1983 werd dit verlengd tot 18 jaar). 3.1.6 barplot(table(dobbel), space = 0) Hoe meer keren je met de dobbelsteen gooit, hoe homogener (gelijker) het resultaat, dus hoe minder de zes balkjes van elkaar zullen verschillen). 3.1.7 a) seq() genereert hier de getallen 1 tot en met 10, terwijl rep deze reeks hier 50 keer herhaalt. Bijgevolg geeft dit onderstaande frequentietabel: > table(data) data 1 2 3 4 5 6 7 8 9 10 50 50 50 50 50 50 50 50 50 50 b) > s2 <- mean((data > s <- sqrt(s2) > m3 <- mean((data > m4 <- mean((data > g1 <- m3 / s^3 > g2 <- m4 / s^4 > g1 [1] 0 > g2 [1] -1.224242 g1 = 0, de verdeling is

- mean(data))^2) - mean(data))^3) - mean(data))^4) 3

dus symmetrish. g2 < 0, de verdeling is dus platykurtisch.

c) barplot(table(data), space = 0) toont net als de frequentietabel dat iedere waarde even vaak voor komt. We hebben dus te maken met een homogene verdeling. Per definitie is die symmetrisch en platykurtisch.

59

3.1.8 a) > vec <- c(1, 3, 4, 5, 7, 7, 8, 9, 10, 11) b) > vec <- c(43, 48, 30, 20, 16, 99, 25, 52, 20, 27) c) > vec <- c(98, 98, 99, 99, 100, 100, 101, 101, 102, 102) d) > vec <- c(0, 2, 4, 4, 6, 8, 16, 18, 20, 22) 3.2.1 > taart <- c(52, 100, 80, 96, 120, 128) > oefeningen <- c(10, 7, 13, 13, 21, 14) a) > cor(taart, oefeningen) [1] 0.5129715 Er is dus een matig sterk positief lineair verband tussen hoeveel chocoladetaart men eet en het aantal oefeningen dat men maakt: Hoe meer chocoladetaart men eet, hoe meer oefeningen men dus maakt. b) > lm(oefeningen ~ taart) Call: lm(formula = oefeningen ~ taart) Coefficients: (Intercept) 4.63025

taart 0.08718

> plot(taart, oefeningen) > abline(4.63025, 0.08718) c) > 0.08718 * 150 + 4.63025 [1] 17.70725 d) > cor(taart, oefeningen)^2 [1] 0.2631398 Iets meer dan een kwart van de variantie in gemaakte oefeningen wordt verklaard door de hoeveelheid chocoladetaart die een student heeft gegeten. 3.2.2 > som afwijkingen <- sum((oefeningen - mean(oefeningen))^2) > voorspellingen <- 0.08718 * taart + 4.63025 > som voorspellingen <- sum((oefeningen - voorspellingen)^2) > (som afwijkingen - som voorspellingen) / som afwijkingen [1] 0.2631398 3.2.3 De correlatiecoëfficiënt van Pearson mag hier niet berekend worden omdat het een ordinale schaal betreft, wat wel mag: Kendalls τ : > cor(psy1, psy2, method = "kendall") [1] 0.7333333 Spearmans rs : > cor(psy1, psy2, method = "spearman") [1] 0.8285714 De sterk positieve rangcorrelaties wijzen er op dat beide psychiaters het in sterke mate met elkaar eens zijn. 3.2.4 > man <- c(45, 343)

60

> vrouw <- c(42, 357) > tabel <- rbind(man, vrouw) a) > lambda.test(tabel) $row [1] 0

$col [1] 0.007731959 De kolom is hier de verklarende variabele, dus λ = 0.007732 b) > a > b > c > d > Q > Q [1]

<<<<<-

tabel[1, 1] tabel[1, 2] tabel[2, 1] tabel[2, 2] unname((a * d - b * c) / (a * d + b * c))

0.05444521

c) > chisq.test(tabel) Pearson’s Chi-squared test with Yates’ continuity correction data: tabel X-squared = 0.1337, df = 1, p-value = 0.7146 d) > sqrt(0.1337 / sum(tabel)) [1] 0.01303402 e) > sqrt(0.1337/ (0.1337 + sum(tabel))) [1] 0.01303291 3.2.5 a) > x <- seq(1, 20) > y <- c(1265, 1234, 1813, 1125, 930, 1050, 893, 753, 822, 836, 1638, 1915, 1539, 1318, 1246, 1238, 1393, 1136, 1298, 1045) > lm(y ~ x) Call: lm(formula = y ~ x) Coefficients: (Intercept) 1174.021

x 4.793

> plot(x, y) > abline(1174.021, 4.793, col = "red") b) > lm(y[1:10] ~ x[1:10]) Call:

lm(formula = y[1:10] ~ x[1:10])

Coefficients: 1492.80

(Intercept) -76.49

x[1:10]

> abline(1492.80, -76.49, col = "blue")

61

> lm(y[11:20] ~ x[11:20]) Call: lm(formula = y[11:20] ~ x[11:20]) Coefficients: (Intercept) 2452.58

x[11:20] -69.42

> abline(2452.58, -69.42, col = "green") c) Simpson paradox. 3.2.6 Er zijn uiteraard oneindig veel mogelijkheden, bijvoorbeeld: > x <- c(10, 20, 30, 40, 50, 60) > y <- c(175, 160, 70, 95, 80, 105) > cor(x, y) [1] -0.6976062 3.2.7 > klein <- c(138, 92, 17) > groot <- c(38, 39, 12) > terrein <- c(27, 53, 82) > tabel <- rbind(klein, groot, terrein) a) > lambda.test(tabel) $row [1] 0.1898305

$col [1] 0.2589641 De rij is hier de verklarende variabele, dus λ = 0.1898 b) > a > b > c > d > Q > Q [1]

<<<<<-

tabel[1, 1] tabel[1, 2] + tabel[1, 3] tabel[2, 1] tabel[2, 2] + tabel[2, 3] unname((a * d - b * c) / (a * d + b * c))

0.259034

c) > chisq.test(tabel) Pearson’s Chi-squared test data: tabel X-squared = 126.6409, df = 4, p-value < 2.2e-16 d) > sqrt(126.6409 / sum(tabel)) [1] 0.5042807 e) > sqrt(126.6409 / (126.6409 + sum(tabel))) [1] 0.4502687 4.1.1

62

a) > dbinom(10, 24, 0.5) [1] 0.1169 b) > pbinom(11, 24, 0.5) [1] 0.4194099 c) > qbinom(0.5, 24, 0.5) [1] 12 4.1.2 a) > integrate(dnorm, mean = 500, sd = 50, lower = 500, upper = Inf) 0.5 with absolute error < 2.4e-05 b) > integrate(dnorm, mean = 500, sd = 50, lower = 400, upper = 500) 0.4772499 with absolute error < 5.3e-15 c) > integrate(dnorm, mean = 500, sd = 50, lower = 400, upper = 600) 0.9544997 with absolute error < 1.8e-11 d) > integrate(dnorm, mean = 500, sd = 50, lower = 500, upper = 500) 0 with absolute error < 0 e) > qnorm(.95, mean = 500, sd = 50) [1] 582.2427 95% haalt een score van 582.2427 of lager. 4.1.3 > integrate(dnorm, mean = 50000, sd = 2850, lower = 55000, upper = Inf) 0.03969129 with absolute error < 0.00011 > plot(function(x) dnorm(x, mean = 50000, sd = 2850), from = 40000, to = 60000) > r <- seq(from = 55000, to = 60000, by = 10) > x <- c(55000, r, 60000) > y <- c(0, dnorm(r, mean = 50000, sd = 2850), 0) > polygon(x, y, col = "purple") 4.1.4 a) > 143.4 - (2.5 * 6.9) [1] 126.15 > 143.4 + (2.5 * 6.9) [1] 160.65 b) > pnorm(126.15, mean = 143.4, sd = 6.9) + (1 - pnorm(160.65, mean = 143.4, sd = 6.9)) [1] 0.01241933 Of: > 2 * pnorm(-2.5) [1] 0.01241933 c) > > > > > > > > > > > > >

links1 <- 143.4 - (3.5 * 6.9) links2 <- 143.4 - (2.5 * 6.9) rechts1 <- 143.4 + (2.5 * 6.9) rechts2 <- 143.4 + (3.5 * 6.9) plot(function(x) dnorm(x, mean = 143.4, sd = 6.9), from = links1, to = rechts2) r <- seq(from = links1, to = links2, by = 0.01) x <- c(links1, r, links2) y <- c(0, dnorm(r, mean = 143.4, sd = 6.9), 0) polygon(x, y, col = "pink") r <- seq(from = rechts1, to = rechts2, by = 0.01) x <- c(rechts1, r, rechts2) y <- c(0, dnorm(r, mean = 143.4, sd = 6.9), 0) polygon(x, y, col = "pink") 63

5.1.1 > data <- read.table("lengtegewicht.txt", header = TRUE) > t.test(data$lengte) One Sample t-test data: data$lengte t = 134.1762, df = 51, p-value < 2.2e-16 alternative hypothesis: true mean is not equal to 0 95 percent confidence interval: 164.7038 169.7074 sample estimates: mean of x 167.2056 > t.test(data$gewicht) One Sample t-test data: data$gewicht t = 51.4641, df = 51, p-value < 2.2e-16 alternative hypothesis: true mean is not equal to 0 95 percent confidence interval: 62.51743 67.59296 sample estimates: mean of x 65.05519 Beide intervallen zijn quasi even breed: > 169.7074 - 164.7038 [1] 5.0036 > 67.59296 - 62.51743 [1] 5.07553 Maar dat van gewicht is een klein beetje groter. Aangezien het betrouwbaarheidsniveau en ook de steekproefgrootte identiek waren voor beide variabelen zal het gewicht in kilogram dus een iets grotere variantie hebben dan de lengte in centimeter, vanwaar het grotere BI. > var(data$lengte) [1] 80.75213 > var(data$gewicht) [1] 83.09192 5.2.1 > data <- read.table("lengtegewicht.txt", header = TRUE) a) > t.test(data$lengte, mu = 170) One Sample t-test data: data$lengte t = -2.2424, df = 51, p-value = 0.02931 alternative hypothesis: true mean is not equal to 170 95 percent confidence interval: 164.7038 169.7074 sample estimates: mean of x 167.2056 64

We verwerpen de nulhypothese, de gemiddelde lengte in de populatie is niet gelijk aan 170cm (gemiddelde = 167.21, SD = 8.99, t(51) = −2.24, p = .029). b) > t.test(data$gewicht, mu = 65, alternative = "greater") One Sample t-test data: data$gewicht t = 0.0437, df = 51, p-value = 0.4827 alternative hypothesis: true mean is greater than 65 95 percent confidence interval: 62.93748 Inf sample estimates: mean of x 65.05519 We behouden de nulhypothese, we kunnen niet stellen dat het gemiddeld gewicht in de populatie hoger is dan 65kg (gemiddelde = 65.06, SD = 9.12, t(51) = 0.044, p = .483). 5.2.2 > t.test(data$cannabis, mu = 0.5, alternative = "less") One Sample t-test data: data$cannabis t = -0.8337, df = 5831, p-value = 0.2022 alternative hypothesis: true mean is less than 0.5 95 percent confidence interval: -Inf 0.5091788 sample estimates: mean of x 0.4905693 We behouden de nulhypothse, we kunnen niet stellen dat het gemiddeld wekelijks cannabisgebruik in deze leeftijdsgroep minder is dan een halve gram gemiddelde = 0.491, SD = 0.834, t(5831) = 0.83, p = .202). 5.3.1 > Mannen <- c(0.39, 0.94, ..., 0.46) > Vrouwen <- c(0.61, 0.54, ..., 0.17) > t.test(Mannen, Vrouwen, alternative = "greater") Welch Two Sample t-test data: Mannen and Vrouwen t = 0.4272, df = 25.321, p-value = 0.3364 alternative hypothesis: true difference in means is greater than 0 95 percent confidence interval: -0.1879014 Inf sample estimates: mean of x mean of y 0.7733333 0.7106250 De nulhypothese kan niet verworpen worden, we hebben geen statistische evidentie dat mannen na het uitgaan meer alcohol in bloed zouden hebben dan vrouwen (t(25) = 0.43, p = .336). 5.3.2

65

> data <- read.table("pils.txt", header = TRUE) > gesplitst <- split(data$Alcohol, data$Herkomst) > t.test(gesplitst$Belgie, gesplitst$Duitsland) Welch Two Sample t-test data: gesplitst$Belgie and gesplitst$Duitsland t(32.679) = 4.6758, p-value = 4.875e-05 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: 0.1144198 0.2908089 sample estimates: mean of x mean of y 5.055556 4.852941 We verwerpen de nulhypothese, het alcoholgehalte is inderdaad verschillend (t(32) = 4.68, p < .001). 5.4.1 > t.test(data$placebo, data$medicatie, paired = TRUE, alternative = "greater") Paired t-test data: data$placebo and data$medicatie t = -0.9043, df = 48, p-value = 0.8148 alternative hypothesis: true difference in means is greater than 0 95 percent confidence interval: -2.272107 Inf sample estimates: mean of the differences -0.7959184 We kunnen de nulhypothese niet verwerpen, de patiënten hebben niet minder symptomen als ze de medicatie nemen (t(48) = −0.90, p = .815). In de steekproef ligt het gemiddeld aantal symptomen trouwens hoger wanneer men medicatie neemt dan wanneer men dat niet doet, dus strikt genomen had je de test zelfs niet moeten uitvoeren. Dokter Hautakivi zal zijn vakantie dus zelf moeten bekostigen. 5.4.2 a) > data <- read.table("IQpartners.txt", header = TRUE) > t.test(data$man, data$vrouw, paired = TRUE) Paired t-test data: data$man and data$vrouw t = -1.1394, df = 137, p-value = 0.2565 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -3.907113 1.050592 sample estimates: mean of the differences -1.428261 We kunnen de nulhypothese niet verwerpen, er is dus geen verschil in IQ-scores tussen mannen en vrouwen uit één koppel (t(137) = −1.14, p = .257). b) Onze onderzoekseenheid is hier niet de persoon maar het koppel.

66

Teken een rechte op een grafische voorstelling. abline ( a, b, lwd = 1, col = "black")

Constante. Richtingscoëfficiënt. Lijndikte. Kleur. JorisPieters,PeterTheuns,AlainIsaac

Maak een barplot. barplot ( height, space = NULL,

Vector met waarden die de hoogte van de balkjes bevat. Spatie tussen de balkjes (standaard wordt een spatie gebruik, wat correct is op nominale en ordinale schaal). Een vector met benamingen bijhorende de balkjes. Toevoegen aan bestaande grafische voorstelling.

names.arg = NULL, add = FALSE)

body ( fun) JorisPieters,PeterTheuns,AlainIsaac

Functie.

boxplot ( x, names = NULL)

Data waarop de boxplot is gebaseerd. Vector met namen voor verschillende boxplots in één afbeelding. JorisPieters,PeterTheuns,AlainIsaac

Vector aanmaken.


Vergelijking van een functie.

Boxplot tonen.


Appendix A: Overzicht functies

c Elementen van de vector gescheiden door een komma.

Vectoren als colommen aan elkaar binden tot een datatabel.

chisq.test ( x)

Tweedimensionale tabel (matrix).

Voeg kolomnamen toe aan een tabel. colnames ( x)


Bereken de χ2 contingentie.



cbind ( ...)

Tabel. JorisPieters,PeterTheuns,AlainIsaac

Overzicht van beschikbare kleurnamen. colors ()

67


( x)

Cosinusfunctie.

Vector. Vector. Type correlatiecoëfficiënt: man”.

cos ( x)

Hoek in radialen.

cov ( x, y)

Vector. Vector. JorisPieters,PeterTheuns,AlainIsaac

Cumulatieve frequenties berekenen. cumsum ( x)

Vector.



Kansdichtheid χ2 -verdeling. dchisq ( x, df)


Kansdichtheid binomiaalverdeling. dbinom ( x, size, prob)


Uitdrukking. Naam van de variabele waarnaar we afleiden (tussen aanhalingstekens). JorisPieters,PeterTheuns,AlainIsaac


Kansdichtheid normaalverdeling. dnorm ( x, mean = 0, sd = 1)


Kansdichtheid F -verdeling.


Afgeleide van een uitdrukking bepalen. deriv ( expr, name)


Populatiecovariantie berekenen (schatter!).

df ( x, df1, df2)

”pearson”, ”kendall” of ”spear-


cor ( x, y, method = "pearson")


Bereken een (rang)correlatieco¨ effici¨ ent.


68


Kansdichtheid t-verdeling. dt ( x, df)


exp ( x)


Exponenti¨ ele functie ex berekenen.

Frequentietabel met frequentie nul voor ontbrekende meetwaarden. factor ( x, levels)


Getal of vector.

Teksteditor openen om een object aan te passen. fix ( x)


Vector. Vector met alle mogelijke waarden.

Naam van het object dat je wil aanpassen. JorisPieters,PeterTheuns,AlainIsaac

Wiskundige functie aanmaken.

Huidige werkmap tonen.


function ()

Help opvragen.


getwd ()

help ( topic)

hist ( x, freq = TRUE, nclass = NULL)


Histogram tonen.

Naam van een functie. Plaats tussen aanhalingstekens indien de parameter speciale tekens bevat.

Vector. Absolute frequenties? (anders relatieve) Aantal klassen; indien weggelaten wordt automatisch gekozen.

Functie. Ondergrens. Bovengrens.

Ga na of een object “missing” is. is.na ( x)


integrate ( f, lower, upper)


Integraal van een functie berekenen.

Object.

69

lambda.test ( table)

length ( x)


Een tabel met twee variabelen.

Lengte van een vector weergeven. Vector. JorisPieters,PeterTheuns,AlainIsaac

library ( package)

Naam van het package

lines ( x, y, col = "black")

Vector met x-coördinaten. Vector met y-coördinaten Kleur. JorisPieters,PeterTheuns,AlainIsaac

Lineaire regressie. lm ( formula)

Y -variabele ˜ X-variabele.

matrix ( data, nrow = 1, ncol = 1)

Vector. Aantal rijen. Aantal kolommen. JorisPieters,PeterTheuns,AlainIsaac

Hoogste waarde. max ( ...)

mean ( x)


Vector.

Rekenkundig gemiddelde berekenen.

min ( ...)


Matrix cre¨ eren op basis van een vector.

Vector. JorisPieters,PeterTheuns,AlainIsaac

Laagste waarde.


Lijnen tekenen op een grafische voorstelling.

Vector.

Cumulatieve kansdichtheid binomiaalverdeling. pbinom ( q, size, prob)


70


Package laden.


Lambda (λ) berekenen (package: rapport)


pf ( q, df1, df2)


Cumulatieve kansdichtheid F -verdeling.


Teken een curve van een wiskundige functie of een puntenwolk.

X-coördinaten of uitdrukking. Y -coördinaten indien x geen uitdrukking is, anders optioneel. Beginpunt. Eindpunt. Kleur. Label X-as. Label Y -as. Teken assen. Toevoegen aan bestaande grafische voorstelling.

Cumulatieve kansdichtheid normaalverdeling. pnorm ( q, mean = 0, sd = 1)


plot ( x, y, from = 0, to = 0, col = "black", xlab = NULL, ylab = NULL, axes = TRUE, add = FALSE)


Cumulatieve kansdichtheid t-verdeling. pt ( q, df)


Vector met x-coördinaten. Vector met y-coördinaten. Kleur.


71


Getal of vector van x-coördinaten. Getal of vector van y-coördinaten. Lijndikte. Kleur.

Grafische voorstelling van een veelterm weergeven. polygon ( x, y, col = "black")


Teken ´ e´ en of meer punten op een grafische voorstelling. points ( x, y, lwd = 1, col = "black")


pchisq ( q, df)


Cumulatieve kansdichtheid χ2 -verdeling.

qbinom ( p, size, prob)

Proportie. Aantal. Kans op succes.

qchisq ( p, df)


Waarde waaronder een gegeven proportie ligt bij een χ2 -verdeling. Proportie. Aantal vrijheidsgraden.

qf ( p, df1, df2)


Waarde waaronder een gegeven proportie ligt bij een F -verdeling.

qnorm ( p, mean = 0, sd = 1)

Waarde waaronder een gegeven proportie ligt bij een t-verdeling. qt ( p, df)

Proportie. Aantal vrijheidsgraden.

rbind ( ...)

rbinom ( n, size, prob)



Genereer random getallen uit een binomiaalverdeling.

rchisq ( n, df)


Aantal random waarden. Aantal. Kans op succes.

Genereer random getallen uit een χ2 -verdeling.

read.table ( file, header = FALSE)


Vectoren als rijen aan elkaar binden tot een datatabel.



Bestandsnaam (tussen aanhalingstekens). Eerste lijn bevat naam variabele

72


Proportie. Gemiddelde. Standaardafwijking.


Proportie. Aantal 1e vrijheidsgraden. Aantal 2e vrijheidsgraden.

Waarde waaronder een gegeven proportie ligt bij een normaalverdeling.

Databestand openen.


Waarde waaronder een gegeven proportie ligt bij een binomiaalverdeling.

rep ( x, times = 1)


Repliceer een object een welbepaald aantal keer. Object. Aantal keren herhalen.

rf ( n, df1, df2)


Genereer random getallen uit een F -verdeling.

Aantal random waarden. Aantal 1e vrijheidsgraden. Aantal 2e vrijheidsgraden.

rnorm ( n, mean = 0, sd = 1)

Aantal random waarden. Gemiddelde. Standaardafwijking. JorisPieters,PeterTheuns,AlainIsaac

Getal of vector. Aantal decimalen.

Voeg rijnamen toe aan een tabel. Tabel.

Genereer random getallen uit een t-verdeling. rt ( n, df)


Random een bepaald aantal getallen uit een vector kiezen. sample ( x, size, replace = FALSE)

Vector met elementen om uit te kiezen. Aantal trekkingen. Met teruglegging.

Eerste getal. Laatste getal. Verschil.

73


Vector.

Reeks getallen met een vast interval aanmaken. seq ( from = 1, to = 10, by = 1)


Populatiestandaardafwijking berekenen (schatter!). sd ( x)


rownames ( x)


round ( x, digits = 0)


Afronden.


Genereer random getallen uit een normaalverdeling.


sin ( x)

Hoek in radialen.

solve ( a)

Vierkante matrix. JorisPieters,PeterTheuns,AlainIsaac

Data in groepen splitsen. split ( x, f)

Te splitsen vector. Vector die de groep weergeeft. JorisPieters,PeterTheuns,AlainIsaac

Bereken de vierkantswortel. sqrt ( x)

stem ( x)

Getal of vector. JorisPieters,PeterTheuns,AlainIsaac

Maak een stamdiagram.


Inverse van een matrix berekenen.

Vector.

sum ( x)


Sommeer de elementen van een vector. Vector.

Vijf-getallen-samenvatting weergeven (frequentietabel voor schaal zonder cijfers). summary ( object)

Een object waarvoor een samenvatting gemaakt moet worden. JorisPieters,PeterTheuns,AlainIsaac

Transponeren van een vector of matrix. t

Vector of matrix.

Betrouwbaarheidsinterval voor het gemiddelde. Vector. Betrouwbaarheidsniveau.

Toets een hypothese omtrent µ0 . t.test ( x, alternative = "two.sided", mu = 0)


t.test ( x, conf.level = 0.95)


( x)

Vector. Alternatieve hypothese, andere mogelijke waarden zijn “less” en “greater”. Waarde die we wensen te toetsen.

74


Sinusfunctie.

Frequentietabel weergeven. table ( ...)

Vector. Vector. Alternatieve hypothese, andere mogelijke waarden zijn “less” en “greater”. Afhankelijke steekproeven? JorisPieters,PeterTheuns,AlainIsaac

t.test ( x, y, alternative = "two.sided", paired = FALSE)


Vergelijken van twee gemiddelden.

Eén of meer objecten.

unname ( obj)

Object. JorisPieters,PeterTheuns,AlainIsaac

Populatievariantie berekenen (schatter!). var ( x)


Verwijder de namen of labels van een object.

Vector.

75

Appendix B: Packages Indien een bepaalde package reeds ge¨ınstalleerd is kan je die laden door library(naam package) in te typen of via Packages → Load package.... Doe je het via het menu dan kan je meteen zien welke packages reeds ge¨ınstalleerd zijn. Indien een package dat je nodig hebt nog niet ge¨ınstalleerd is, ga dan naar Packages → Install package(s)... en selecteer dewelke die je wenst. Het kan zijn dat je eerst gevraagd zal worden om een “CRAN mirror” te kiezen. Dat zijn de servers waarop alle packages staan. Theoretisch kan je best een server dicht bij huis kiezen, maar indien het geen gigantisch grote packages betreft zal het verschil in downloadsnelheid niet merkbaar zijn. De werking van een package kan afhankelijk zijn van andere packages. Daarom is het mogelijk dat bij het installeren van een package er automatisch nog andere packages ge¨ınstalleerd worden waardoor het proces wat langer duurt. Met de instructie .libPaths() kan worden opgevraagd waar R de packages bijhoudt.

76

Appendix C: Databestanden Alle databestanden die gebruikt werden in deze cursus zitten als bijlage in dit document. Door hieronder op de bestandsnaam te dubbelklikken kan je ze meteen openen, maar voor gebruik in R is het gemakkelijker als je ze bewaart. In Acrobat Reader klik je hiervoor op het paperclip-symbool, en dan zie je opnieuw de lijst met bestandsnamen. Selecteer alle bestanden, klik op “Save” en kies een map naar keuze. Indien je een andere pdf-viewer gebruikt dan Acrobat Reader is het mogelijk dat die niet in staat is om met bijlagen te werken. In dat geval moet je de databestanden handmatig downloaden. afval.txt cannabis.txt erytrocyten.txt hautakivi.txt IQpartners.txt kleur.txt kleuters.txt lengtegewicht.txt les.txt neerslag.txt pils.txt studenten.txt voltijdsdeeltijds.txt

77

Inleiding tot R bij Statistiek voor de Gedragswetenschappen

Recommend Documents