Enkele vaardigheden in R Dit document is een aanvulling op de Handleiding R die te vinden is op de website bij het vak Biostatistiek 2. De Handleiding is een algemene introductie, in dit document worden enkele specifieke vaardigheden toegelicht die je bij dit vak regelmatig zult moeten gebruiken.
1
Random getallen genereren
(Zie ook de paragraaf “Kansverdelingen en Toevalsgetallen” uit de Handleiding R.) Regelmatig zul je random getallen moeten genereren uit een gegeven verdeling. Voor standaardverdelingen zijn hiervoor functies ingebouwd in R. Deze hebben allemaal dezelfde structuur: de naam van de functie is rDIST, waarbij DIST vervangen moet worden door de naam van de verdeling. Voorbeelden zijn gegeven in de volgende tabel: binom binomiaal pois poisson geom geometrisch unif uniform exp exponentieel norm normaal Deze functies hebben allemaal als eerste argument n, het aantal random getallen dat uit de betreffende verdeling getrokken moet worden. De output is een vector ter lengte n. Je moet bij al deze functies ook nog een of meerdere andere argumenten opgeven, namelijk de parameters die de verdeling vastleggen. Kijk voor details in de help van de betreffende functie (type op de command line bijvoorbeeld help(rbinom)). Voorbeeld: de volgende twee commando’s zijn equivalent: 1 2
randvec = rnorm(1000,3,2) randvec = rnorm(n=1000,mean=3,sd=2)
en hebben allebei als resultaat, dat de variabele randvec een vector is ter lengte 1000, waarin elk element een (onafhankelijke) trekking is uit de normale verdeling, met verwachting 3 en standaardafwijking 2 (dus let op: het derde argument is de standaardafwijking, niet de variantie!). Opmerking: Als je K random getallen wilt genereren uit de binomiale verdeling met parameters n en p, gebruik je de aanroep rbinom(n=K, size=n, prob=p). Het eerste argument n van de functie rbinom is dus niet de parameter van de binomiale verdeling!
2
Histogram maken
Het histogram van een rijtje getallen (een vector) geeft inzicht in de verdeling van die getallen. Als de getallen trekkingen zijn uit eenzelfde kansverdeling, dan geeft het histogram een beeld van de bijbehorende kansdichtheidsfunctie. (Voor een onderbouwing van deze bewering, zie het dictaat dat hoort bij dit vak.) 1
Het histogram van een vector x kun je tekenen met de functie hist. Deze functie heeft, zoals blijkt uit help(hist), vele argumenten. De argumenten in tabel 1 zijn echter de belangrijkste. Tabel 1: Argumenten voor een histogram x breaks
probability
De vector waarvan je het histogram wilt zien (verplicht) Geeft aan waar de grenzen tussen de verschillende bins (histogramcategorie¨en) moeten komen. Kan een getal zijn: breaks=20 geeft aan dat R moet proberen om 20 bins te maken van gelijke breedte. Kan ook een vector zijn, dit zijn dan de grenzen van de bins. Als je breaks als een vector specificeert, moeten alle gegevens binnen de aldus gedefinieerde bins vallen. De aanroep hist(c(0,1,2,3,4),breaks=c(1,2,3)) zal bijvoorbeeld een foutmelding geven. Een true/false waarde, waarmee je kunt aangeven of de verticale as geschaald moet worden zodat de totale oppervlakte van de balkjes in het histogram gelijk is aan 1. Standaard gebruikt hist de instelling probability=FALSE, dan geeft de verticale as aan hoeveel getallen uit de vector x er in de betreffende bin vallen. Als je een histogram maakt van trekkingen uit een kansverdeling, is de instelling probability=TRUE meestal zinvoller, omdat je dan de kansdichtheidsfunctie in hetzelfde plaatje kunt tekenen en deze met het histogram kunt vergelijken.
Een handigheidje: je mag namen van argumenten afkorten als het geen dubbelzinnigheid oplevert. Zo kun je het argument probability=TRUE ook schrijven als prob=TRUE. Voorbeeld: de volgende commando’s produceren een histogram van 500 random trekkingen uit de exponenti¨ele verdeling met parameter λ = 4. Er zijn 50 bins en de verticale as wordt geschaald zodat de oppervlakte 1 is. 1 2
randexp = rexp(500,rate=4) hist(randexp, breaks=50, prob=TRUE)
In de praktijk weet je niet welke verdeling goed bij je data past. Indien je een idee hebt uit welke kansverdeling de getallen zijn getrokken, is het vaak zinvol om de kansdichtheidsfunctie en het histogram in hetzelfde figuur weer te geven. De dichtheidsfunctie is dDIST (met DIST weer als boven; zie de Handleiding R). Als we zouden proberen te plotten met het commando plot, wordt het histogram weer verwijderd en wordt een nieuw assenstelsel getekend. Om de dichtheid weer te geven in hetzelfde assenstelsel als het histogram, geven we het low level commando lines. We krijgen dan als vervolg van het bovenstaande voorbeeld: 3 4 5
xpoints = seq(0,4,by=0.1) ypoints = dexp(xpoints,rate=4) lines(xpoints,ypoints)
2
3
Plots maken en be¨ınvloeden
Met de functie-aanroep plot(x,y) wordt een nieuw assenstelsel aangemaakt. De vectoren x en y, die even lang moeten zijn, worden hierin tegen elkaar geplot als losse punten (x[i],y[i]). Een uitgebreide uitleg over plotten is opgenomen in de Handleiding R. Let bij het lezen van die handleiding met name op het verschil tussen het high level commando plot en het low level commando lines. Er zijn diverse opties die je extra mee kunt geven aan de functie plot, waarmee je de vorm van de plot kunt aanpassen. Zie hiervoor help(plot), help(par) en de Handleiding R. Een aantal veelgebruikte opties worden ook weergegeven in tabel 2. Tabel 2: Opties voor grafieken xlab="tekst" ylab="tekst" main="titel" xlim=c(-2,5) type=..
lty=..
pch=..
schrijft tekst langs de x-as schrijft tekst langs de y-as schrijft de titel titel boven de plot stelt het bereik van de x-as in van −2 tot 5 (analoog voor ylim) vertelt hoe de plot eruit moet zien: "l" verbindt de punten met lijnstukken, "p" tekent alleen punten (standaard voor plot), "b" tekent beide [zie help(plot) voor meer mogelijkheden] regelt hoe een geplotte lijn wordt weergegeven. Mogelijke instellingen zijn "solid", "dashed", "dotted", "dotdash", "longdash", "twodash". regelt de weergave van geplotte punten. Zie help(points), sectie pch values, voor de mogelijkheden.
Deze opties kun je ook aan hist meegeven. In de help van hist kun je zien dat een aantal opties hier een standaardwaarde krijgt.
4
Simulatie
In de statistiek gebruik je simulaties meestal om de verdeling van de een of andere grootheid te onderzoeken. In de opgavenseries zul je af en toe een simulatiestudie moeten uitvoeren. Hier wordt aan de hand van een voorbeeld uitgelegd hoe je dat kunt aanpakken. Als voorbeeld zullen we de verdeling onderzoeken van de stochast Y , die gedefinieerd is als Y :=
n X
Xi ,
i=1
de som van n onafhankelijke stochastische grootheden Xi (i = 1, . . . , n), waarbij de Xi ’s alternatief verdeeld zijn met parameter p: er geldt P(Xi = 1) = p en P(Xi = 0) = 1 − p. Uit de kansrekening is bekend dat Y zelf binomiaal verdeeld is met parameters n en p (ga dit na), maar we zullen dit voor nu even vergeten en de verdeling van Y door simulatie onderzoeken. We geven hier eerst de code van een geschikt simulatieprogramma voor dit probleem, en daarna lichten we de code regel voor regel toe. Hierbij wordt steeds aangegeven in 3
hoeverre de aanpak algemeen toepasbaar is. Na de bespreking wordt dan nog eens het algemene schema van een simulatieprogramma gegeven. Het doel van het simulatieprogramma is om verschillende realisaties van Y te genereren en op basis daarvan de verdeling te onderzoeken. Een geschikt programma voor het boven geschetste probleem is als volgt: 1 2 3 4 5 6 7 8 9 10 11 12 13
n = 25 #aantal Xi die we optellen p = 0.3 #succeskans M = 1000 #aantal random iteraties resultaat = rep(0,M) for( i in 1:M ) { x = rbinom(n,1,p) y = sum(x) resultaat[i] = y } hist(resultaat,prob=TRUE,breaks=seq(-0.5,n+0.5,by=1)) xpunten = 0:n ypunten = dbinom(xpunten,n,p) lines(xpunten,ypunten,lty="dashed")
In de regels 1 en 2 van dit programma worden de parameters van het probleem in variabelen opgeslagen. Het is verplicht om getalwaarden toe te kennen aan de parameters, maar welke waarde dat is, kun je zelf kiezen. Hier zijn de waarden n = 25 en p = 0.3 gekozen. Het is handig om de waarden op te slaan in variabelen, zodat je later de waarde van de parameters gemakkelijk kunt veranderen. De kern van het programma bevindt zich in de regels 6 t/m 8. In regel 6 wordt de variabele x gevuld met een vector ter lengte n; het i’de element van deze vector bevat een realisatie van Xi . In deze toepassing worden de realisaties van de n alternatief(p) verdeelde stochasten gegenereerd door gebruik te maken van de functie rbinom: er worden n trekkingen gegenereerd uit de binomiale verdeling met parameters 1 en p, wat precies de alt(p)-verdeling is. (n is dus niet een parameter van deze binomiale verdeling.) P In regel 7 wordt de realisatie van Y = ni=1 Xi berekend uit de realisaties van de Xi ’s die zijn opgeslagen in de vector x. In regel 8 wordt deze opgeslagen in de vector resultaat, hierop komen we verderop terug. Het doel van de simulatie is om zeer veel realisaties van Y te genereren. In de regels 6 t/m 8 wordt steeds ´e´en realisatie gegenereerd, dus worden deze regels met de for lus M keer herhaald. In elke iteratie wordt de realisatie y opgeslagen in de i’de positie van de vector resultaat. Door deze aanpak is aan het einde elke realisatie beschikbaar. De variabele resultaat moet van te voren wel gedefinieerd worden, anders zal regel 8 een foutmelding geven. Dit gebeurt op regel 4. Op deze regel wordt tevens genoeg geheugenruimte gereserveerd voor alle waarden die we in de for lus gaan opslaan in resultaat (de aanroep rep(0,M) geeft een vector met M nullen). Zouden we niet van te voren deze geheugenruimte reserveren, dan zou R in elke iteratie van de for lus de vector resultaat ´e´en element langer maken en daarvoor opnieuw geheugenruimte moeten vrijmaken, wat relatief veel tijd kost – het programma wordt daardoor trager, wat vooral bij grote M te merken is. 4
Na de for lus bevat de vector resultaat M realisaties van de stochast Y . De verdeling van Y kan dan gevisualiseerd worden met behulp van een histogram. Hieroverheen kan dan ook de theoretisch bekende kansdichtheid van Y worden weergegeven. Dit gebeurt in de regels 10 t/m 13. (Merk op dat in regel 10 de breaks zo zijn gekozen dat er een aparte staaf in het histogram komt voor elke (integer)waarde die Y kan krijgen.) De vraag hoe de waarde van M te kiezen blijft over. Hoe groter M wordt gekozen, hoe beter het beeld is dat van de verdeling van Y verkregen wordt. Dus in beginsel moet M zo groot mogelijk worden genomen. De beperkende factor is de tijd die de computer nodig heeft om de simulatie-iteraties uit te voeren. Daarom moet M niet te groot worden gekozen. Het is dus de kunst om een goede waarde van M te kiezen; uiteraard hangt dit af van de toepassing. Begin bijvoorbeeld met M = 100 om je programma te testen. Daarna neem je M steeds groter, totdat je merkt dat het resultaat niet meer veel verandert. Samenvattend zul je dus in een simulatiestudie veelal de volgende structuur aanhouden: 1 2 3 4 5 6 7 8 9 10 11
#definitie variabelen M = 1000 #aantal simulatie-iteraties resultaat = rep(0,M) for( i in 1:M ) { # - genereer met een geschikte rDIST een realisatie van een # stochastische variabele # - bereken de grootheid waarvan je de verdeling wilt onderzoeken, # op basis van die realisatie resultaat[i] = grootheid_die_je_berekend_hebt } hist(resultaat,prob=TRUE)
5 5.1
Fitten Lineaire regressie
Met de functie lm kun je regressie toepassen voor een Lineair Model. Veronderstel dat we een n-tal waarnemingen (xi , yi ) hebben, waarvan de co¨ordinaten zijn opgeslagen in vectoren x en y (de waarnemingen zijn dus (x[i], y[i])). Met lm kun je een kleinstekwadraten aanpassing aan de data laten uitvoeren, voor een model Yi = f (xi , θ) + ei , waarbij de stochastische storingstermen ei onafhankelijk zijn, verwachting 0 hebben en gelijke (maar onbekende) variantie σ 2 . De functie f (xi , θ) moet lineair zijn in de parameters θ. Deze functie wordt in de aanroep van de R-functie lm gespecificeerd in het argument formula. De volgende voorbeelden dienen om het gebruik van formula toe te lichten. model Yi = α + βxi + ei Yi = α + β1 xi + β2 x2i + ei Yi = βxi + ei Yi = β1 cos(xi ) + β2 exp(wi ) + ei 5
formula y ~ x y ~ x + x^2 y ~ -1 + x y ~ -1 + cos(x) + exp(w)
Een aantal verdere opmerkingen: • De symbolen in formula kunnen ofwel wiskundige functies zijn, ofwel namen van variabelen. De vectoren x en y, en in het laatste voorbeeld ook w, moeten dus van te voren zijn gedefinieerd. • Standaard neemt R een intercept (constante term) op in het model. Als je dit niet wilt, voeg je -1 toe aan de instelling van formula (lees dit als ‘niet de intercept’). • Voor de meeste resultaten die de functie lm uitvoert, zijn geen verdere aannames nodig over de precieze verdeling van de meetfouten ei . Waar wel extra aannames nodig zijn, neemt R aan dat de verdeling normaal is. • Het model moet lineair zijn in de regressieparameters. Om een model te fitten dat niet lineair is in de parameters, moet je ofwel een transformatie uitvoeren om zo een lineair model te verkrijgen, ofwel de functie nls gebruiken om niet-lineaire regressie uit te voeren. Als je de functie lm aanroept vanaf de command line, wordt een deel van de output naar het scherm geschreven. Je kunt de output ook opvangen in een variabele, bijvoorbeeld door de aanroep z=lm(y~x). Dan bevat z een lijst met allerlei resultaten van de regressie, die zijn geplaatst in de elementen van de lijst. Het element genaamd helementi van de lijst z kan worden opgeroepen met z$helementi. Nog meer informatie kun je krijgen met de aanroep sz=summary(z). De variabele sz bevat dan ook weer een lijst. Enkele elementen uit de lijsten z (de uitvoer van lm) en sz (de uitvoer van summary(z)), die je vaak zult gebruiken, zijn opgenomen in de volgende tabel. z$coefficients De kleinste-kwadratenschatters θˆ voor de regressieparameters θ ˆ z$residuals De ongeschaalde residuen yi − f (xi , θ) sz$coefficients Een matrix met daarin voor elke regressieparameter: de kleinstekwadratenschatting, de standaardfout, de waarde van de toetsingsgrootheid van de t-toets voor de nulhypothese dat de parameter gelijk is aan 0, en de bijbehorende tweezijdige overschrijdingskans √ ˆ waarbij wi het opgegeven sz$residuals De geschaalde residuen wi (yi −f (xi , θ)), gewicht is van het i-de meetpunt (zie onder; als je geen gewichten hebt opgegeven is wi = 1) 1 Pn sz$sigma Schatting voor de standaardafwijking σ ˆ: σ ˆ 2 = n−p i=1 wi (yi − 2 ˆ f (xi , θ)) (p is het aantal regressieparameters). De correlatiematrix van de kleinste-kwadratenschatters voor de regressieparameters, alleen beschikbaar als je de aanroep sz=summary(z,corr=T) doet. Je mag de namen van de elementen afkorten, bijvoorbeeld tot z$coef. De helptekst over lm (type help(lm)) beschrijft welke zaken de variabele z nog meer zal bevatten; help(summary.lm) beschrijft het object sz. (Een handigheidje: als je op de command line typt z$ en daarna ´e´en- of tweemaal op de tab-toets drukt, worden alle mogelijke aanvullingen van dit commando weergegeven. Dit aanvullen van een commando met de tab-toets werkt ook voor alle andere R-commando’s.) sz$corr
6
Uitbreidingen en variaties Wanneer bekend is dat de varianties van de verschillende meetfouten niet allemaal hetzelfde zijn, kan een gewogen regressie worden uitgevoerd. Beschouw opnieuw het model Yi = f (xi , θ) + ei , met f weer een lineaire functie van de parameters θ, en de ei onafhankelijk verdeeld met verwachting 0 en ongelijke varianties σi2 . De σi voldoen aan σi = si σ, met bekende si en onbekende σ. (De verhouding tussen de σi is dus bekend.) De kleinste kwadratenmethode wordt dan bijgesteld: er wordt een gewicht wi = 1/s2i aan het i’de meetpunt toegekend. De grootheid n X
2
wi (Yi − f (xi , θ))
=
i=1
n X Yi − f (xi , θ) 2 i=1
si
wordt geminimaliseerd naar θ. In R kan dit worden gedaan door het argument weights van de functie lm te specificeren. Heb je bijvoorbeeld vier metingen met bekende si die zijn opgeslagen in een vector svector, dan geef je een commando als lm( y~x, weights=c(1/svector^2) ). In het standaard regressiemodel wordt aangenomen dat de meetfouten alle dezelfde variantie σ 2 hebben, die een onbekende waarde heeft. (In de zojuist beschreven gewogen regressie is de waarde van σ ook onbekend.) Om de waarde van een uitdrukking waar σ in voor komt te kunnen berekenen, wordt de waarde van σ geschat met σ ˆ . In de situatie dat de variantie van de meetfouten wel bekend is, bijvoorbeeld door externe meting, is de procedure die R uitvoert dus niet correct. In deze situatie moet je de resultaten handmatig corrigeren. Daarvoor is het nodig om te weten hoe de resultaten afhangen van σ, en om te weten welke waarde R heeft berekend voor σ ˆ . Dit laatste kun je vinden in de uitvoer van de functie summary (zie boven). Wanneer de varianties van de meetfouten bekend zijn, kan ook een chikwadraattoets worden uitgevoerd om de ‘goodness of fit’ van de kleinste-kwadratenaanpassing te evalueren. De waarde van de toetsingsgrootheid wordt niet standaard door R berekend. Je kunt dit echter zelf doen door de som van de kwadraten van de residuen te berekenen met sum((z$residuals)^2), als z de output van lm bevat. (Bij een gewogen regressie sommeer je de kwadraten van de geschaalde residuen met de code sum((summary(z)$residuals)^2).) Met behulp van de functie pchisq kan dan de overschrijdingskans van de realisatie van U worden bepaald.
5.2
Niet-lineaire regressie
Om te fitten met een modelfunctie f (x, θ) die niet-lineair is in de parameters, gebruik je de R-functie nls (nonlinear least squares). Als argument voor deze functie moet ook weer een formula worden opgegeven, waarmee de modelfunctie wordt gespecificeerd. Ook de parameters moeten hierin worden opgenomen. Wil men bijvoorbeeld de functie p exp(−p x) aan de data fitten, dan gebruikt men formula=y~p*exp(-p*x). Opnieuw moeten de vectoren x en y van te voren zijn gedefinieerd. Bij de niet-lineaire kleinste-kwadratenmethode moeten startwaarden voor de parameters worden opgegeven. Dit doet men door het argument start op te geven. In het voorbeeldje hierboven kan het argument bijvoorbeeld worden start=list(p=0.5). Indien meerdere parameters zijn gebruikt in de modelfunctie, bijvoorbeeld alpha en xmin, dan worden de startwaarden bijvoorbeeld gespecificeerd als start=list(alpha=2.1,xmin=4). 7
De output van nls is niet, zoals wel bij lm, gemakkelijk te verwerken. Als je echter eerst de functie summary toepast op de output van nls, krijg je een lijst die je wel kunt lezen. Deze bevat bijvoorbeeld de parameter-matrix coefficients, en de vector residuals. Zie de help voor meer mogelijkheden.
8
