Transcript Document
Statistik Brogaarden 20. og 21. januar 2014 1 Introduktion Mig Kasper K. Berthelsen, statistiker [email protected] people.math.aau.dk/~kkb Min arbejdsplads Institut for Matematiske Fag, AAU www.math.aau.dk Statistikgruppen 2 professorer, 8 lektorer, 2 adjunkter, 4-6 ph.d. studerende. 2 Matematiske uddannelser ved AAU To-fags-uddannelse (”Gymnasielæreruddanelsen”) Fx. matematik og fysik, matematik og dansk… Et-fags-uddannelse (”Anvendt matematik”) 4 år med matematik og 1 år med ”noget andet”, Ofte er det ene år på en ingeniør-retning. Matematik-Økonomi – 5-årig kandidatuddannelse Matematik-Teknologi – 5-årig kandidatuddannelse Seneste optage er 60-65 studerende Kandidatproduktion 5-10 studerende… må forventes at vokse 3 Tilbud til gymnasieelever Studerende for en dag Kontakt sekretær Lisbeth Grubbe [email protected] Numbers bloggen numb3rs.math.aau.dk Masser af inspiration i mere end 200 blog-indlæg Studiepraktik www.studiepraktik.aau.dk 4 Statistik – hvad er det? Statistik er kunsten at drage en generel konklusion på baggrund af ufuldstændig information. Det er uinteressant at finde ud af om der er (statistisk signifikant) forskel på udfaldet af to valg. Jeg kan jo bare tælle stemmerne…. Har jeg kun lavet en rundspørge, så er det straks mere interessant. I statistik består data af to komponenter En systematisk variation (”signal”) En tilfældig variation (”støj”) Statistik handler om at fjerne støj fra data. 5 Eksempel på (u)statistisk problemstilling Matematik B, Studentereksamen 2012 Ifølge hjemmesiden givblod.dk er fordelingen af blodtyper i den danske befolkning som følger: Tabellen nedenfor viser, hvorledes de 950 patienter i en bestemt lægeklinik fordeler sig på blodtyperne. Lægeklinikken vil undersøge nulhypotesen: Lægeklinikkens patienter har samme blodtypefordeling, som den danske befolkning. 6 Statistik: Religionskrig Der findes to hovedretninger indenfor statistik Frekventiel statistik Her er sandsynligheder baseret på frekvenser Den klassiske metode Bayesiansk statistik Baserer sig på subjektive sandsynligheder Moderne? Mere naturlig? Der har tidligere været en vældig krig mellem de to retninger… 7 Mange forskellige statistiske discipliner Survey/sprøgeskema Overlevelsesanalyse Longitudinelle analyser (gentagede målinger over tid) Tidsrækkeanalyse (fx aktiekurser) Rumlig statistik osv. osv. Anvendelsesområder: Biostatistik / biometri Økonometri Psykometri osv osv 8 Hvad så med mig…? Jeg beskæftiger mig med rumlig statistik. Mest punktprocesser: ”Støj” = placering Mindre geostatistik: ”Støj” = målte værdier 9 Hvad så med mig…? Simulationsbaseret inferens Særligt Markov kæde Monte Carlo Formål: Skabe en sekvens af afhængige stokastiske variable, der har de rette egenskaber i det lange løb… 10 En population og en stikprøve Population Konkret: Stemmeberettigede i Danmark Abstrakt: Alle målinger af lysets hastighed Stikprøve Vi udvælger (tilfældigt) elementer fra populationen. Kan gøres på mange måder. Ide: Vi vil gerne udtale os om hele populationen med udgangspunkt i stikprøven. Princip: Vi skal gøre dette med tanke på, at stikprøven er tilfældig. Vi kunne have været uheldige… 11 Stikprøveudtagning Tilfældig Alle elementer i populationen har lige stor sandsynlighed for at blive udvalgt. Eksempel: CPR-registret. Stratificeret Populationen inddeles i undergrupper (strata). Man udtager en tilfældig stikprøve fra hver gruppe. Nyttig metode, hvis man vil sikre at alle delgrupper er repræsenteret i stikprøven. Eksempel: Sammenligning af hjemløse og resten. 12 Stikprøveudtagning Cluster sampling. Man vælger et antal ”steder” i befolkningen og sampler i nærheden. Eksempel: Tilfældige veje i en kommune udvælges, hvorefter alle på de vej bliver spurgt. Problem: Systematiske fejl I den virkelige verden opstår mange statistiske fejl allerede i indsamlingsfasen… Ofte introduceres systematisk fejl – såkaldt bias. 13 The Literary Digest Poll (1936) Ikke-biased stikprøve Demokrater Republikanere Population Folk, der har telefon og/eller bil og/eller læser Digest. Demokrater Population Biased stikprøve Republikanere Ikke-biased, repræsentativ stikprøve fra hele populationen. Biased, ikke repræsentativ stikprøve af folk, der har telefon og/eller bil og/eller læser Digest. Literary Digest resultat: Alfred Landon slår Frenklin Roosevelt stort. Faktiske resultat: Landskredssejr til Roosevelt. Andre slags bias Formulering af spørgsmål har betydning for svar: In favour of new gasoline tax: 12% yes In favour of new gasoline tax to reduce US dependence on foreign oil: 55% yes In favour of new gasoline tax to reduce global warming: 59% yes 15 Andre slags bias Rækkefølgen af spørgsmål. Under den kolde krig blev følgende to spørgsmål stillet: A. B. Do you think the U.S. should let Russian newspaper reporters come here and report back whatever they want? Do you think Russia should let American newspaper reporters come in and report back whatever they want? 36% svarer ja til A. hvis det er første spørgsmål. 73% svarer ja til A. hvis det er andet spørgsmål. 16 Andre slags bias I en amerikansk undersøgelse afhang svarene i et telefoninterview af, den interviewedes forestilling om interviewerens etniske tilhørsforhold. Is the US society fair to everyone? 14% / 31% I medicinske forsøg: Alle der oplever bivirkninger dropper ud. Konklusion: Ingen bivirkninger… Medicinsk vs Kirurgisk behandling. Svage patienter udsættes ikke for den kirurgiske behandling. 17 Sandsynlighed: Opvarmning Udfald Resultatet af et ”eksperiment” kaldes et udfald. Eksempler: Eksperiment: Vælg en partileder / mål lysets hastighed Udfald: Lars / 299791 km/s (stikprøver fra hvilke populationer…?) Hændelse En hændelse er en mængde af udfald. Eksempler: Vælge en kvinde / Hastighedsmåling er ml. 299790 km/s og 299793 km/s 18 Sandsynlighed Sandsynlighed Sandsynligheden for en hændelse A er andelen af gange eksperimentet resulterer i hændelsen A i det lange løb. Notation P(A) betegner sandsynligheden for hændelsen A. Det behøver ikke være sådan: Hva’ nu hvis… P(Bayern München vinder CPL) P(Det regner i morgen) Eksempler på subjektive sandsynligheder 19 Sandsynlighed: Egenskaber og regneregler 1) 2) 3) 0 ≤ P(A) ≤ 1 P(A) = 0 - hændelsen A indtræffer aldrig. P(A) = 1 - hændelsen A indtræffer hver gang. P( ikke A) = 1 – P(A) Hvis A ikke indtræffer, så må ”ikke A” nødvendigvis indtræffe Hvis hændelserne A og B ikke kan indtræffe samtidigt gælder: P( A U B ) = P(A) + P(B) 20 Betinget sandsynlighed Betinget sandsynlighed Hvis A og B er mulige udfald, så gælder P B |A PA B PA Hvilket kan omskrives til multiplikationsreglen: P(A ∩ B) = P(A)P(B | A) 21 Uafhængighed Uafhængighed To hændelser A og B er uafhængige hvis og kun hvis P ( A B ) P ( A) P ( B ) hvilket kan omskrives til P( B | A ) = P(B) 22 Stokastisk variabel Stokastisk variabel Antag vi kan knytte en talværdi til hvert udfald af et eksperiment. Hvert eksperiment fører således til et tilfældigt tal. Dette tilfældige tal kaldes en stokastisk variabel. X 0 1 5 23 Diskret stokastisk variabel (SV) En stokastisk variabel X er diskret, hvis den kun kan tage et tælleligt antal værdier. Typisk 0, 1, 2, 3,… Lad P(k) betegne sandsynligheden for at den stokastiske variabel X tager værdien k. Dvs. P(1) = ”sandsynligheden for y tager værdien 1”. P(x) skal opfylde: 0 ≤ P(x) ≤ 1 for alle x. Salle xP(x) = 1 24 Kontinuert stokastisk variabel Hvis y er en kontinuert stokastisk variabel kan den tage alle værdier i et interval. Vi angiver sandsynligheden for at X falder i et interval [a ; b] ved et areal under en kurve. Tæthedsfunktion f(x) 0.0 0.1 0.2 0.3 0.4 P(1 ≤ X ≤ 2) = Areal -3 -2 -1 0 1 2 3 25 Tæthedsfunktionen 1) 2) (Sandsynligheds)Tæthedsfunktion f(x) f ( x) 0 for alle x f ( x ) dx 1 arealet under kurven f er 1 b 0.4 for a b 0.3 a f ( x ) dx 0.2 0.1 P (a x b) 0.0 3) -3 -2 -1 0 1 2 3 26 Middelværdi og varians for SV Middelværdi E[ X ] xP ( x ) xf ( x ) dx x 2 2 Varians V ( X ) E [( X ) ] 2 x ( x ) P ( x) 2 2 ( x ) f ( x ) dx 2 27 Normalfordelingen Normalfordelingen Klokkeformet og (fuldstændigt) karakteriseret ved middelværdi og standardafvigelse . Notation: x ~ N(,2) betyder at y er en kontinuert stokastisk variabel, der er normalfordelt med middelværdi og varians 2. Tæthedsfunktionen for normalfordelingen er f ( x) 1 2 2 x 2 exp 2 2 Egenskaber: Symmetrisk omkring f(x) > 0 for alle x. 95% 1.96 +1.96 28 Standardafvigelsen σ når X~N(μ,σ2) Cirka 68% af all observationer ligger indenfor en standard afvigelse fra middelværdien P ( X + ) 68 % Cirka 95% af alle observationer ligger indenfor 1.96 standard afvigelser fra middelværdien P ( 1 . 96 X + 1 . 96 ) 95 % Cirka 99.7% af alle observationer ligger indenfor 3 standard afvigelser fra middelværdien P ( 3 X + 3 ) 99 , 7 % Chebychevs ulighed (Tjebysjov?) Lad X være en stokastisk variabel med middelværdi og varians 2>0. Da gælder følgende ulighed P (| X | k ) 1 k 2 Eksempel: k = 2 P (| X | 2 ) 1 4 Dvs. sandsynligheden for at afvige mere end to standardafvigelser fra middelværdien er altid mindre en 25%. 30 Interessante størrelser I statistik optræder masser af stokastiske variable. Ofte er vi interesserede i at udregne en eller flere af følgende tre størrelser Middelværdi Variansen Sandsynlighed = E[X] 2 = Var(X) = E[(X-)2] P(X ∈ A) Dette kan være svært - eller umuligt. En løsning er (computer) simulationer. 31 Simpel Monte Carlo X stokastisk variabel ~ P Middelværdi: E[X] = Antag X1, X2,…, Xn ~ P og uafhængige. Udregn stikprøvegennemsnit X 1 n n Xi i 1 Da er X et Monte Carlo estimat af . Kan udvides til afhængige X1, X2,…, Xn Hvis X1, X2,… er en Markov kæde, så kaldes det MCMC. MCMC = Markov chain Monte Carlo. 32 Middelværdier er (næsten) alt! Antag vi er interesseret i en sandsynlighed . Vi har Bernoulli variabel X som vi kan simulere P(X = 1) = og P(X = 0) = 1 - Middelværdien: E [ X ] xP X x 0 (1 ) + 1 x Simuler uafh.: X1, X2,…, Xn ~ Ber Monte Carlo estimat af : X 1 n n Xi i 1 33 Den Centrale Grænseværdi Sætning (CLT) (Central limit theorem) Sætning: Lad X1, X2,…, Xn, er være n uafhængige stokastiske variable fra samme fordeling med middelværdi og varians 2. Da gælder, at når stikprøvestørrelsen n øges, så vil fordelingen af X N ( 0 ,1) n Tommelfingerregel: n = 30 er nok til en god tilnærmelse. Praktisk omskrivning: 2 X N , n Dvs. Monte Carlo estimatet er tilnærmet normalfordelt. Simulation i R 200 50 100 0 Man kan simulere en normalfordeling x = rnorm(n=1000, mean=198, sd=10) mean(x) ## stikprøvegennemsnit [1] 198.3801 sd(x) ## stikprøvestandardafvigelse Histogram of x [1] 10.02518 hist(x) ## histogram Frequency 170 180 190 200 210 220 230 x 35 Lidt flere plot i R 0.04 0.02 0.00 hist(x,freq=FALSE) curve(dnorm(x,mean=198,sd=10), 150,250,add=TRUE) Density Histogram of x 170 180 190 200 210 220 x 0.8 0.4 0.0 Den kumulerede fordelingsfunktion (F(x) = P(X≤x)) curve(pnorm(x,mean=198,sd=10), 150,250) pnorm(x, mean = 198, sd = 10) 160 180 200 220 240 x 36 Monte Carlo eksempel i R Antag X ~N(198,100) Find P(X<190) Simuler stikprøve fra normalfordeling: x = rnorm(n=1000,mean=198,sd=10) mean(x<190) ## falsk = 0 / sand = 1 [1] 0.207 Korrekte svar findes vha. pnorm(q=190,mean=198,sd=10) ## kumulerede ford. [1] 0.2118554 37 For at illustrere CLT kan vi nemt gentage Monte Carlo simulationen 2500 gange: > xbar = replicate(2500,mean(rnorm(n=1000,mean=198,sd=10)<190)) > hist(xbar,breaks=25,freq=FALSE) 38 c2-fordelingen Antag Z1,Z2,…,Zk er k uafhængige standard normalfordelte stokastiske variable, dvs. Zi ~ N(0,1). Definer k W Z 2 i Z 1 + Z 2 + ... + Z k 2 2 2 i 1 Da følger W en c2-fordeling med k frihedsgrader. Tæthedsfunktionen for c2-fordeling med k frihedsgrader er 39 antager kun positive værdier er højreskæv facon er givet ved antal frihedsgrader (df = degrees of freedom) har middelværdi = df og varians 2 = 2df. 0.15 0.10 df = 10 df = 10 0.05 c2-fordeling… df = 5 0.00 c2-fordelingen 0 10 20 30 40 40 R R er et open source statstikprogram og programmeringssprog introduceret i 1993. Seneste version er 3.0.2 R kan downloades på www.r-project.org R er i udgangspunktet uden peg-og-klik Mere end 2000 pakker (udvidelser a la et plugin) I det følgende tager vi udgangspunkt i Windows versionen. Der eksisterer versioner til Mac og Linux. For at få en smartere brugerflade anvender vi RStudio RStudio Sådan ser RStudio typisk ud første gang man starter det. Nederste vestre vindue er hvor man snakker direkte med R vha. tekst-kommandoer. RStudio – lidt opsætning Det er nyttigt at ændre R’s standard-mappe. Vælg Tools → Options Under ‘Default working directory..’ vælg den mappe hvor I vil gemme filer relateret til R (fx. data) R hjælp Man kan få hjælp vha. ?<kommando> > ?sum Man kan få RStudio til at gætte vha. Tab-knappen Man kan også søge efter hjælp vha. > help.search("plot") Statistisk test: Motivation (based on a true story…) Setup: Vi vil undersøge om der er sammenhæng mellem køn og om man gennemfører sit studie på normeret tid! Vi har spurgt 2000 (fiktive) AAU kandidater. Opsummering af data i en kontingenstabel Er der en sammenhæng ml. køn og rettidighed? 45 Hypotesetest Vi vil afgøre spørgsmålet vha. et såkaldt hypotesetest. Denne testteori er udvikle af Neyman og Pearson i 30erne. Grundlæggende ide: Vi inddeler verden i to Nul-hypotesen (H0) Alternativ-hypotesen (H1) (Der er ingen sammenhæng) (Der er en sammenhæng) Princip: Nul-hypotesen er sand indtil det modsatte er ”bevist”. Alle udregninger foretages under antagelse af H 0. Tvivlen skal altid komme nul-hypotesen til gode. 46 Forventede antal under H0 Nu tager vi udgangspunkt i at H0 er sand – dvs. ingen sammenhæng! Hvilke antal havde vi forventet, hvis H0 var sand? Uden sammenhæng burde andelen af ”På normeret tid” være den samme blandt mænd og kvinder, dvs. 64.25%. Da der er både 1000 mænd og kvinder ville vi forvente 642.5 ”rettidige” kandidater. 47 Observerede vs. Forventede antal En sammenligning af forventede og observerede antal: Observeret antal Forventet antal Er de observerede antal for langt fra de forventede til at vi kan tro på H0? 48 Observerede vs. Forventede antal Observeret antal Forventet antal O3 E3 Vi måler ”afstanden” mellem observerede og forventede vha. c 2 O i E i 2 i Ei O 1 E 1 2 E1 + O 2 E 2 2 E2 + O 3 E 3 2 E3 + O 4 E 4 2 E4 49 Resultat og Konklusion Resultat c 2 283 357 ,5 357 ,5 2 + 717 642 ,5 642 ,5 2 + 432 357 ,5 2 + 357 ,5 568 642 ,5 2 48 ,33 642 ,5 Så ”afstanden” en 48,33… er det for stor? Tommelfinger-regel (for 2x2 tabel): Afstanden er for stor hvis c2 > 3.84. Konklusion: Afstanden er for stor! Vi tror ikke på ”H0: Ingen sammenhæng”. Med andre ord: Der er en sammenhæng! Mere præcist: Der er en statistisk signifikant sammenhæng. 50 Forventede antal: Uafhængighed De forventede antal kan også udregnes med udgangspunkt i definitionen på uafhængige hændelser: Lad P(i,j) være sandsynligheden for at en tilfældig observation havner i celle (i,j). Lad N være det totale antal observationer. Det forventede antal i celle (i,j) er N· P(i,j). 51 Forventede antal: Uafhængighed Hvis der er uafhængighed så P(i,j) = P(i)P(j), dvs. det forventede antal er N· P(i)P(j). Antag Ni og Nj er antal observation i hhv. række i og kolonne j. Vi kan nu estimere de marginale sandsynligheder. Dvs. P(i) = Ni / N og P(j) = Nj / N Det forventede antal i celle (i,j) er da N· P(i)P(j) = N·P(i)(Ni / N)(Nj / N) = Ni Nj / N 52 Tabeller i R matrix(c(283,717,432,568),2,2,byrow=TRUE) [,1] [,2] [1,] 283 717 [2,] 432 568 as.table(matrix(c(283,717,432,568),2,2,byrow=TRUE)) A B A 283 717 B 432 568 tabel = as.table(matrix(c(283,717,432,568),2,2,byrow=TRUE)) addmargins(tabel) A B Sum A 283 717 1000 B 432 568 1000 Sum 715 1285 2000 53 c2-test i R c 2 O i Ei 1 2 2 i Ei > chisq.test(tabel) Pearson's Chi-squared test with Yates' continuity correction data: tabel X-squared = 47.6809, df = 1, p-value = 5.016e-12 > chisq.test(tabel)$statistic X-squared 47.68088 > chisq.test(tabel)$p.value [1] 5.015596e-12 54 Hvor kommer de 3,84 fra? For at forstå, hvor den kritiske værdi 3,84 kommer fra, så skal vi første forstå, hvor 2x2 tabellen kommer fra. A1 A2 B1 O11 O12 R1 B2 O21 O22 R2 C1 C2 n Der er n observationer, hvor Oij samtidigt tilhører population Aj og Bi. …men hvordan er data indsamlet? 55 Hvor kommer de 3,84 fra? Vores 2x2 tabel kan opstå på en af tre måder: 1. 2. 3. Vi har indsamlet n tilfældige observationer og kategoriseret dem mht. variablene A og B. Vi har indsamlet R1 observationer fra population B1 og tilsvarende indsamlet R2 observationer fra population B2. Indsamlet n observationer betinget af, at de marginale antal er givet ved R1, R2, C1 af C2. Asymptotisk, dvs. i grænsen når n vokser, så følger X2 i alle tre tilfælde en c2-fordelingen med 1 frihedsgrad. 56 Hvor kommer de 3,84 fra? For lille n er fordelingerne under de tre setup forskellige. 1. 2. 3. Vektoren (O11, O12, O21, O22) er multinomialfordelt, med sandsynligheder (p11, p12, p21, p22). Under H0 gælder p11/p12 = p21/p22. O11 og O21 er uafhængige og binomialfordelte med antal parametre R1 og R2. Under H0 har de samme sandsynlighedsparameter. Under H0 følger O11 en såkaldt hypergeometrisk fordeling. Husk: Asymptotisk, dvs. i grænsen når n vokser, så følger X2 i alle tre tilfælde en c2-fordelingen med 1 frihedsgrad. 57 Hypergeometriske fordeling En hypergeometrisk fordeling beskriver følgende situation: Sandsynligheden for k succeser i n trækninger uden tilbagelægning fra en endelig population af størrelse N med K succeser. Hvis X følger en hypergeometrisk fordeling (parametriseret som ovenfor), hvis P(X k) K n N K n k N n 58 Hvor kommer de 3,84 fra? Antag at H0 er sand – dvs. ingen sammenhæng! Vigtigt! Forestil jer, at vi nu udfører 10.000 nye undersøgelser af sammenhængen mellem køn og studietid. Trick: ”Omrøring” Svarer til situation 3! Hvor kommer de 3,84 fra? Antag at H0 er sand – dvs. ingen sammenhæng! Vigtigt! Forestil jer, at vi nu udfører 10.000 nye undersøgelser af sammenhængen mellem køn og studietid. Vi får 10.000 nye tabeller! Vi får 10.000 nye c2-værdier! De 3,84 er valgt, så for ki-i-anden 5% af tabeller overskrider c2-værdien 3,84. Konsekvens? 0 1000 3000 5000 7000 0 5 10 Her forkaster vi H0. 15 60 c2-fordelingen Histogrammet passer med en c2-fordeling (med én frihedsgrad): 0.3 4000 0.4 6000 0.5 ki-i-anden 0.2 Areal 5% 0.0 0.1 2000 0 0 5 Forkast 10 H0. 15 0 5 Forkast H0. 10 15 61 Omrøring i R ## Genskab data ud fra tabel data = TableData(2,2,c(283,717,432,568)) ## Nulstil variabel X2 X2 = 0 ## Omrøring og efterfølgende test 1000 gange for(i in 1:1000){ data$A = sample(data$A) data$B = sample(data$B) tabel = table(data) X2[i] = chisq.test(tabel)$statistic } 62 Grafisk opsummering hist(X2,freq=FALSE,breaks=50) curve(dchisq(x,df=1),0,25,add=TRUE) 0.8 0.4 0.0 Density Histogram of X2 0 5 10 15 > quantile(x=X2,probs=0.95)X2 95% 3.65922 63 Lidt asymptotik Et c2-test for uafhængighed i en r x c tabel er asymptotisk c2-fordelt med (r-1)(c-1) frihedsgrader. Asymptotikken består i at det totale antal observationer vokser. Tommelfingerregel: Hvis Ei≥5 for alle celler, så har asymptotikken ”slået til”. Dvs. hvis Ei≥5 for alle celler, så anvender vi en c2fordeling til at finde de kritiske værdier. 64 Type I fejl Hvis vi afviser H0, selvom den er sand, så har vi begået en Type I fejl. Sandsynligheden for at begå en Type I fejl betegnes signifikansnivauet og betegnes a. Når vi bruger 3,84 er så er sandsynligheden for at begå en Type I fejl (i princippet) 5%. 3,84 kaldes en kritisk værdi. Vi kan nemt undgå at begå en Type I fejl…! Hvordan? 65 Antag at H0 ikke passer! Antag at den sande fordeling af svar passer med følgende tabel: Mænd 40% 60% Mænd 33 68 Kvinder 60% 40% Kvinder 55 44 Vi skaber 10000 tabeller og 10000 c2-værdier 0 200 400 600 800 0.2242 0 Forkast 10 H020 . Nu begår vi en anden slags fejl…! 30 40 (0,2242) 66 Type II fejl Hvis vi ikke afviser H0, selvom den er falsk, så har vi begået en Type II fejl. I eksemplet er sandsynligheden for Type II fejl = 22.42%. Hvordan undgår man Type II fejl? Sandsynligheden for at afvise en falsk H0 kaldes styrken. - afhænger naturligvis af hvad ”sandheden” er. At få det korrekte signifikansniveau er ”nemt” – kunsten er at finde et test med en høj styrke. 67 Simulere H1 i R 0.00 0.04 0.08 ## Nul-stil variabel X2.H1 = 0 for(i in 1:10000){ ## Simuler tilfældig tabel med givne ssh’er data = rTableData(2,2,c(0.4,0.6,0.6,0.4),200) ## Udfør chi-i-anden test og gem X^2-værdi X2.H1[i] = chisq.test(table(data))$statistic } Histogram of X2.H1 ## Tegn histogram hist(X2.H1,freq=FALSE, breaks=50) Density 0 10 20 X2.H1 30 68 P-værdi En mere nuanceret metoder benytter P-værdien. P-værdien er sandsynligheden for at observere en mindst lige så ekstrem teststørrelse ”næste gang”, hvis nul-hypotesen er sand. Hvad der menes med ”mere ekstrem” afgøres af alternativhypotesen: Jo mere teststørrelsen peger i retning af alternativ-hypotesen jo mere ekstrem er den. I studenter-eksemplet: Jo større c2, jo mindre tror vi på H0. Vi afviser H0, hvis p-værdien < 5%. 69 c2-test vha. P-værdi P-værdien er sandsynligheden for mere kritiske værdier, hvis H0 er sand c2 , df = (r – 1)(c – 1) Grønne areal = P-værdien c2 Teststørrelsen, fx c2 = 48,33 P-værdien findes. vha. computer/lommeregner 70 Eksempel: Køn og rettidighed Teststørrelsen er c2 = 48,33 P-værdien er P = 0,0003. Konklusion: Da P-værdien er mindre end 0.05 afviser vi H0 Dvs. vi accepterer at der er en sammenhæng mellem køn og hvor rettidigt man afslutter sit studie. P-værdien 48,33 71 Fortolkning af P-værdi Nyttig og uformel fortolkning P-værdien er et udtryk for, hvor godt data passer med H0hypotesen. Jo mindre P-værdien er… jo mindre sandsynligt er det at observere end mindst lige så ekstrem værdi næste gang (hvis H0 er sand). jo mere ekstrem er teststørrelsen (hvis H0 er sand) jo mindre troværdig er H0. Beslutningsregel (hvis man ikke vil tænke selv): Vi afviser H0, hvis P-værdien er mindre end fx 5%. Andre bruger 1% eller 10% - kommer an på omstændighederne. 72 Dokumenterede misforståelser… Citater fra gymnasiebog i samfundsfag i forbindelse med c2test: Tekst der handler om P-værdien: ”Ved at gennemføre testen fås et resultat på 0,000144,…” ”Hvad betyder værdien på 0,000144? Værdien er sandsynligheden for at der ikke er sammenhæng. Der er således en sandsynlighed på 0,99856 for at sammenhængen mellem køn og partivalg er signifikant.” *Det skal understreges at andre allerede har gjort forlaget/forfatterne opmærksom på ovenstående fejl, og at der i den forbindelse er publiceret en rettelse på forlagets hjemmeside. 73 Pointer P-værdien er ikke sandsynligheden for at H0 er sand. Det er for simpelt at sammenligne P-værdien med fx 5%. Husk at tænke over, om vores antagelser opfyldt (fx Ei≥5). Spørgsmål 1: Hvis H0 er sand, hvilken fordeling følger P-værdien? Spørgsmål 2: Hvordan finder man en P-værdi uden en ”tabel”. 74 Fordelingen af P-værdien Hvis H0 er sand, så er P-værdien ligefordelt mellem 0 og 1. Simulation under H0 af 10.000 2x2 tabel med hver 10.000 svar. Histogram for de tilhørende 10.000 P-værdier: P-værdi 600 Tilsvarende histogram, hvis hver tabel består af kun 200 svar 0 200 200 600 400 1000 P-værdi 0 0.0 0.0 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 1.0 1.0 75 Monte Carlo P-værdi 0 68 Kvinder 55 44 ki-i-anden. MC P-værdi 0.0014 2000 4000 6000 8000 Ingrediens: Simulation af data under H0. Opskrift Simuler data under H0. For hvert simuleret data-sæt finder vi c2. Monte Carlo P-værdi Andelen af simulerede c2værdier > obs c2-værdi. 33 0 Mænd 5 10 15 Tabel P-værdi: 0.0027 (= 0.0014) 76 Monte Carlo p-værdi i R data = TableData(2,2,c(33,68,55,44)) tabel = table(data) X2.data = chisq.test(tabel)$statistic ## Teststørrelse for data X2 = 0 for(i in 1:10000){ ## 10.000 gange omrøring data$A = sample(data$A) data$B = sample(data$B) tabel = table(data) X2[i] = chisq.test(tabel)$statistic } 77 Ny Teststørrelse Med Monte Carlo kan vi konstruere nye teststørrelser. Fx Testsørrelsen er Maxi ( Ei – Oi ) H0 sand H0 falsk max(E-O) 0 0 1000 Kritisk værdi 500 1000 2000 3000 2000 max(E-O) 0 5 10 15 0 5 10 15 20 25 Sandsynlighed for Type II fejl 78 Pointer Med en Monte Carlo tilgang kan vi ”nemt” finde en P-værdi. Nemt at konstruere en ny teststørrelse og et nyt test med det korrekte signifikansniveau. Motivationer for at undersøge andre teststørrelser: Øge styrken? Resultat af en bedre model. 79 Mål for sammenhæng (a la korrelation) Race For Imod Total 600 Hvid 600 0 600 160 400 Sort 0 400 400 400 1000 Total 600 400 1000 Race For Imod Total Hvid 360 240 Sort 240 Total 600 P-værdi = 1 Forskel i andel hvide 360 600 240 0 . 60 0 . 60 0 . 0 Ingen sammenhæng 400 Forskel i andel hvide 600 600 P-værdi = 0 0 400 1 . 00 0 . 00 1 . 00 Stærk sammenhæng Forskellen i andele er et udtryk for styrken af sammenhængen mellem række og søjle variabel. 80 Mål for sammenhæng Eksempler på 2x2 tabeller og tilhørende mål for sammenhæng. 25 25 30 20 35 15 40 10 45 5 50 0 25 25 20 30 15 35 10 40 5 45 0 50 0.0 0.2 0.4 0.6 0.8 1.0 Bemærk: Målet er ikke symmetrisk: 80 20 80 40 40 10 20 10 0.3 0.0 81 Mars vs Snickers! Hvad kan du bedst li’? Mars eller Snickers? Vi har spurgt 100 mænd og 100m kvinder. En udregning viser, at X2 0,08 svarende til p-værdi = 0.8875 Mål for sammenhæng: 0.49-0.51 = -0.02 Spørgsmål: Er der en statistisk signifikant sammenhæng / forskel? Har den observerede forskel en praktisk betydning? 82 Mars vs Snickers! Flere folk! Vi spørger nogle flere! 200 mænd og 200 kvinder: Vi udregning viser at X2 = 0,16. Sammenh.: -0.02 Statistisk signifikant forskel? Praktisk forskel? 10.000 mænd og 10.000 kvinder: Vi udregning viser at X2 = 8,00. Konklusion? 83 Pointer P-værdien er ikke et udtryk for styrken af en sammenhæng. Der er en forskel på praktisk signifikans og statistisk signifikans. P-værdien er et udtryk for, hvor meget vi tror på H0 givet data. 84 Tilbage til Kandidat-undersøgelsen Kan det virkelig passe at kvinder er dårligere til at gennemføre på normeret tid??? Vi har faktisk spurgt om en ekstra ting: Hvilket fakultet har du studeret ved? TekNat (En af de teknisk naturvidenskabelige uddannelser) Samf (En af de samfundsvidenskabelige uddannelser) Vi laver nu en tabel for hvert af de to fakulteter! 85 Et Universitet - To Tabeller c2 0,024 c2 0,361 Konklusion? 86 Forklaret sammenhæng Konklusion: Om man gennemfører sit studie til tiden… Afhænger af fakultet! Afhænger ikke af køn! Hvorfor konkluderede vi først at køn har en betydning? TekNat-studier er nemme at gennemføre til tiden Samf-studier er svære at gennemføre til tiden. Piger foretrækker de svære studier. Drengene foretrækker de nemme studier. Og husk: Jeg har fundet på det hele… 87 Betinget uafhængighed To hændelser A og B er betinget uafhængige givet hændelsen C, hvis og kun hvis P( A B | C ) P( A | C )P(B | C ) Grafisk fremstilling (for variable) C A Fak B Køn Tid 88 Pointer Vi kan finde sammenhænge ml. to variable… …men i virkeligheden er der ingen sammenhæng, da de er betinget uafhængige! Ovenstående problem kan ikke løses af en statistiker Der er ingen indlysende måde, at ”opdage”, at man mangler en vigtig variabel. Løsningen kræver input fra en person med indsigt i den praktiske problemstilling. 89 Tabeller med få observationer Indtil nu har vi antaget at c2-fordelingen beskriver X2teststørrelsens fordeling godt. For 2x2 tabeller, hvor de observerede antal er små er c2fordelingen en dårlig tilnærmelse. I disse tilfælde er det almindeligt, at anvende Fishers eksakte test. Dette test baserer sig på sitution 3 fra tidligere hvor O11 følger en hypergeometrisk fordeling. P-værdien er summen af sandsynligheder for tabeller, der selv har en sandsynlighed, der er lig eller mindre den for den observerede tabel. 90 Fishers eksakte test 12 6 5 10 Data: Fordelingen af x2-teststørrelsen under H0. Sort streg: Sædvanlige kritiske værdi. Rød streg: X2-teststørrelsen for data. Bemærk: Forventede værdier er alle mindst 5. Med sædvanligt test ville det effektive signifikansniveau være 1,56%! 91 Fishers eksakte test i R > fisher.test(tabel) Fisher's Exact Test for Count Data data: tabel p-value = 0.08441 alternative hypothesis: true odds ratio is not equal to 1 ... > chisq.test(tabel, correct=FALSE) Pearson's Chi-squared test data: tabel X-squared = 2.4275, df = 1, p-value = 0.05642 92 En lidt større tabel Spørgsmål: Er der sammenhæng mellem køn og den måde man stemmer på? To variable: Køn: Mand / kvinde Partiforhold: Demokrat/ Uafhængig / Republikaner Vi er interesserede i fordelingen af stemmer, ikke de absolutte antal. 93 Relative fordeling Tabel over stemme fordelingen Stemmefordelingen blandt: Kvinder: Mænd: Alle: Vi ser at stemmefordelingen er forskellig Er forskellen statistisk signifikant? 94 c2-test vha. P-værdi P-værdien er sandsynligheden for mere kritiske værdier, hvis H0 er sand c2 , df = (r – 1)(c – 1) P-værdien c2 Teststørrelsen, fx c2 = 16.2 P-værdien findes. fx vha. software 95 Frihedsgrader Hvorfor har en 2x3 tabel 2 frihedsgrader? Antag vi kender alle række- og søjletotaler. Hvis vi kender antallet i bare to celler, så kan vi finde resten af antallene. Vi har frihed til at vælge to antal – derefter er resten givet! Partiforhold Demokrat Uafhængig Republikaner Total Kvinde 573 516 - 1511 Mand - - - 1260 959 991 821 2711 96 Residual: Motivation c2-testet kan afsløre, at data passer dårligt med nulhypotesen om statistisk uafhængighed. c2-testet siger intet om hvordan data passer dårligt. Det kunne fx være fordi: Et lille antal celler afviger meget. Et stort antal celler afviger lidt. Et residual siger noget om, hvor meget den enkelte celle afviger fra det forventede. 97 Residual Et (råt) residual for en celle er forskellen mellem Oi og Ei. Et standardiseret residual for en celle er zi Oi Ei se Oi Ei E i 1 rækkeandel 1 søjleandel Her er se standardfejlen, hvis H0 er sand. Dvs. det standardiserede residual måler antal se som forskellen mellem Oi og Ei afviger fra 0. z ligger omkring 0 med en standardafvigelse på 1. For store stikprøver er z ca. normalfordelt. 98 Residual: Eksempel For cellen ’Kvinde’ og ’Demokrat’ har vi z Oi Ei E i 1 rækkeandel 1 søjleandel 573 522 . 9 522 . 9 1 0 . 346 1 0 . 545 4 .0 Søjleandel: 1511/2771 0.545 Rækkeandel: 959/2771 0.346 99 0.3 0.2 0.1 Da z er cirka normalfordelt med middelværdi 0 og standardafvigelse 1, så er 4.0 ret ekstremt. I SPSS vælges ’Adjusted Standardized’ under ’Residuals’ 0.0 0.4 Residual: Eksempel fortsat -3 -2 -1 0 1 2 3 Det ses at det specielt er blandt demokrater, at afvigelsen mellem forventede og observerede værdier er stor. 100 Stadardiserede residualer i R > tabel = as.table(matrix(c(573,516,422,386,475,399),2,3, byrow=TRUE)) > tabel A B C A 573 516 422 B 386 475 399 > chisq.test(tabel)$stdres A B C A 4.015090 -1.940777 -2.145867 B -4.015090 1.940777 2.145867 101 Endnu en anvendelse af c2-testet Indkomstfordelingen for personer over 15 år i 2007 ifølge Danmark Statistik. I = indkomst i 1000kr I<50 50≤I <100 100≤I <150 150≤I <200 200≤I <300 300≤I <400 400≤I <500 500≤I % af befolkningen 6.4 9.3 17.8 12.3 24.3 18.0 6.6 5.3 I en markedsanalyse har vi spurgt 1000 personer i det lokale supermarked om deres mening om 3D-fjernsyn. Vi er nu kommet i tvivl om gruppen af adspurgte er repræsentativ. NB: Tåbelig måde at indsamle data på… 102 Repræsentativ stikprøve? Observerede antal (fra data) I = indkomst i 1000kr Antal i stikprøven 50≤I <100 100≤I <150 150≤I <200 200≤I <300 300≤I <400 400≤I <500 500≤I 98 88 199 136 210 179 52 38 Forventede antal (baseret på Danmarks Statistik) I = indkomst i 1000kr Antal i stikprøven I<50 I<50 50≤I <100 100≤I <150 150≤I <200 200≤I <300 300≤I <400 400≤I <500 500≤I 64 93 178 123 243 180 66 53 Igen måler vi afstanden vha. c 2 O i E i 2 i Ei .... Den kritiske værdi er 14,07 103 Test i R obs = c(98,88,199,136,210,179,52,38) p = c(0.064,0.093,0.178,0.123,0.243,0.180,0.066,0.053) chisq.test(x=obs,p=p) Chi-squared test for given probabilities data: obs X-squared = 33.8848, df = 7, p-value = 1.81e-05 Da p-værdi < 0.05 afviser vi H0. Kritisk værdi: > qchisq(0.95,df=7) [1] 14.06714 104 Goodness-of-Fit Resultat c 2 O i E i 2 i 33 ,88 Ei Da 33,87 > 14,07 så er afvigelsen for stor! Vi konkluderer derfor, at stikprøven ikke er repræsentativ. Samme undersøgelse, men i en ny butik: I = indkomst i 1000kr Antal i stikprøven I<50 50≤I <100 100≤I <150 150≤I <200 200≤I <300 300≤I <400 400≤I <500 500≤I 60 99 190 116 248 173 69 45 O i E i 2 3 , 56 Resultat c i Ei Er dette en repræsentativ stikprøve? 2 105 Testen Teststørrelsen er c2-fordelt med 7 frihedsgrader (antal grupper minus en). Hvis vi vælger et signifikansniveau på 5%, så er den kritiske værdi 14.07. Da vi har c2 = 33.87 > 14.07, så afviser vi påstanden om at stikprøven er repræsentativ. Hvis vi ikke havde afvist påstanden, så er det ikke det samme som, at stikprøven er repræsentativ… Under alle omstændigheder har vi kun forholdt os til om indkomstfordelingen matcher befolkningens. 106 Pointe Vi kan ikke bevise, at en stikprøve er repræsentativ. Vi kan ”bevise”, at stikprøven ikke er repræsentativ. 107 Hvorfor er c2-testen c2-fordelt? Goodness-of-fit for k=2 kategorier. Antag vi har n total antal observationer n1 antal observationer i kategori 1. p forventede andel observationer i kategori 1. I dette tilfælde er Goodness-of-fit teststørrelsen: X 2 ( n1 np ) 2 + ( n n1 ( n np )) n (1 p ) np (1 p )( n1 np ) np (1 p ) ( n1 np ) 2 2 + p ( np n1 ) 2 np (1 p ) 2 np (1 p ) 108 Hvorfor er c2-testen c2-fordelt? Bemærk, at n1 ~ B(n,p), dvs. E[n1] = np og Var[n1] = np(1-p) CLT giver n1 np np (1 p ) Dvs. X 2 ( n1 np ) N ( 0 ,1) 2 np (1 p ) c , df 1 2 109 Hva’ så hvis Ei < 5? 0.4 0.6 0.8 1.0 Histogram of X2 0.2 Antag p1 = 0.01, p2 = 0.99 og n = 100 Histogram for X2. Sande pi’er anvendt. 0.0 0 5 10 15 20 25 110 Asymptotik Chi-i-anden testen følger asymptotisk en c2-fordeling. Dvs. jo større n er, mere nærmer fordelingen sig en c2fordeling. Eksempel: Goodness-of-fit test med k=2 kategorier og p1 = p2 = ½. For n=16 ser fordelingen af X2 Kritisk værdi 3.84 P(x2≥3.84) = 0.076 Dvs. signifikansniveauet er 7.6% og ikke 5%. Bemærk: Ei = 8. 111 c2-test - generelt setup n observationer inddelt i k kategorier pi er sandsynligheden for i’te kategori pi er kendt Oi er det observerede antal i i’te kategori Ei = npi er det forventede antal observationer i i’te kategori. c2- teststørrelsen er X 2 O i E i 2 i Ei X2 er asymptotisk c2- teststørrelsen med k-1 frihedsgrader. Ei ≥ 5 for alle i er nok til at ”sikre asymptotikken”. 112 Goodness-of-fit for parametrisk fordeling Antag X1,…,Xn er uafhængige observationer fra en fordeling, der er specificeret ved s parametre q1,q2,…,qs. Antag Xi’erne kan inddeles i k kategorier. Lad pjq) betegne sandsynligheden for at Xi havner i j’te kategori. Hvis q er kendt er Ej = n pjq) og X2 følger en c2-fordeling med k-1 frihedsgrader. 113 Goodness-of-fit - ukendt q Setup kan betragtes som et multinomial-eksperiment, med sandsynligheder p1(q),…, pk(q) Antag qˆ1 , qˆ2 ,..., qˆs er maksimum likelihood estimater, opnået under multinomial-fordelingen. Antag Ei = npi(qˆ ) I dette tilfælde er X2 asymptotisk c2-fordelt med k-s-1 frihedsgrader. Estimeres q ud fra den oprindelige likelihood gælder c2k-s-1 ≤sd X2 ≤sd c2k-1 ,hvor ≤sd betyder ”stokastisk domineret” (A. Dasgupta, Asymptotic Theory of Statistics and Probability, 2008) 114 Eksempel: Hardy-Weinberg ligevægt Observationerne kan inddeles i k=3 kategorier. Med n observationer er de forventede antal E(AA) = n p2 E(Aa) = n 2p(1-p) E(aa) = n (1-p)2 p er ukendt. Et estimat af p udfra multinomial fordelingen er pˆ 2 O ( AA ) + O ( Aa ) 2 O ( AA ) + O ( Aa ) + O ( aa ) 2 O ( AA ) + O ( Aa ) 2n Antal frihedsgrader er således k-s-1 = 3-1-1 = 1 (og ikke 2) 115 Hardy-Weinberg Øverst: Ei udregnet vha. sande p. 0.0 Histogram for X^2, p kendt, df=2 0.4 Antag p = 0.95 og n = 1500 Begge histogrammer: Fordelingen af X2 (10,000 gentagelser) 0.2 Nederst: Ei udregnet vha. estimeret p. Effektivt sig.niv. med df=2 er 1,4% (ikke 5%) 15 10 5 Histogram for X^2, p ukendt, df=1 X2 0.0 0.4 0.8 0 0 5 15 10 116 Test af uafhængighed: Frihedsgrader Ved test af uafhængighed er der k = rc kategorier. Beregningen af de forventede værdier involvere estimation af (c-1) + (r-1) marginale sandsynligheder. Antal frihedsgrader er derfor k – s – 1 = rc – (r – 1) – (c – 1) – 1 = (r-1)(c-1). 117 Repetition: Normalfordelingen Normalfordelingen Karakteriseret ved middelværdi og standardafvigelse . 2 Notation: y ~ N(, ) betyder at y er kontinuert stokastisk variabel, der er normalfordelt med middelværdi og varians 2. Tæthedsfunktionen for normalfordelingen er f ( y) 1 2 2 y 2 exp 2 2 Egenskaber: Symmetrisk omkring f(y) > 0 for alle y. 95% 1.96 + 1.96 Goodness-of-fit for normalfordelingen Setup: X1,…,Xn stikprøve fra normalfordeling m. ukendt og 2. 2 2 ˆ ˆ x n er estimater af og ˆ Antag x i n og i 2. Vi kan inddele R i k intervaller så pi = 1/k. Her gælder c2k-3 ≤sd X2 ≤sd c2k-1 119 Goodness-of-fit for normalfordelingen Monte Carlo simulation af x2-teststørrelsen med k=6: Bemærk at den empiriske fordeling ligger mellem to c2fordelinger med hhv. 3 og 5 frihedsgrader. 120 Parametriske test c2testene vi har set på indtil nu er eksempler på såkaldt ikke-parametrisk test. Testene er ikke-parametriske, da de groft sagt ikke omhandler parametre, men kun egenskaber (enten uafhængighed eller goodness-of-fit). I det følgende vil vi betragte et parametrisk problem med tilhørende parametriske test. 121 Parametriske hypoteser Videnskabelige hypotese Vi påstår at havenisser i gennemsnit tjener mere end 42 kroner i timen. Data Vi har spurgt 36 havenisser om deres timeløn. Gennemsnittet var x = 43,2. Antagelser Lønningerne er normalfordelte og uafhængige Standardafvigelsen er kendt, og = 3,2 betegner (den sande ukendte) populationsmiddelværdi. 122 Kunstige havenisser 0 2 4 6 8 Frequency 12 Skab et ”kunstigt” normalfordelt data-sæt der passer med x = 43,2 og s=3.6. > x = rnorm(36) > x = (x-mean(x))*3.6/sd(x)+43.2 > mean(x) Histogram of x [1] 43.2 > sd(x) [1] 3.6 hist(x) Ser jo meget 35 40 45 50 normalfordelt ud…? x 123 Hypoteser og Teststørrelse Nul-hypotesen H0: ≤ 42 (alt. notation H0: = 42 ) Alternativ-hypotesen H1: > 42 Teststørrelse Stikprøvegennemsnittet: x Al tvivl skal komme H0 til gode, så vi antager at = 42. Dvs. under H0 har vi: 2 3, 2 x ~ N 42 , 36 124 Afgørelse vha. P-værdi Hypoteser: H0: ≤ 42 2 3, 2 Dvs. under H0 har vi: x ~ N 42 , 36 Jo større x , jo mere kritisk for H0. Dvs. P-værdi = P( x ≥ 43,2) = 0.0122. (Unuanceret) konklusion: Vi afviser H0 da P-værdi < 0.05 vs H1: > 42 P-værdi 42 43,2 P-værdi i R: > 1-pnorm(43.2,mean=42,sd=3.2/sqrt(36)) [1] 0.01222447 125 Afgørelse vha. kritisk værdi Beslutningsregel: Vi afviser H0, hvis x er større end den kritiske værdi. Vi ønsker et signifikansniveau på 5%. Kritisk værdi finder vi vha. R > qnorm(0.95,mean=42,sd=3.2/6) [1] 42.87726 Konklusion Da x = 43,2 > 42,9 afviser vi H0. x a = 5% 42 42,9 126 Valg af alternativ-hypotese Videnskabelige hypotese Vi påstår at havenisser i gennemsnit tjener mere end 42 kroner i timen. Hvorfor vælger vi ikke? H0: ≥ 42 vs H1: < 42 P-værdi P( x < 43,2 ) = 0.988 Konklusion: Vi kan ikke afvise H0… P-værdi 42 43,2 127 Pointer Stærkere konklusion: Vi afviser H0 Svag konklusion: Vi kan ikke kunne afvise H0. Det styrede valget af alternativ-hypotese før. Vi har lystigt regnet på ting og sager under antagelse af, at H0 er sand. Hvorfor ikke H1? Det er typisk ikke muligt, da H1 er mindre ”skarp”… HUSK: Vi skal kontrollere at antagelserne er opfyldt Her: Er stikprøven normalfordelt. 128 En tilbagevendende svaghed H0-fordelinegn Antag sandheden er = 42.1 H0: = 42 vs H1: ≠ 42 1 . 96 n Afviser H0 I virkeligheden vil gennemsnit falde her: Vi vil med mere end 95% ssh 3 .2 afvise H0 hvis 2 1 . 96 0 .1 n Dvs. hvis 2 3 .2 n 2 1 . 96 16000 0 .1 3 .2 42 Afviser H0 1 . 96 3 .2 n 42 42.1 Afviser H0 129 Pointer Med nok data ender man næsten sikkert med at afvise H0. Afviser vi H0 er vi ret sikre på, at H0 er forkert. Men det er ikke det samme som, at data afvigelse fra H0 har nogen praktisk betydning. 130 Stikprøve-gennemsnittet Lad de stokastiske variable X1, X2,…,Xn være en tilfældig stikprøve fra en population m. middelværdi og varians 2. Stikprøve-gennemsnittet af disse SV er X 1 n n Xi i 1 x er et punktestimat for . Den forventede værdi og varians for stikprøvegennemsnittet er 2 V X og E X n Ubiased / central / middelret estimator Konsistent estimator Konfidensinterval for middelværdien - Opvarmning 2 Da X ~ N ( , n ) gælder følgende: P 1 . 96 X + 1 . 96 0 . 95 n n Dvs. med 95% sandsynlighed ligger (den stokastiske variabel) X i det faste interval 1 . 96 n . Det kan omskrives til P X 1 . 96 X + 1 . 96 0 . 95 n n Dvs. det stokastiske interval X 1 . 96 95% sandsynlighed det faste tal . n indeholder med Konfidens-interval for middelværdi 0,95 0.4 0,025 0.0 0.1 0.2 0.3 0,025 -3 2.5% falder nedenfor intervallet -2 -1 0 1 2 3 x x x x x x 95% falder indenfor intervallet 2.5% falder over intervallet x Approksimativt 95% af stikprøve middelværdierne kan forventes at falde indenfor intervallet 1.96 , + 1.96 n n Omvendt, cirka 2.5% kan forventes at være under 1.9 6 n og 2.5% kan forventes at være over + 1.9 6 n . Så 5% kan forventes at være udenfor intervallet. . Konfidens-interval for middelværdi 0,95 0.4 0,025 0.0 0.1 0.2 0.3 0,025 -3 -2 -1 0 1 2 3 x * Approksimativt 95% af intervallerne omkring stikprøvex 1 . 96 n middelværdien kan forventes at indeholde den faktiske værdi af populations middelværdien, . x 1 . 96 x x x x + 1 . 96 *5% af sådanne intervaller omkring x x x 95% falder indenfor intervallet * x stikprøve middelværdien kan forventes ikke at inkludere den faktiske værdi af populations middelværdien. Konfidensinterval for middelværdien - når X er normal-fordelt eller stikprøven er stor Vi har altså P X 1 . 96 X + 1 . 96 0 . 95 n n Hvis vi erstatter estimatoren X (”et tilfældigt tal”) med estimatet x (”et fast tal”) får vi konfidensintervallet: For en stikprøve der enten er stor eller fra en normalpopulation er et 95% konfidensinterval for middelværdien når variansen er kendt x 1 . 96 Bemærk at estimatoren X er n er ersattet med estimatet x. Hypotesetest og konfidensinterval Hypoteser H0: 0 vs H1: ≠ 0 Teststørrelse x Kritiske værdier: 0 ± 1.96 /√n Konsekvens: x 0 + 1 . 96 Hvis 0 1 . 96 n n ß x 1 . 96 0 x + 1 . 96 Hvis n n Ej afvis H0 Ej afvis H0 136 Konfidensinterval for andele Antal x~binom(n,p) Dvs. x er antal succeser i n uafhængige forsøg. pˆ x er et estimat af p. n ca pˆ ~ N ( p , p (1 p ) n ) Der gælder (jf. CLT) Et tilnærmet 95% konfindensinterval for p: pˆ 1, 96 pˆ (1 pˆ ) n 137 Monte Carlo p-værdi Simuler data under H0 og udregn teststørrelse n gange. Lad x være antal teststørrelser der er mere ekstreme end den observerede. pˆ x n er et estimat af p-værdien. Et tilnærmet 95% konfindensinterval for p-værdien: pˆ 1, 96 Afvigelsen1,96 pˆ (1 pˆ ) n pˆ (1 pˆ ) n kaldes også Monte Carlo fejlen. 138 Monte Carlo stikprøvestørrelse Konfidensintervallet er pˆ D hvor fejlmarginen er Isoler n giver D 1, 96 pˆ (1 pˆ ) n 1, 96 ˆ ˆ n p (1 p ) D 2 Hvis D = 0.001 og p = 0.05 fås 2 1, 96 n 0 . 05 0 . 95 182476 0 , 001 En fejlmargin på 1 promille kræver næsten 200.000 simulationer. 139 Students t fordeling Antag populationen er normalfordelt med middelværdi og varians 2. Gammel viden: Hvis vi kender variansen 2, så kan vi bruge: X ~ N 0 ,1 n Ny viden: Hvis vi ikke kender variansen 2, så kan vi erstatte 2 med stikprøve-variansen s2: X s n ~ t n 1 ”følger en t-fordeling med n-1 frihedsgrader”. Students t fordeling: definition Lad Z~N(0,1) - standard normalfordeling X2 ~ c2(k) - c2-fordeling med k frihedsgrader. Z og X2 er uafhængige Lad T Z X 2 k Så følger T en t-fordeling med k frihedsgrader. Students t fordeling: egenskaber t fordelingen er klokkeformet Standard normal og symmetrisk og defineret ved antal frihedsgrader (df). t, df=20 Middelværdien er altid lig 0. t, df=10 Variansen af t er større end 1, 0 men går mod 1, når antallet af frihedsgrader vokser. t fordelingen er fladere og har ”tykkere haler” end en standard normal fordelingen. t fordelingen går mod standard normal fordelingen nå antallet af frihedsgrader vokser. Konfidensinterval for når er ukendt t-fordelingen Defintion: Et (1-a)100% konfidensinterval for når er ukendt (og man antager en normalfordelt population): x ta 2 s n hvor t a 2 er værdien i t-fordelingen med n-1 frihedsgraders, hvor sandsynligheden for at t er højere end denne værdi, er a. a/2 ta/2 t-test Hypoteser H0: 0 vs Teststørrelse t x 0 s n H1: ≠ 0 x 0 se se = s/√n er standardfejlen. Kritiske værdier: ± ta/2 144 t-test i R t.test(x,mu=42,alternative=”two.sided”) One Sample t-test data: x t = 2, df = 35, p-value = 0.05331 alternative hypothesis: true mean is not equal to 42 95 percent confidence interval: 41.98194 44.41806 sample estimates: mean of x 43.2 145 Parametriske vs Ikke-parametriske metoder Parametriske metoder Udgangspunktet er en statistisk model. Hypoteser omhandler parametrene i modellen. Undersøgelsen bygger på modelantagelser og er ”en præcis løsning, til et approksimeret problem”. Ikke-parametriske metoder Undersøgelsen bygger ikke på antagelser om specielle fordelinger og er ”en approksimeret løsning til et præcist problem”. Har lavere styrke en parametriske metoder. 146 Ikke normalfordelt stikprøve Udgangspunkt: Stikprøve x1,..,xn fra én (symmetrisk) population. Alternativ til t-test af H0: 0. Lad m betegne medianen i populationen. Hypoteser H0: 0 HA: ≠0 Antagelser: Fordeling er symmetrisk. Observationerne er indbyrdes uafhængige Observationerne er skala variable 147 Ranks Mange ikke-parametriske test benytter sig af ranks. Mindste tal tildeles rank 1, næstmindste rank 2 osv. Hvis flere observationer tager samme værdi tildeles de et gennemsnitsrank. Eksempel: Antag vi har følgende observ.: (fodtegn er rank) 31 42 7? 7? 105 146 157 De to 7-taller ”burde” have fået rank 3 og 4 Derfor får de rank (3+4)/2 = 3.5 Resultat 31 42 73.5 73.5 105 146 157 148 Wilcoxon Signed-rank Test Beregning: 1. Udregn afvigelser fra m0: di = xi - m0. 2. Find ranks for |di| (den absolutte/numeriske værdi af di). Dvs. mindste afvigelse får rank 1, næstmindste afvigelse får rank 2 osv. 3. Find sum af ranks af |di| hvor di > 0 - betegnes T + sum af ranks af |di| hvor di < 0 - betegnes T - Wilcoxon signed-rank teststørrelse: T=T + Store og små værdier af T er kritiske. 149 Wilcoxon Signed-rank: Store stikprøver Når n er større end 10 kan man anvende en normalfordelingsapproksimation. Under H0 gælder der for teststørrelsen T: T n ( n + 1) T 4 24 Vi kan nu standardisere z n ( n + 1)( 2 n + 1) T T T Under H0 gælder z ~N(0,1), dvs. z følger en standard normalfordeling. 150 Eksempel Data: Modtagne opkald i løbet af en time. Vi observeret data for 25 timer. Hypotese: Medianen er 149 H0: m = 149 vs HA: m ≠ 149 Sum af positive ranks T+ = 163.5 xi 151 144 123 178 105 112 140 167 177 185 129 160 110 170 198 165 109 118 155 102 164 180 139 166 182 di 2 -5 -26 29 -44 -37 -9 18 28 36 -20 11 -39 21 49 16 -40 -31 6 -47 15 31 -10 17 33 |di| 2 5 26 29 44 37 9 18 28 36 20 11 39 21 49 16 40 31 6 47 15 31 10 17 33 Rank 1 2 13 15 23 20 4 10 14 19 11 6 21 12 25 8 22 16.5 3 24 7 16.5 5 9 18 di > 0 1 15.0 di < 0 2 13 0 23 20 4 10 14 19 11 6 21 12 25 8 22 16.5 3 24 7 16.5 5 9 18 163.5 161.5 151 Eksempel - fortsat H0: m = 149 vs HA: m ≠ 149 Sum af positive ranks T+ = 163.5 Mellemregninger T n n + 1 25 25 + 1 4 162 . 5 4 T n ( n + 1)( 2 n + 1) 24 33150 25 ( 25 + 1)( 2 25 + 1) 24 37 . 165 24 Teststørresle Z T T T 163 . 5 162 . 5 0 . 027 37 . 165 Konklusion Da Z 0.027 1.96 kan vi ikke afvise H0. p-værdi = 0.979. 152 Wilcoxon Signed-rank i R > wilcox.test(messages$Messages, mu=149) Wilcoxon signed rank test with continuity correction data: messages$Messages V = 163.5, p-value = 0.9893 alternative hypothesis: true location is not equal to 149 Samme resultat som før. 153 Mere avancerede test I det følgende vil vi kort beskæftige os med Simpel lineær regression Til analyse af en lineær sammenhæng mellem to variable Variansanalyse Til at undersøge om middelværdien er den samme i et antal grupper. 154 Kriminalitet og uddannelse i Florida: Er der en sammenhæng? crime = read.table("fl-crime.csv", header=TRUE,sep=";",dec=",") plot(crime$C ~ crime$HS) Plot af ”kriminalitet” (y) mod ”Andel high school” (x): 80 40 0 crime$C 120 55 60 Er der en sammenhæng? 65 70 crime$HS 75 80 85 Scatterplot Y Et scatterplot er et plot af to variable: x : forklarende variabel (xi,yi) yi (percent high school) y : respons variabel (crime rate) For den i’te observation har vi xi (crime rate for i’te distrikt) yi (% high school for i’te distrikt) Data: (x1,y1), (x2,y2),…, (xn,yn) xi x Forventet respons: En ret linje Den rette linje a + bx beskriver den forventede (dvs. middel) respons: y UK: Expected E[y] = a + bx E[y] = a + bx Eksempel: E[y] = 210 + 2,5x b Fortolkning: Antag x = 40 (% high school), 1 så er den forventede crime rate 210 + 2,5·40 310 a Hvis x øges med 1, så øges x den forventede værdi af y med 2,5. Hvis x = 0 , så er den forventede værdi af y = 210. Fejlleddet De enkelte datapunkter (xi,yi) ligger ikke præcist på regressionslinjen. yi Afvigelsen mellem punkt og linjen betegnes fejlleddet ei. y (xi,yi) a + bx ei Regressionsmodel: yi = a + bxi+ ei Bemærk: n fejlled e1, e2, ..., en. Flere detaljer og antagelser på næste slide… xi x Simpel lineær regressionsmodel y i b 0 + b 1 xi + e i •Y •X •β •β0 •β1 •iid •ε •εi e i iid N ( 0 , ) 2 - den afhængige variabel. - den uafhængige variabel – faste - det græske bogstav ”beta” - skæringspunkt med y-aksen - hældningskoefficient - UK: independent, identically distributed = uafhængig, identisk fordelte - det græske bogstav ”epsilon” - det eneste stokastiske element i modellen Lineær regressionsmodel: Figur Model: yi = a + bxi+ ei Om fejlledene ei antager vi: Normalfordelt Middelværdi nul Konstant standardafvigelse Dvs. punkterne ligger usystematisk spredt omkring en ret linje, hvor variationen er konstant. Yi b 0 + b 1 x i + e i Y Fordelingen af yi omkring regressionslinjen. i.i.d. normalfordelte fejlled X x1 x2 x3 x4 x5 Kontinuert forklarende variabel x Visuelt check af antagelser Lav et scatter plot y √ y % x x √ y x % y x En tilnærmet linje En estimeret regressionslinje er givet ved: yˆ = a + bx y Her er a et estimat af a b et estimat af b ”y hat” er estimat af E(y) Afstanden fra punktet til den estimerede regressionslinje kaldes residualet ei = yi - yˆ i . (xi,yi) E[y] = a + bx yi ei yˆ = a + bx yˆ i xi x Mindste kvadraters metode y Summen af de kvadrede residualer betegnes: n SSE y i yˆ i 2 i 1 (xi,yi) n e 2 yi i ei i 1 UK: Sum of Squared Errors. SSE kan skrives som E[y] = a + bx yˆ = a + bx yˆ i n SSE 2 y a + bx i i i 1 Vi vælger a og b, så SSE er mindst mulig. Dette kaldes mindste kvadraters metode. xi x Forklaret og uforklaret afvigelse yi’s afvigelse fra y kan opdeles i to: y yˆ a + bx y i Uforklaret afvigelse yˆ i Totale afvigelse Forklaret afvigelse y x xi x 164 Multipel determinations koefficient Den totale variation i y’erne: TSS (Total Sum of Squares) Den uforklarede del af variationen i y’erne: SSE 2 y y i i i yi yˆ i 2 2 i ei (Sum of Squared Errors) Den forklarede del af variationen i y’erne: SSR 2 ˆ y y i i (Sum of Squars for Regression) 165 Multipel determinations koefficient Der gælder TSS SSE + SSR Dvs. Forklarede var. = Uforklarede var. + Forklarede var. Determinationskoefficienten R 2 SSR TSS TSS SSE TSS Fortolkning: Andelen af den totale variation, der er forklaret. 166 Hypotesetest af b Nul-hypoteser: H0: b = 0 Alternativ-hypoteser: Ha: b 0 Ha: b > 0 Teststørrelse b t se Ha: b < 0 Hvis H0 er sand, så følger t en tfordeling med df=n-2 frihedsgrader Fortolkning af H0: β = 0 Er der en lineær sammenhæng mellem X og Y? H0: β1 = 0 Ha: β1 ≠ 0 ingen lineær sammenhæng lineær sammenhæng Følgende er eksempler, hvor H0 accepteres. Konstant Y Usystematisk variation Y Y X Ikke-lineær sammenhæng Y X X Simpel lineær regression i R model = lm(C~HS,data=crime) Summary(mode) Call: lm(formula = C ~ HS, data = crime) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -51.8018 25.1447 -2.060 0.043588 * HS 1.5010 0.3611 4.156 0.000101 *** Konklusion: Andel af folk der har gennemført high school (HS) påvirker omfanget af kriminalitet (C). Regressionslinje i SPSS 80 100 60 40 20 0 plot(crime$C~crime$HS) abline(model) crime$C 55 60 65 70 crime$HS 75 80 85 Multipel lineær regression Vi tilføjer graden af urbanisering: model2 = lm(C ~ HS + U, data = crime) summary(model2) Coefficients: Estimate Std. Error t value (Intercept) 56.8006 29.3996 1.932 HS -0.5413 0.4906 -1.103 U 0.6734 0.1276 5.279 Pr(>|t|) 0.058 . 0.274 1.82e-06 *** Umiddelbart er det kun urbaniseringen og ikke uddannelse, der har en betydning for kriminaliteten. 171 Proxy 85 75 Urbanization 65 Education 55 High school er en proxy (”erstatnings variabel”) for urbanisering: Crime rate plot(crime$HS ~ crime$U) crime$HS 0 20 40 60 80 100 crime$U 172 Variansanalyse (ANOVA Analysis of Variance ) Setup: Kun kategoriske forklarende variable Eksempel: Y: Månedlige forbrug (Amount spent - amtspend) X1: Shoppestil (Shopping style - style) Hver anden uge: Biweekly (B) Hver uge: Weekly (W) Ofte: Often (O) Spørgsmål: Påvirker ’style’ forbruget? Grafisk overblik 200 300 400 500 600 700 plot(grocery$amtspent~grocery$style) grocery$amtspent Biweekly Weekly Often grocery$style 174 Omkodning vha. Dummies For at kunne anvende en MLR model må den kategoriske style variabel omkodes til dummy variable: To binære dummy variable: XB og XW Style XB XW Biweekly 1 0 Weekly 0 1 Often 0 0 Bemærk: k kategorier omkodes til k-1 dummy variable Model: Y a + b B x B + b W xW + e Hypotesen Model: Y a + b B x B + b W xW + e E[Y | Style = B] = a + bB E[Y | Style = W] = a + bW E[Y | Style = O] = a Bemærk: bB og bW angiver hvordan Bi-weekly og Weekly adskiller sig fra Often. Often er referencekategori. Hypotese: Middelværdien er den samme for alle styles: H0: bB bW 0 H1: bB 0 og/eller bW 0 Afgøres vha. et F-test. ANOVA i R model3 = lm(amtspent~factor(style),data=shopping) summary(model3) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 378.52 11.73 32.279 <2e-16 *** styleWeekly 26.03 13.45 1.936 0.0537 . styleOften 28.25 17.34 1.629 0.1042 anova(model3) Analysis of Variance Table Response: amtspent Df Sum Sq Mean Sq F value Pr(>F) style 2 39579 19789.4 2.0558 0.1295 Residuals 348 3349883 9626.1 Bayesiansk statistik Indtil nu har vi tænkt på sandsynligheder som andelen af succeser ”i det lange løb”: #succeser / #forsøg → P(succes). I bayesian statistik er sandsynligehder subjektive! Eksempler Sandsynligheden for at to virksomheder fusionerer Sandsynligheden for at en aktiekurs stiger Sandsynligheden for at det regner i morgen Typisk vil vi udtale os om en parameter q, fx , 2 or . Hvordan gøres det med subjektive sandsynligheder? 178 Bayesiansk statistik: Prior Bayesianske ide: Be beskriver vores “viden” om parameteren q vha. tæthedsfunktion (q). Denne er kendt som a priori foirdelingen (eller bare ‘prioren’) – idet den beskriver situationen inden vi har set data. Eksempel: Antag q er sandsynligheden for succes (0q1. Prioren beskriver de værdier vi tror q har: 179 Bayesian statistik: Posterior Lad x betegne vores data. Den betingede fordeing af q givet x betegnes posterior fordelingen: q | x f x | q q g x Her betegner f(x|q) datas fordeling betinget af q. Eksempel: Lad x betegne antal succeser i n forsøg. Betinget af q, følger x en binomialfordeling: n x n x f ( x | q ) q (1 q ) x 180 Bayesiansk statistik Vi observerer n = 100 forsøg med x = 30 successer, dvs. x/n = 0.3 Posterioren – vores ”viden” efter at have set data: Gråt område: A priori fordelingen Linje: A posteriori fordelingen Bemærk at de tre a posteriori fordelinger ligner hinanden. 181 Bayesiansk statistik: Matematikken bag Som prior har vi brug en såkaldt beta-fordeling med parametre a > 0 and b > 0: (q ) (a + b ) (a ) ( b ) q a 1 (1 q ) b 1 for 0 q 1 Posterioren er da (q | x ) (a + b + n ) (a + x ) ( b + n x ) q a + x 1 (1 q ) b + n x 1 en beta fordeling med parametre a+x og b+nx. 182