Matematikcentrum Matematisk statistik Lunds universitet MASB11 HT14, lp2
Laboration 4
Statistiska test 2015-01-09
Del I: Styrkefunktion
Del II: Standardtest
2
Syftet med laborationen är att ni ska bekanta er med olika typer av hypotestest och dess egenskaper:
Ni ska träna på de grundläggande begreppen inom hypotesprövning (t.ex.
signifikansnivå och styrka) samt vilka slutsatser man kan dra från analysen (del I av handledningen).
Ni skall bekanta er med lite av de funktioner som finns i R vad det gäller olika grundläggande statistiska test, främst parametriska men i viss mån även icke- parametriska test (del II av handledningen). Till varje parametriskt test ges också motsvarande icke-parametriska alternativ.
Del I: Grundläggande begrepp vid hypotesprövning
Denna del av laborationen fokuserar på de grundläggande egenskaperna hos ett hypotestest, felrisker, signifikansnivå samt styrka. Det finns också tre videoinspelningar på kursens hemsida som presenterar begreppen.
Exempel: Muntorrhet
Läkemedel kan ge en nedsatt salivkörtelfunktion, vilket är en riskfaktor för karies och andra sjukdomar i munhålan. På 7 slumpmässigt utvalda patienter som alla fick samma medicin mätte man under 5 minuter den så kallade tuggstimulerade saliven. Normal mängd saliv under dessa förhållanden är 1 ml/min och muntorrhet anses föreligga när mängden saliv understiger 0.7 ml/min. Som modell antog man att salivmängden är normalfördelad med väntevärde µ och standardavvikelse σ, där σ anses vara 0.5 ml/min. Intressanta frågeställningar är t.ex.:
Stöder data vår misstanke att medicinen sänker salivproduktionen?
Om medicinen ger upphov till en genomsnittlig salivproduktion på 0.8 ml/min, hur troligt är det att vi kommer att missa den nedsatta salivproduktionen med vårt test?
Hur många patienter ska vi mäta på om vi vill att testet ska upptäcka en nedsatt salivproduktion på 0.7 ml/min med sannolikheten 0.95?
På kursens hemsida hittar ni data i filen Saliv.RData. Kopiera också ner rutinerna Hypotes.R och Styrkefkn.R från kursens hemsida. Kortfattade svar till frågorna som ställs i uppgifterna finns i slutet på Del I av handledningen.
1. De grundläggande begreppen med rutinen Hypotes.R
a) Först vill man undersöka om data från de 7 patienterna stöder vår misstanke att medicinen sänker salivproduktionen Ställ upp lämpliga hypoteser.
b) Beräkna medelvärdet av mätningarna.
c) Använd rutinen Hypotes.R för att illustrera testets kritiska område då testet utförs på signifikansnivå α = 0.05. Det aktuella kommandot är hypotes1(σ, n, µ0, α, H1-riktn), så om hypoteserna är H0: µ = 1; H1: µ < 1 blir kommandot hypotes1(0.5, 7, 1, 0.05,
’<’).
d) Rutinen markerar det kritiska området och anger ett värde k som är gränsen till området. Hur har k beräknats?
e) Använd ditt beräknade medelvärde för att utföra testet. Vad är din slutsats om H0? f) Vad är din konkreta tolkning av signifikansnivån α = 0.05 i detta exempel?
g) Undersök hur det kritiska området (gränsen k) ändras då du ändrar signifikansnivån till α = 0.01. Vad blir din slutsats nu?
3
2. Testets styrka och styrkefunktion med rutinerna Hypotes.R och Styrkefkn.R
a) Antag nu att genomsnittlig salivutsöndring i riskgruppen i själva verket är 0.8. Då är förstås H0: µ = 1 falsk och vi vill att vårt test ska upptäcka detta och förkasta denna hypotes till förmån för hypotesen H1: µ < 1. Sannolikheten att testet verkligen klarar av detta kallas för testets styrka i punkten 0.8. Använd rutinen Hypotes.R för att illustrera testets styrka i punkten 0.8. Kommandot är nu
hypotes1(σ, n, µ0, α, H1-riktn, sant µ), så i detta fall skriver du hypotes1(0.5, 7, 1,0.05, ’<’, 0.8).
Rutinen ger dig ytterligare en figur som, förutom signifikansnivån α (felrisk av typ I), även visar β (felrisk av typ II). Vad är den konkreta tolkningen av β i detta exempel?
Hur hänger β ihop med testets styrka?
b) Mer generellt, testets styrka i punkten µ, är P(H0 förkastas då µ är verklig
genomsnittlig salivutsöndring i riskgruppen). Observera att styrkan beror på värdet µ.
I detta exempel gäller att ju mindre µ är i förhållande till µ0 = 1 desto större är chansen att testet ska upptäcka att H0 inte gäller. Därför är det intressant att studera styrkan som en funktion av µ, denna funktion betecknas ofta S(µ), (eller h(µ)).
Rutinen Styrkefkn.R ritar upp styrkefunktionen, kommandot är styrka1(σ, n, µ0, α, H1- riktn, sant µ) om du vill rita upp funktionen och markera ett speciellt µ-värde. Använd alltså kommandot styrka1(0.5, 7, 1,0.05, ’<’, 0.8), vilket ger dig den tidigare figuren plus styrkan som en funktion av µ.
Uppskatta utifrån styrkefunktionen hur stor sannolikheten är att vi med vårt test kommer att upptäcka att en grupp som bör klassas som muntorra (µ= 0.7) har en sänkt salivproduktion. Hur många patienter bör vi mäta på om vi med sannolikheten 0.95 verkligen ska upptäcka att muntorra har en sänkt salivproduktion?
Svar:
1a) H0: µ = 1; H1: µ < 1
1d) k=µ0 - λα σ/√n = 1-1.6445 ∙0.05/√7=0.689
1e) Eftersom medelvärdet < 0.689 förkastas H0 på nivå 0.05
1f) Det är 5% risk att vi påstår att riskgruppen har en sänkt salivproduktion när den i själva verket är normal
1g) H0 kan ej förkastas på nivå 0.01
2a) β=P(ej förkasta H0 då verklig genomsnittlig salivproduktion i riskgruppen är 0.8). β=1- S(0.8), d.v.s. β är 1- styrkan i punkten 0.8
2b) Styrkan i punkten 0.7 är S(0.7) vilket enligt figuren kan uppskattas till 0.48. Det krävs n=30 patienter för att styrkan ska vara 0.95 i punkten 0.7.
4
Del II: Statistiska standardtest i R
Kortfattade svar till frågorna som ställs i uppgifterna finns i slutet på handledningen.
På kursens hemsida hittar ni de fem datamaterial som används i denna del av handledningen.
Datamaterialen som används är Albumin, Median, Klorofyll, Dammar och Gråsparvar.
1 Test av medelvärde/median från en population.
Använd datamaterialet Albumin.RData.
a) t-test. En blandning av blodserum innehåller exakt 42 g albumin per liter. Två laboratorier (A och B) får göra sex bestämningar var av koncentrationen. Vi vill undersöka om det finns någon systematisk avvikelse från det sanna värdet (42 g/l) (tvåsidig mothypotes). I R gör man två t-test med nedanstående kommandon i kommandofönstret (glöm inte att spara det som, t.ex. lab3.R) och kör dem.
t.test(AlbuminA, mu=42) t.test(AlbuminB, mu=42)
Tolka utskriften. Hur stora är p-värdena och vad blir slutsatserna? Vad blir konfidensintervallen för medelkoncentrationerna?
b) Test av median. Tillämpning av Binomialfördelningen. Använd datamaterialet Median.RData. I materialet finns total längd för ett stickprov om 16 gråsparvar. Vi vill testa om populationens medianlängd är större än 160 mm. Detta kan göras genom att undersöka hur många gråsparvar i stickprovet som är längre än 160 mm. Under H0 är dessa antal Bin(n, 0.5). Detta innebär att vi kan använda ett test för p i en
binomialfördelning. Vi behöver då n (=16) och beräkna x = antal som är längre än 160 mm och testa H0: p=0.5:
antal <- length(Median) stora <- sum(Median > 160) binom.test(x=stora, n=antal, p=0.5)
Lägg märke till att det p-värde vi får gäller ett tvåsidigt test. Om vi vill ha ett ensidigt test får vi använda ett av kommandona
binom.test(x=stora, n=antal, p=0.5, alternative=”less”) binom.test(x=stora, n=antal, p=0.5, alternative=”greater”)
Om vi vill testa om medianen är större än 160, vilket alternativ ska vi då använda?
Vad blir resultatet om vi vill ha ett ensidigt test? Hur stort är p-värdet och vad blir slutsatsen?
c) Rang-test. Alternativ variant av mediantestet är att använda ranger istället. När vi bara undersöker hur stor andel av värdena som ligger över resp. under den önskade medianen förlorar vi informationen om hur långt från medianen de ligger. Om de värden som ligger över 160 dessutom ligger längre från 160 än vad de värden som ligger under 160 gör tyder det ju på att medianen borde vara större än 160. Detta utnyttjar man i ett Wilcoxon signed rank-test:
5 wilcox.test(Median, mu=160)
wilcox.test(Median, mu=160, alternative=”less”) wilcox.test(Median, mu=160, alternative=”greater”)
Vad blir slutsatsen nu? Om vi vill testa om medianen är större än 160, vilket ensidigt alternativ ska vi välja? (Ignorera varningarna om cannot compute exact test with ties.) 2 Test vid jämförelse av två populationer.
Använd datamaterialet Klorofyll.RData.
a) t-test (oberoende grupper). Alger fick växa under ljusa respektive mörka
förhållanden. Undersök med ett t-test om det finns skillnader i förväntad klorofyllhalt mellan de två grupperna genom att använda t.test igen men nu genom att lägga till gruppvariabeln (jämför Lab1 där vi delade upp överlevnaden på hannar och honor med ~):
t.test(Klorofyll$Klorofyll ~ Klorofyll$Grupp)
Tolka utskriften. Om vi inte säger något annat förutsätter t.test att varianserna i de två grupperna är olika och kompenserar för det. Om vi vet (eller antar) att varianserna är lika kan vi utnyttja det:
t.test(Klorofyll$Klorofyll ~ Klorofyll$Grupp, var.equal=TRUE)
Blev det någon skillnad? För säkerhets skull är det bäst att verkligen testa om varianserna är olika:
var.test(Klorofyll$Klorofyll ~ Klorofyll$Grupp)
Är varianserna olika? Vilket av de två t-testen ska vi använda? Slutsatser?
Vi kan få en grafisk beskrivning av skillnaderna med en boxplot:
plot(Klorofyll$Klorofyll ~ Klorofyll$Grupp) Verkar det stämma med slutsatsen ovan?
b) Rangsummetest. Gör nu om testet i 2(a) som ett rangsummetest genom wilcox.test(Klorofyll$Klorofyll ~ Klorofyll$Grupp)
Hur stort är p-värdet och vad blir slutsatsen? Jämför resultatet med t-testet som ni gjorde i (a)-uppgiften.
3 Test vid matchade data.
Använd datamaterialet Dammar.RData.
a) t-test. Jämför om det finns någon skillnad i kvävebelastning mellan vår och sommar. Vi måste nu tala om för t.test att vi har två variabler som är parade:
t.test(Dammar$N_belast_V, Dammar$N_belast_S, paired=TRUE)
6
Finns det några skillnader mellan vår och sommar?
b) Teckentest. Gör nu motsvarande analys som ett icke-parametriskt test med ett teckentest som räknar hur ofta vår-värdet överstiger motsvarade sommar-värde:
antal <- nrow(Dammar)
stora <- sum(Dammar$N_belast_V > Dammar$N_belast_S) binom.test(x=stora, n=antal, p=0.5)
Finns det några skillnader? Jämför resultatet med t-testet som ni gjorde i (a)-uppgiften.
c) Rangsummetest. Utnyttja nu storleken på värdena med ett rangsummetest istället:
wilcox.test(Dammar$N_belast_V, Dammar$N_belast_S, paired=TRUE) Vad blir slutsatsen nu?
4 Test av kategoridata.
Använd datamaterialet Gråsparvar.RData.
a) χ2-test. Vi vill nu testa om vi har en jämn könsfördelning i materialet. Om man gör en enkel tabell på variabeln kön får man följande resultat:
antal <- table(Grasparvar$Kön) andel <- prop.table(antal) cbind(antal, andel) antal andel Hane 87 0.6397059 Hona 49 0.3602941
Vi kan utföra testet på flera olika sätt. För det första kan vi göra et χ2-test där vi jämför den observerade fördelningen med den förväntade under hypotesen om jämn
könsfördelning. Detta görs med chisq.test(antal)
Hur stort är p-värdet och vad blir slutsatsen?
Om man inte säger något annat testar R om andelarna är lika. Vill man t.ex. testa om det kan vara 60 % hanar och 40 % honor får man göra så här:
chisq.test(antal, p=c(0.6,0.4)) Vad blir slutsatsen nu?
b) Test av ett proportionstal. Vi kan också göra samma test med hjälp av binomialfördelningen:
antal <- nrow(Grasparvar)
hannar <- sum(Grasparvar$Kön=="Hane") binom.test(x=hannar, n=antal, p=0.5)
7
Gör testet. Ändrar sig slutsatsen? En fördel här är att vi på köpet fick ett
konfidensintervall för andelen. Men å andra sidan fungerar det bara när vi har två kategorier. Har vi fler måste vi använda χ2-testet.
c) Test av två proportionstal. Testa nu om andelen som överlever är lika stor bland honor som bland hanar. Detta görs med med hjälp av ett χ2-test på en korstabell:
tabell <- table(Grasparvar$Kön, Grasparvar$Överlevnad) prop.table(tabell, margin=1)
chisq.test(tabell)
Hur stort är p-värdet och vad blir slutsatsen?
Svar:
1. a) A: p=0.033; B: p=0.081; Vi kan påvisa en skillnad för A men inte för B.
95% KI – A: (42,06, 42.94); B: (35.69, 42.51) b) p=0.038; Vi kan påvisa att Md>160
c) p=0.00024; Vi kan påvisa att Md>160.
2. a) Varianserna är inte olika (p=0.48); t-test: p=0.009; Vi kan påvisa en skillnad i klorofyll.
b) p=0.015; Vi kan påvisa en skillnad. Samma slutsats men inte lika signifikant.
3. a) p=0.262 Vi kan inte påvisa någon skillnad b) p=0.070; Vi kan inte påvisa någon skillnad c) p=0.195; Vi kan inte påvisa någon skillnad.
4. a) p=0.0011; Vi kan påvisa en avvikelse från jämn könsfördelning; p=0.345: könsfördelningen skulle kunna vara 60/40.
b) p=0.0014: Vi kan påvisa en avvikelse från jämn könsfördelning c) p=0.112: Vi kan inte påvisa någon skillnad i andelen som dör
