1. No category

Notater i ST2304
H. T. L
Innhold
1 Fordelingsfunksjonene i R
1.1 α-kvantilen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
3
2 Fisher test for ubalanserte modeller
2.1 Test mellom alternative modeller . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
4
3 AIC
4
4 Test av hvor godt en multinomisk modell passer dataene
4.1 Kontigenstabeller . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
5
5 Bevis for at D er kjikvadratfordelt
6
6 Feilforplantning med delta metoden
7
7 Generaliserte lineære modeller
7.1 Binomial respons – logit . . .
7.2 Normalfordelt – probit . . . .
7.3 Eksponentsialfordelt – cloglog
7.4 Poissonfordeling . . . . . . . .
.
.
.
.
7
7
8
8
9
8 Devians
8.1 Test av nøstede modeller basert på devians . . . . . . . . . . . . . . . . . . . . .
9
9
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
9 Over- og underdispersjon
10
9.1 Test for overdispersjon . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
9.2 Ekstra parameter for overdispersjon . . . . . . . . . . . . . . . . . . . . . . . . . 10
10 Numerisk maksimering av likelihood
11
10.1 Lage loglikelihoodfunksjon . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
10.2 Numerisk testing av ulike modeller . . . . . . . . . . . . . . . . . . . . . . . . . . 12
A Vedlegg A: Kule R-funksjoner
i
B Ting som må på regelarket
i
C R-kode til grunn for figurer
i
1
1
Fordelingsfunksjonene i R
For hver av de ulike sannsynlighetstetthetene har R fire tilhørende funksjoner: sannsynlighetstetthet, kummulativ funksjon, kvantilfunksjon og en funksjon som generer tilfeldige tall etter
fordelingen. Hver av disse funksjonene er spesifisert ved bokstavene, henholdsvis, d, p, q og r
forann et spesifikt ord som spesifiserer fordeling. Disse ordene er gitt i tabell 1. I figur 2 er den
kummulativ og tetthetsfunksjonen vist for normalfordelingen.
Tabell 1: Ulike sannsynlighetsfordelinger, deres funksjoner i R og eventuelle nødvendige ekstraargumenter. Ekstraargumenter i parantes er ikke
nødvendige, men svært vanlige.
Fordeling
Binomial
Hypergeometrisk
Geometrisk
Poisson
Eksponential
Normal
Kjikvadrat
Fisher
T-fordelingen
Weibull
Funksjon
-binom
-hyper
-geom
-pois
-exp
-norm
-chisq
-f
-t
-weilbull
Ekstraargument
size, prob
m, n, k
prob
lambda
(rate)
(mean, sd)
df
df1, df2
df
shape
Figur 1: Kummulativ funksjon og sannsynlighetstetthet for normalfordelingen med µ = 0 og σ = 1.
2
Figur 2: Kvantilfunksjon for normalfordelingen med µ = 0 og σ = 1. Som
du kan se er denne den inverse av den kummulative.
1.1
α-kvantilen
Hva som menes med α-kvantilen til en stokastisk variabel er litt forskjellig ulike steder. I R
defineres α-kvantilen alltid som nedre halvdel, mens i mange andre lærebøker er α-kvantilen
definert slik at sannsynlighetsmassen i øvre hale til fordelingen er lik α.
Ta f.eks. at vi skal utføre Fishers test for ubalanserte modeller. Her forkaster vi H0 for tilstrekkelige høye verdier for F , se likning 5. En fisherfordeling er plottet i figur 3 hvor det grå området
tilsvarer 5 % sannsynlighet for at testobservatoren vår, F , tar en slik verdi. Dette betyr at for
F = 2.35 er det 5 % sannsynlighet for at H0 er riktig, og altså velger vi å forkaste H0 for høyere
verdier for F enn dette.
Figur 3: Fisherfordelingsfunksjon med df1 = 10 og df2 = 20. Det grå
området på figuren tilsvarer en 5 % sannsynlighet for at H0 er sann og
grensen mellom det grå og det ikke-grå området er derfor forkastningsgrensa
for H0 på et 0.05-nivå.
3
2
Fisher test for ubalanserte modeller
Hvis det er ubalansert deisgn gjelder det ikke lenger at totalvariansen kan bli dekomponert
til komponenter som korresponderer til variasjon i ulike variabler. Men det gjelder fortsatt at
totatvariasjonen kan bli dekomponert til
n
X
SSDT =SSDmodell + SSDres
n
n
X
X
(Yi − Y )2 =
(ŷi − Y )2 +
(Yi − ŷi )2
i=1
i=1
(1)
(2)
i=1
hvor ŷi er forventet respons til den i’te variabelen basert på den estimerte modellen. Av dette
følger det at
SSDT
SSDmodell
SSDres
,
og
(3)
σ2
σ2
σ2
er kjkvadratfordelt med n − 1, p − 1 og n − p frihetsgrader under nullhypotesen der ingen av
forklaringsvariablene er tatt med:
H0 : Yi = µ + ei
(4)
som kan testes med den fisherfordelte testobservatoren
F =
SSDmodell
/(p − 1)
σ2
SSDres
/(n − p)
σ2
(5)
som har p − 1 og n − p frihetsgrader. H0 blir forkastet ved tilstrekkelige høye verdier for F .
Dette indikerer at en stor del av den totale variasjonen er forklart av modellen.
2.1
Test mellom alternative modeller
Når vi legger til en variabel til en modell vil residual sum of squares alltid synke. I slike tilfeller
kan H0 mot H1 testes med følgende testobservator
F =
SSDres,H0 −SSDres,H1
/(p1
σ2
SSDres
/(n − p1 )
σ2
− p0 )
(6)
som er F -fordelt med p1 − p0 og n − p1 frihetsgrader under nullhypotesen at den ekstra termen
har ingen effekt på responsen. Vi kan forkaste H0 for tilstrekkelige høye verdier av F .
Denne testen kan utføres med funksjonen drop1() i R. Denne brukes gjerne slik
drop1 ( f u l l m o d e l , t e s t=”F” )
Variabler kan også testes ved å legge til én og én med funksjonen add1() i R. Denne brukes
gjerne slik
add1 ( reduced , ˜+x1+x2 , t e s t=”F” )
3
AIC
En modell med flere variabler vil alltid passe bedre til dataene enn en modell med færre variabler.
Samtidig vil en modell med færre variabler ofte predikere dataene bedre enn en med mange
variabler. Dette er fordi variansen til en modell med mange variabler øker. Dette forholdet
mellom hvor god en modell er til å predikere data og hvor godt den passer kan bli estimert
4
med AIC (Akaikes Information Criterion). En lav AIC-verdi indikerer at modellen er god til å
predikere nye data. Det vil si at den passer de eksisterende dataene godt, samtidig som den ikke
er overtilpasset. AIC-verdien er gitt ved
AIC = 2k − 2 ln(L)
(7)
der k er lik antall parametre i funksjonen og L er maksimalverdien til likelihoodfunksjonen.
4
Test av hvor godt en multinomisk modell passer dataene
Anta at X1 , X2 . . . , Xk er multinomisk fordelt med parametre n, p1 , p2 , . . . , pk . Da gjelder følgende
Var(Xi ) = npi (1 − pi )
E(Xi ) = npi
Cov(Xi , Xj ) = −npi pj
(8)
og da er også
X − npi
p i
npi (1 − pi )
(9)
tilnærmet N(0,1) og vi har også at
D=
k
X
(Xi − npi )2
(10)
npi
i=1
er tilnærmet kjikvadratfordelt med k − 1 frihetsgrader gitt at alle E(Xi ) = npi ≥ 5.
Denne siste testobservatoren kan brukes til å tese nullhypotesen på formen
H0 : p1 = p1,0 , p2 = p2,0 , . . . , pk = pk,0
(11)
Hvor vi forkaster H0 hvis D > χ2α,k−1 . Dette kan gjøres med funksjonen chisq.test() i R som
ofte brukes på formen
c h i s q . t e s t ( c ( x1 , x2 , x3 ) , p=p )
Ofte vil vi i stedet teste om det er et spesielt matematisk forhold mellom pi ’ene. Dette vil si at
p1 , p2 , . . . , pk er funksjoner av et mindre antal parametre s, altså at pi = pi (θ1 , θ2 , . . . , θs ). Dette
kan testes med testobservatoren
D=
k
X
(Xi − np̂i )2
i=1
np̂i
der
p̂i = pi (θ̂1 , θ̂2 , . . . , θ̂s )
som er tilnærmet kjikvadratfordelt med k − 1 − s frihetsgrader gitt at np̂i > 5 for alle i.
4.1
Kontigenstabeller
Gitt at vi har en 2 × 2 kontigenstabell over 61 pasienter
Pirenzepine
Trithiozine
Total
Frisk
23
18
41
Så er sannsynligheten for hver kategori gitt ved
5
Syk
7
13
20
Totalt
30
31
61
(12)
Pirenzepine
Trithiozine
Margsannsynlighet
Her vil testobservatoren
Frisk
p11 = pq
p21 = p(1 − q)
p
Syk
p12 = (1 − p)q
p22 = (1 − p)(1 − q)
1−p
Margsannsynlighet
q
1−q
2 X
2
X
(Xi,j − np̂i,j )2
D=
np̂i,j
(13)
i=1 j=1
være kjikvadratfordelt med k−1−s = 1 frihetsgrader. Dette brukes i R i funksjonen chisq.test()
i den følgende koden
x<−matrix ( c ( 2 3 , 7 , 1 8 , 1 3 ) , 2 , 2 , byrow=T)
chisq . test (x)
Enhver r × c kontigenstabell vil ha k − 1 − s = (r − 1)(c − 1) frihetsgrader.
5
Bevis for at D er kjikvadratfordelt
Vi skal vise at
D=
k
X
(xi − npi )2
npi
i=1
(14)
er kjikvadratfordelt med k − 1 frihetsgrader for k = 2.
La oss anta at X1 , X2 , . . . , Xk er multinomisk fordelt med parametre p1 , p2 , . . . , pk og n.
For k = 2 vil det kun være to verdier for p og kun to verdier for X, dvs. at p2 = 1 − p1 og at
X2 = n − X1 . Ved å sette dette inn i uttrykket for D får vi at
D=
2
X
(xi − npi )2
i=1
npi
=
(x1 − np1 )2 (x2 − np2 )2
+
np1
np2
(x1 − np1 )2 ((n − x1 ) − n(1 − p1 ))2
=
+
np1
n(1 − p1 )
2
(x1 − np1 ) (1 − p1 ) + (np1 − x1 )2 p1
=
np1 (1 − p1 )
2
(x1 − np1 ) − (x1 − np1 )2 p1 + (np1 − x1 )2 p1
=
np1 (1 − p1 )
2
(x1 − np1 ) − (−1)2 (np1 − x1 )2 p1 + (np1 − x1 )2 p1
=
np1 (1 − p1 )
2
(x1 − np1 )
=
np1 (1 − p1 )
"
#2
x1 − np1
= p
np1 (1 − p1 )
(15)
(16)
(17)
(18)
(19)
(20)
(21)
Fra tidligere vet vi at for en binomisk fordeling med variablene X, p og n er Z tilnærmet
normalfordelt når
X − np
Z=p
(22)
np(1 − p)
6
Vi ser at ved å sette inn likning 22 i likning 21, får vi at
D = Z2
(23)
Vi vet at en sum med i ledd, der hvert ledd er kvadrert, av standardnormalfordelte variabler
kjikvadratfordelt med i frihetsgrader. Her er Z 2 en slik sum.
Altså er D kjikvadratfordelt med k − 1 frihetsgrader.
Vi får k − 1 frihetsgrader fordi D vi har 1 flere estimerte variabler i D enn antall ledd i summen.
Dermed må vi trekke fra 1.
6
Feilforplantning med delta metoden
Gitt at vi kjenner variansen til X og Y , kan variansen til Z = f (X, Y ) finnes med
Var Z ≈
∂f
∂x
2
Var X +
∂f
∂y
2
Var Y +
∂f
∂x
∂f
∂y
2Cov(X, Y )
(24)
Kovariansen kan finnes i R ved å bruke vcov(modell) på den tilpassede modellen.
7
Generaliserte lineære modeller
Grunnlaget for vanlig lineærregresjon er at vi antar at responsvariabelen er normalfordelt med
konstant varians lik σ 2 . I mange tilfeller stemmer ikke dette. En løsning på dette problemet
kan være å log-transformere forklaringsvariabelen, men dette er bare en delvis løsning fordi
forklaringsvariablene kan være ≤ 0, og log(x ≤ 0) er ikke definert. Generaliserte lineære modeller
er en løsning på dette problemet.
7.1
Binomial respons – logit
Anta at en responsvariabel er bionomialfordelt. I stedet for å jobbe med sannsynlighetsverdiene
0 ≤ p ≤ 1, kan vi jobbe med oddsen som er gitt ved
p
1−p
(25)
Fordelen er at vi ikke får noe øvre grense da 0 ≤ p/(1 − p) ≤ +∞. Ved å ta ln-verdien til oddsen
får vi utvidet dette gyldighetsintervallet til alle reelle verdier, −∞ ≤ ln(p/(1 − p)) ≤ +∞. Denne
nye funksjonen av p definerer vi som logit-funksjonen av p
logit p = ln
p
1−p
(26)
Vi antar at
logit p = β0 + β1 x1 + . . . + βk xk
(27)
Hvor høyre side er den lineære prediktoren, η. Ved å løse for p får vi at
p=
1
1
=
−η
−(β
+β
0
1 x1 +...+βk xk )
1+e
1+e
7
(28)
7.2
Normalfordelt – probit
Her antar vi at tidspunktet T da en hendelse intreffer er standardnormalfordelt med gjennomsnitt µ og standardavvik σ. Ut fra denne antakelsen vil sannsynligheten p for at en jente har
hatt menstruasjon ved alderen x være gitt ved
p = P (T ≤ x)
(29)
Da T er standardnormalfordelt med forventning µ og standarddavvik σ kan vi skrive om dette
til
T −µ
x−µ
x−µ
p=P
=φ
(30)
≤
σ
σ
σ
der φ er den kummulative standardnormalfordelingen. Likning 30 kan skrives om til
φ−1 (p) =
x−µ
σ
(31)
Vi definerer φ−1 () til probit(), og skriver om likning 31 til
probit p = β0 + β1 x
(32)
hvor vi har at
µ
1
β1 =
(33)
σ
σ
Likningene 33 kan omskrives slik at de kan benyttes til å estimere µ og σ på følgende måte
β0 = −
µ̂ = −
β̂0
σ̂ =
β̂1
1
β̂1
(34)
Denne modellen kan bli estimert i R ved å spesifisere den alternative probit-link funksjonen på
følgende måte
glm ( y˜x , b i n o m i a l ( l i n k=” p r o b i t ” ) )
7.3
Eksponentsialfordelt – cloglog
Anta at vi ønsker å modellere hvordan en hendelse intreffer over et tidsintervaller av ulik lengde
t og vi vil se hvordan sannsynligheten for at en hendelse intreffer basert på forklaringsvariablene
x1 , x2 , . . . , xk og lengden på tidsintervallet. Hvis raten av hendelsessintreffelse er λ og den er
konstant med hensyn på tid, vil tiden før intreffelse for et hvert tilfelle følge en eksponentialfordeling med parameter λ og sannsynligheten for at en hendelse skjer blir da
p = P (X ≤ 1) = 1 − P (X = 0) = 1 −
e−λt (λt)0
= 1 − e−λt
0!
(35)
Det blir også antatt at de ulike responsvariablene har en additiv effekt på logaritmen til λ
ln λ = β0 + β1 x1 + · · · + βk xk
(36)
Liking 35 kan skrives om til λ = −tln(1 − p), og ved å sette denne inn i 36, får vi
ln(−ln(1 − p)) = β0 + β1 x1 + · · · + βk xk + ln t
(37)
Ved å definere cloglog p = ln(−ln(1 − p)), får vi dermed følgende modell
cloglog p = β0 + β1 x1 + · · · + βk xk + ln t
8
(38)
Dette er altså en generalisert lineær modell med binomial respons og en cloglog link-funksjon
med en offset-term ln t. Offset-termen indikerer hvordan sannsynligheten p forandrer seg med
ulik lengde på tidsintervallene. På matematisk form kan vi skrive om modellen i likning 38 til
cloglog p = µ + αi + ln t
der λ = eµ+αi
(39)
Denne modellen kan bli spesifisert i R med følgende kode
glm ( y˜x , f a m i l y=b i n o m i a l ( l i n k=” c l o g l o g ” ) , o f f s e t=l o g ( t ) )
7.4
Poissonfordeling
Vi antar at
Yi ∼ pois(λi )
(40)
ln λi = β0 + β1 x1
(41)
og bruker modellen
8
Devians
Devians er et mål på hvor godt en modell passer dataene og spiller en liknende rolle som residual
sum of squares gjør for lineære modeller. Deviansen er definert som
D = 2(ln Lf ull − ln L)
(42)
hvor ln L er maksimal log-likelihood for tilpasset funksjon og ln Lf ull er maksimal log-likelihood
for den fulle modellen som har like mange parametre som det er observasjoner. D er tilnærmet
kjikvadratfordelt med n − p frihetsgrader.
Da deviansen til en binomialt fordelt responsvariabel er gitt så lenge gjennomsnittet er spesifisert, er det mulig å teste hvor godt en modell passer. Gitt H0 at dataene passer modellen er
deviansen D tilnærmet kjikvadratfordelt med n − p frihetsgrader. Hvis deviansen er stor nok
forkaster vi H0 . Deviansen kan finnes i R ved å bruke funksjonen summary(modell) og se på
“Residual deviance”.
8.1
Test av nøstede modeller basert på devians
Devians kan også brukes til å teste ulike nøstede modeller opp mot hverandre. Dette kan gjøres
ved å anta at H1 er en utvidelse av H0 ved å legge til parametre. La p1 og p0 være antallet
parameterestimater under H1 og H0 . Da blir forandringen i devians
D0 − D1
(43)
tilnærmet kjikvadratfordelt med p1 − p0 frihetsgrader under H0 . En slik test kan utføres ved
hjelp av R-koden
drop1 ( modell , t e s t=” Chisq ” )
En annen måte å teste det samme på er som følger. Gitt at H0 er den forenklede modellen og
H1 er den fulle vil da
D = 2(ln L1 − ln L0 )
(44)
være kjikvadratfordelt med p1 − p0 frihetsgrader. Dette kan altså testes med følgende R-kode
9
p c h i s q (D, d f =1, l o w e r . t a i l =FALSE)
Sagt på tabellspråk kan vi forkaste H0 når D er større enn den øvre α-kvantilen i kjikvadratfordelingen.
9
Over- og underdispersjon
Over- og underdispersjon er at den faktiske variansen i en forklaringsvariabel er henholdsvis
større eller mindre enn det som er antatt i modellen. Overdispersjon er vanlig i modeller hvor
observasjon av en hendelse øker sannsynligheten for å gjøre enda en observasjon. Et eksempel
på observasjoner hvor det vil forekomme overdispersjon er ved observasjon av et flokkdyr, her
vil observasjoner av et individ øke sannsynligheten for at man snart ser et nytt individ. Observasjonene er ikke uavhengige og variansen vil øke. Underdispersjon er det samme, bortsett fra
at den virkelige variansen er lavere enn den i modellen.
Problemet med over- og underdispersjon er altså at den faktiske variansen er annerledes enn
det modellen benytter. Dette kan løses ved å gange en faktor, ϕ, med variansen i modellen. Vi
får et estimat av ϕ ved å dele Residual deviance på antall frihetsgrader.
9.1
Test for overdispersjon
En stor devians kan være et tegn på overdispersjon. Dette testes med en kjikvadratfordeling
med n−p frihetsgrader under nullhypotesen at det ikke er noe overdispersjon. Gitt at vi har fått
oppgitt “Residual deviance”, kalt d med frihetsgrader df kan p-verdien beregnes på følgende
måte
p c h i s q ( d , d f=df , l o w e r . t a i l =F)
Sagt i tabellspråk: Vi forkaster altså H0 hvis deviansen er større enn den øvre α-kvantilen i
kjikvadratfordelingen.
9.2
Ekstra parameter for overdispersjon
Problemet med over- og underdispersjon er altså at den faktiske variansen er annerledes enn det
modellen benytter. Dette kan løses ved å gange en faktor, ϕ, med variansen i modellen. Denne
faktoren kan estimeres ved
D
ϕ̂ =
(45)
n−p
og denne faktoren kan testes med observatoren
D0 − D1
ϕ
(46)
som er omtrent kjikvadratfordelt med p1 − p0 frihetsgrader. I praksis må ϕ bli estimert og vi
bruker dermed i stedet
(D0 − D1 )/(p1 − P0 )
(47)
ϕ̂
som er tilnærmet F -fordelt med p1 − p0 og n − p1 frihetsgrader.
Forå tilpasse en modell med tilleggsparameteren ϕ brukes funksjonene family=quasibinomial()
eller family=quasipoisson() når modellen tilpasses med glm(). F -tester mellom ulike modeller fås ved å bruke drop1() eller add1() med argumentet test="F".
10
10
Numerisk maksimering av likelihood
For enkelte maksimeringer av likelihoodfunksjoner er det ikke mulig å finne parametrene algebraisk. Dette må derfor gjøres numerisk. I R brukes funksjonen optim() for å finne minimumspunktet i en funksjon. Arbeid med loglikelihood virker best, derfor gir vi negativ loglikelihood
til optim. optim() har også det tillegskravet at den vil ha de variablene den optimere som en
vektor, her gitt som p, og dette er det første argumentet i funksjonen. Gitt at vi har de observerte variablene x, og en funksjon som returnerer maksimal likelihood for vår modell lnL(p,x),
blir kallet til optim seende slik ut
optim ( par=c ( i n t 1 , i n t 2 ) , lnL , x=x )
der int1 og int2 er startverdiene for optimeringsfunksjonen.
Det kan ofte være nødvendig å spesifisere øvre og nedre grenser for optim. Når dette gjøres
brukes optimeringsteknikken "L-BFGS-B" og dette kan spesifiseres på følgende måte
optim ( c ( 1 , 1 ) , lnL , x=x , l o w e r=c ( 0 , 0 ) , upper=c ( I n f , I n f ) , method=”L−BFGS−B
”)
Standardfeilen til disse estimatene kan finnes fra kvadratroten av diagonalen av Hessian matrisen. Hessianmatrisen er en matrise av alle de andrederiverte. Disse må løses med funksjonen
solve(). Diagonalen vil gi deg de ulike variansene som du kan finne standardfeilen til. F.eks.
hvis du har estimert 3 parametre vil du få en 3 × 3 hessianmatrise hvor [1,1] gir første standardfeil, [2,2] gir nummer 2 og [3,3] gir nummer 3. Skrevet i R-kode får du standardfeilene ved
å skrive
f i t<−optim ( c ( 1 , 1 ) , lnL , x=x , h e s s i a n=TRUE)
sqrt ( diag ( solve ( f i t $ hessian ) ) )
Korrelasjonen kan finnes ved
cov2cor ( solve ( f i t $ hessian ) )
10.1
Lage loglikelihoodfunksjon
Log-likelihoodfunksjonen vi gir videre til optim() har vi sjeldent fra før av, og derfor må vi
gjerne skrive den selv. La oss si at vi har en generalisert lineær modell med en poissonfordelt
responsvariabel Yi som har faktoren λ lik
λi =
K
1 + ec(xi −x0 )
der c, K og x0 skal estimeres. Gitt λi beskrevet ovenfor blir log-likelihoodfunksjonen lik
yi
n
X
λi −λi
ln(L(c, K, x0 )) =
ln
e
yi !
i=1
Vi kan skrive dette uttrykket i R-kode, som følger kravene til optim(), på følgende måte
lnL<−f u n c t i o n ( par , x , y ) {
x0<−par [ 1 ]
c<−par [ 2 ]
K<−par [ 3 ]
lambda<−K/(1+ exp ( c ∗ ( x−x0 ) ) )
−sum ( d p o i s ( y , r a t e=lambda , l o g=TRUE) )
}
11
(48)
(49)
Siste linjen i koden er litt mystisk, men jeg forstår det slik: Yi er poissonfordelt med parameter λ. Likelihoodfunksjonen er jo gitt ved å gange sammen sannsynligheten for alle Y , og
denne er jo gitt av fordelingsfunksjonen med raten λ for alle y. log=TRUE gjør at vi får loglikelihoodfunksjonen og vi bruker sum() fordi å addere log-transformerte variabler er det samme
som å gange sammen utransformerte variabler. Minustegnet skyldes at optim() finner minimumspunkt og ikke maksimum.
10.2
Numerisk testing av ulike modeller
Når vi utvider en modell H0 til en mer generell H1 kan vi kjøre en statistisk test på denne
modellen ved å analysere maksimal likelihoodfunksjonene. Her må H0 vær et spesialtilfelle av
H1 . Vi har at dette forholdet
2(ln L1 − ln L0 )
(50)
er tilnærmet kjikvadratfordelt med p1 − p0 frihetsgrader.
12
A
Vedlegg A: Kule R-funksjoner
anova, step(),drop1(),add1(),predict,chisq.test(x,p=p),optim, glm, prop.test(), power.prop.test()
B
Ting som må på regelarket
1. Kravene bak de enkleste fordelingene
2. Hvilken vei noe skal være før det blir forkastet
3. Om man skal bruke lower.tail=TRUE/FALSE i ulike fordelinger
4. Husk også at du må repetere hvordan du bruker statistikktabellen
5. Eksempel på Hardy-Weinberg testing, s. 4 handout 2
6. Likelihoodfunksjon til logit p, s. 2 handout 4
C
R-kode til grunn for figurer
#s a n n s y n l i g h e t s t e t t h e t f u n k s j o n e r
x<−s e q ( − 4 , 4 , 0 . 0 1 )
p l o t ( −1 , −1 , ylim=c ( 0 , 1 ) , xlim=c ( min ( x ) ,max( x ) ) , main=” ” , y l a b=”
S a n n s y n l i g h e t ” , x l a b=”x” )
l i n e s ( x , dnorm ( x ) , c o l=” r e d ” )
l i n e s ( x , pnorm ( x ) , c o l=” b l u e ” )
#l i n e s ( x , qnorm ( pnorm ( x ) ) , c o l =”g r e e n ” )
l e g e n d ( −4 ,1 ,
c ( ”Kummulativ , p−f u n k s j o n ” , ” Tetthet , d−f u n k s j o n ” ) ,# p u t s
text in the legend
l t y=c ( 1 , 1 ) , #g i v e s t h e l e g e n d a p p r o p r i a t e symbols ( l i n e s )
lwd=c ( 2 . 5 , 2 . 5 ) , c o l=c ( ” b l u e ” , ” r e d ” ) )#g i v e s t h e l e g e n d l i n e s
t h e c o r r e c t c o l o r and width
#s a n n s y n l i g h e t s t e t t h e t f u n k s j o n e r
x<−s e q ( − 4 , 4 , 0 . 0 1 )
p l o t ( −1 , −1 , xlim=c ( 0 , 1 ) , ylim=c ( min ( x ) ,max( x ) ) , main=” ” , y l a b=”x” , x l a b=
” Sannsynlighet ”)
l i n e s ( pnorm ( x ) , qnorm ( pnorm ( x ) ) , c o l=” b l a c k ” )
legend (0 ,4 ,
c ( ” K van ti l , q−f u n k s j o n ” ) ,# p u t s t e x t i n t h e l e g e n d
l t y=c ( 1 ) , #g i v e s t h e l e g e n d a p p r o p r i a t e symbols ( l i n e s )
lwd=c ( 2 . 5 ) , c o l=c ( ” b l a c k ” ) )#g i v e s t h e l e g e n d l i n e s t h e
c o r r e c t c o l o r and width
#F o r k a s t n i n g med f i s h e r f o r d e l i n g
c u r v e ( d f ( x , d f 1 =10 , d f 2 =20) , from =0, t o =5,main=” F o r k a s t n i n g med
f i s h e r f o r d e l i n g e n ” ) #l a g e r kurven
n e d r e<−q f ( 0 . 0 5 , 1 0 , 2 0 , l o w e r . t a i l =F) #f i n n e r f o r k a s t n i n g s v e r d i
o v r e<−5
c o r d . x <− c ( nedre , s e q ( nedre , ovre , 0 . 0 1 ) , o v r e )
c o r d . y <− c ( 0 , d f ( s e q ( nedre , ovre , 0 . 0 1 ) , 1 0 , 2 0 ) , 0 )
polygon ( c o r d . x , c o r d . y , c o l= ’ g r e y ’ ) #t e g n e r g rå t t område
a x i s ( 1 , a t=nedre , round ( nedre , d i g i t s =2) ) #l e g g e r t i l 2 . 3 5 på aksen
i