library(openxlsx)
library(mosaic)
library(vcd)
library(dplyr)
Varning
Den hÀr labben förutsÀtter att följande paket finns installerade:
openxlsxmosaicvcddplyr
Paket kan installeras via kommandot install.packages('packagename'), dÀr 'packagename' Àr namnet pÄ paketet, t.ex 'mosaic'. DÀrefter behöver paketen ocksÄ aktiveras med funktionen library()
Introduktion
I den hÀr datorlabben kommer vi att beskriva samband mellan tvÄ kategoriska variabler, samt samband mellan en numerisk och kategorisk variabel. Ni kommer ocksÄ lÀra er om betingade fördelningar i R och hur man jÀmför fördelningar grafiskt. Vi kommer ocksÄ att arbeta med tidsserier.
đȘ Uppgift 0.1
Se till att paketen ovan Àr installerade och aktiverade innan du fortsÀtter med resten.
đȘ Uppgift 0.2
Skapa en mapp Lab3 i din kursmapp SDA1 (sĂ„som ni gjorde i Lab 1). Ladda ner Quarto-filen för denna lab genom att högerklicka hĂ€r och vĂ€lj âSpara lĂ€nkâ eller nĂ„got likande frĂ„n menyn som dyker upp. Spara filen i din nya Lab3 mapp. Ăppna Quarto-filen i RStudio. Du kan nu fortsĂ€tta med denna laboration direkt i Quarto-dokumentet, dĂ€r du ocksĂ„ fyller i svaren pĂ„ dina laborationsövningar. Du kan alltsĂ„ lĂ€mna den hĂ€r webbsidan nu och fortsĂ€tta arbetet i RStudio.
VĂ€l inne i RStudio med öppnat Quarto-dokument i âEditorâ kan ni gĂ„ över till âSource modeâ genom att klicka pĂ„ âSourceâ i det vĂ€nstra hörnet av din âEditorâ. Source mode Ă€r detaljerat och bra att skriva kod i eftersom man har full kontroll pĂ„ dokumentet, men det Ă€r svĂ„rt att fĂ„ en översikt av dokumentet. Prova nu att gĂ„ över till âVisual modeâ genom att klicka pĂ„ âVisualâ i det vĂ€nstra hörnet av din âEditorâ. Vi rekommenderar att ni mestadels arbetar i Visual mode, men att gĂ„ över till Source mode nĂ€r man verkligen vill fĂ„ till nĂ„gon detalj som Ă€r svĂ„r att Ă€ndra i Visual mode. ÂŽ
đȘ Uppgift 0.3
Klicka pÄ knappen Render i Editor-fönstret för att kompilera filen till en webbsida (html). Webbsidan kommer antingen att visas i Viewer-fönstret i RStudio eller i webblÀsaren pÄ din dator. Om din webbsida visas i webblÀsaren rekommenderar vi att du Àndrar instÀllningarna i RStudio sÄ webbsidan visas i Viewer-fönstret. Du stÀller in detta pÄ menyn Tools/Global Options⊠och sen vÀljer du Viewer Pane vid Show output preview in:
đȘ Uppgift 0.4
Quarto sÀkerstÀller att man befinner sig i korrekt arbetsmapp nÀr koden i dokumentet exekveras. Med korrekt arbetsmapp menas den mappen ni sparade .qmd filen i. Vill du komma Ät dataseten via load() kommandot mÄste dataseten vara sparade i samma mapp.
Ett vanligt arbetssĂ€tt Ă€r att man jobbar i RStudio i en separat .R fil, dĂ€r man kan testa att köra kod innan den kopieras över till Quarto dokumentet. Den .R filen kommer inte att ha samma âworking directoryâ som Quarto filen. Du mĂ„ste dĂ„ anvĂ€nda setwd() funktionen i .R filen för att stĂ€lla in âworking directoryâ. Notera att man ocksĂ„ vĂ€lja att skriva kod i Console som finns lĂ€ngst ner i RStudio. DĂ€r mĂ„ste du ocksĂ„ stĂ€lla in rĂ€tt âworking directoryâ. Det Ă€r inte rekommenderat att anvĂ€nda Console eftersom koden inte sparas dĂ€r. Du kanske kommer pĂ„ nĂ„got jĂ€ttesmart som du glömmer att kopiera över till Quarto dokumentet och kan senare inte hitta det du skrev.
RStudios Editor Àndras beroende pÄ vilken sorts fil du har öppen
Notera att ikonerna i Editor Àr annorlunda nÀr du har ett Quarto-dokument öppet jÀmfört med tidigare nÀr vi hade en fil med ren R-kod (dvs en .R fil) öppen.
Skapa en ny .R som du döper till Lab3_test_code.R och stĂ€ll in âworking directoryâ till den nya mappen genom att följa anvisningarna frĂ„n Lab 1.
đȘ Uppgift 0.5
Ladda ner SmartPhones.RData (hÀr), FevChildren.RData (hÀr) och CAPM_data.RData (hÀr). Filerna kommer automatiskt att sparas ner, oftast i download mappen pÄ datorn ni sitter vid. Kopiera över filerna till Lab3 mappen ni skapade ovan.
Det finns olika sÀtt att lÀsa in .RData filer genom load() funktionen. Ett sÀtt Àr att skriva t.ex load("SmartPhones.RData"). Notera att om du gör detta i Lab3_test_code.R förutsÀtter det att SmartPhones.RData ligger i den arbetsmappen man angivit genom setwd() funktionen. Ett annat sÀtt att ladda in en .RData fil Àr att lÀsa in den direkt frÄn webben med en kombination av load() funktionen och url() funktionen.
load(file = url("https://github.com/StatisticsSU/SDA1/blob/main/datorlab/lab3/SmartPhones.RData?raw=true"))
Smartphones_from_URL <- SmartPhonesArgumentet till url() funktionen Àr en strÀng som innehÄller lÀnken till kurshemsidans git repository dÀr SmartPhones.RData finns lagrad. I raden efter har vi sparat dataframen i en variabel som heter SmartPhones_from_URL sÄ att du i nÀsta uppgift kan jÀmföra mot den du lÀser in frÄn din lokala fil som du sparade ovan i i Lab3 mappen.
đȘ Uppgift 0.6
AnvÀnd load() funktionen (utan url() funktionen) för att lÀsa in SmartPhones.RData lokalt frÄn din dator i Quarto dokumentet. JÀmför den lokalt inlÀsta filen med den du lÀste in frÄn webben (finns sparad i SmartPhones_from_URL).
# Write your code here
Tips
För att jÀmföra dataseten hÀr rÀcker det att anvÀnda str() och head().
Det Àr viktigt att ni vet hur man lÀser in .RData filer lokalt. Det kan vara sÄ att github Àr tillfÀlligt nere, eller att ni temporÀrt jobbar utan tillgÄng till internet. I resten av labben sÄ lÀses filerna in frÄn webben, men nu vet ni hur ni kan göra om webben av nÄgon anledning skulle vara nere.
1. Simultan och marginella fördelningar utifrÄn en korstabell
Kategoriska variabler Àr variabler vars utfall Àr kategorier. I Lab 1 lÀrde vi oss att uttrycka fördelningen för en kategorisk variabel via en frekvenstabell uttryckt i andelar. En fördelning av en endaste variabel, till exempel \(y\), kallas för marginalfördelningen för \(y\). NÀr vi betraktar tvÄ kategoriska variabler, till exempel \(x\) och \(y\) , sÄ finns det olika fördelningar vi kan tÀnka oss:
- Marginella fördelningar, en fördelning för vardera variabel \(x\) och \(y\).
- Simultan fördelning för variablerna \(x\) och \(y\) (endast en fördelning).
- Betingade fördelningar, en fördelning för vardera variabel \(x\) och \(y\), betingat pÄ den andra.
Notera att dessa tre typer av fördelningar finns oavsett om variablerna Àr numeriska eller kategoriska. NÀr bÀgge variablerna Àr kategoriska rÀknas dessa fördelningar utifrÄn en sÄ kallad korstabell (se nedan). Att kÀnna till de olika fördelningstyperna Àr mycket viktigt och den kunskapen tillÀmpas pÄ alla statistikkurser framöver. Att förstÄ dessa, samt kunna urskilja dom frÄn varandra, kommer att underlÀtta förstÄelsen av statistik sÄvÀl som sannolikhetsteori som ni kommer stöta pÄ under Del 2 av kursen.
För att illustrera dessa koncept i ett praktiskt problem ska vi undersöka om det föreligger ett samband mellan vilken smartphone man föredrar och vilken Älderskategori man tillhör. Data finns i SmartPhones.RData.
load(file = url("https://github.com/StatisticsSU/SDA1/blob/main/datorlab/lab3/SmartPhones.RData?raw=true"))
str(SmartPhones)'data.frame': 2498 obs. of 2 variables:
$ age.group: chr "50+" "50+" "50+" "35-49" ...
$ brand : chr "Samsung" "Samsung" "Others" "Samsung" ...head(SmartPhones) # Prints first few observations age.group brand
1 50+ Samsung
2 50+ Samsung
3 50+ Others
4 35-49 Samsung
5 50+ Motorola
6 35-49 iPhoneLÄt oss börja med att skapa en korstabell. En korstabell Àr en tabell som innehÄller frekvensen för varje utfall av variablerna. Notera att om vi har tvÄ kategoriska variabler, sÀg \(x\) och \(y\), sÄ Àr utfallen alla parvisa kombinationer av utfallen av de respektive variablerna.
tally(~ brand + age.group, data = SmartPhones) age.group
brand 18-34 35-49 50+
Asus 25 17 8
Casio 25 11 6
HTC 26 11 10
iPhone 354 219 284
LG 87 83 166
Motorola 33 30 71
Others 6 6 44
Samsung 264 252 351
Sony 48 33 28Eftersom age.group har 3 utfall och brand har 9 utfall, finns det totalt \(3\cdot9=27\) olika utfall i korstabellen.
Formula-syntaxen (nÀmndes i Lab 1) ~ brand + age.group talar om att vi vill göra en frekvenstabell av brand och age.group. Vi kan ocksÄ anvÀnda tecknet & istÀllet för +.
đȘ Uppgift 1.1
Gör en korstabell genom att anvÀnda & istÀllet för +. SlÀng Àven om ordningen pÄ variablerna, dvs age.group före brand. FörÀndras frekvensantalen?
# Write your code hereVi kan lÀgga till radsummor och kolumnsummor till korstabellen genom argumentet margins = TRUE.
tally(~ brand + age.group, data = SmartPhones, margins = TRUE) # Show row and column totals age.group
brand 18-34 35-49 50+ Total
Asus 25 17 8 50
Casio 25 11 6 42
HTC 26 11 10 47
iPhone 354 219 284 857
LG 87 83 166 336
Motorola 33 30 71 134
Others 6 6 44 56
Samsung 264 252 351 867
Sony 48 33 28 109
Total 868 662 968 2498đȘ Uppgift 1.2
Svara pÄ följande frÄgor:
Hur mÄnga deltog i undersökningen?
Hur mÄnga föredrog en Casio smartphone?
Hur mÄnga i Äldersgruppen 18-34 föredrog en Samsung telefon?
Tabeller Àr ofta enklare att utlÀsa nÀr de anges i relativa frekvenser (alternativt i procentform). Detta kan Ästadkommas med argumentet format.
tally(~ brand + age.group, data = SmartPhones, margins = TRUE, format = "proportion") age.group
brand 18-34 35-49 50+ Total
Asus 0.010008006 0.006805444 0.003202562 0.020016013
Casio 0.010008006 0.004403523 0.002401922 0.016813451
HTC 0.010408327 0.004403523 0.004003203 0.018815052
iPhone 0.141713371 0.087670136 0.113690953 0.343074460
LG 0.034827862 0.033226581 0.066453163 0.134507606
Motorola 0.013210568 0.012009608 0.028422738 0.053642914
Others 0.002401922 0.002401922 0.017614091 0.022417934
Samsung 0.105684548 0.100880705 0.140512410 0.347077662
Sony 0.019215372 0.013210568 0.011208967 0.043634908
Total 0.347477982 0.265012010 0.387510008 1.000000000tally(~ brand + age.group, data = SmartPhones, margins = TRUE, format = "percent") age.group
brand 18-34 35-49 50+ Total
Asus 1.0008006 0.6805444 0.3202562 2.0016013
Casio 1.0008006 0.4403523 0.2401922 1.6813451
HTC 1.0408327 0.4403523 0.4003203 1.8815052
iPhone 14.1713371 8.7670136 11.3690953 34.3074460
LG 3.4827862 3.3226581 6.6453163 13.4507606
Motorola 1.3210568 1.2009608 2.8422738 5.3642914
Others 0.2401922 0.2401922 1.7614091 2.2417934
Samsung 10.5684548 10.0880705 14.0512410 34.7077662
Sony 1.9215372 1.3210568 1.1208967 4.3634908
Total 34.7477982 26.5012010 38.7510008 100.0000000Fördelningen ovan i de celler som inte inkluderar Total Àr den simultana fördelningen för variablerna, dvs:
tally(~ brand + age.group, data = SmartPhones, margins = FALSE, format = "proportion") age.group
brand 18-34 35-49 50+
Asus 0.010008006 0.006805444 0.003202562
Casio 0.010008006 0.004403523 0.002401922
HTC 0.010408327 0.004403523 0.004003203
iPhone 0.141713371 0.087670136 0.113690953
LG 0.034827862 0.033226581 0.066453163
Motorola 0.013210568 0.012009608 0.028422738
Others 0.002401922 0.002401922 0.017614091
Samsung 0.105684548 0.100880705 0.140512410
Sony 0.019215372 0.013210568 0.011208967En cell i den simultana fördelningen anger andelen av det totala antalet som hamnar i det utfallet. Exempelvis Àr \(\approx 0.088\) (\(\approx 8.8\%\) ) av alla 2498 deltagare i Äldersgruppen 35-49 och föredrar en iPhone. En fördelning summerar alltid till 1, vilket bekrÀftades ovan nÀr vi inkluderade Total via margins = TRUE (cellen lÀngst ner till höger).
đȘ Uppgift 1.3
Svara pÄ följande frÄgor:
Hur stor andel av alla föredrog en Casio smartphone?
Hur stor andel av alla tillhörde Äldersgruppen 18-34 samt föredrog en Samsung telefon?
Vad Ă€r problemet med pĂ„stĂ„endet âCa 14% av de tillfrĂ„gade i Ă„ldersgruppen 18-34 föredrog en iPhoneâ?
đȘ Uppgift 1.4
AnvÀnd tally() funktionen sÄsom ni gjorde i Lab 1 för att berÀkna marginalfördelningen för vardera variabel. AnvÀnd argumentet format = "proportion". KÀnner ni igen dessa fördelningar frÄn tidigare i labben?
# Write your code here.Marginalfördelningarna kan alltsÄ fÄs via kolumntotalen och radtotalen i simultanfördelningen via argumentet margins = TRUE.
đȘ Uppgift 1.5
Vad summerar varje marginalfördelning till? Du bör kunna svara pÄ den hÀr frÄgan utan att koda.
2. Betingade fördelningar utifrÄn en korstabell
I sista delfrĂ„gan pĂ„ Uppgift 1.3 frĂ„gades vad problemet med pĂ„stĂ„endet âCa 14% av de tillfrĂ„gade i Ă„ldersgruppen 18-34 föredrog en iPhoneâ var. PĂ„stĂ„endet hade varit korrekt om formuleringen var âCa 14% av alla tillfrĂ„gade tillhörde Ă„ldersgruppen 18-34 och föredrog en iPhoneâ. För att rĂ€kna ut den korrekta andelen i det felaktiga pĂ„stĂ„endet mĂ„ste vi enbart betrakta de tillfrĂ„gade som tillhör Ă„ldersgruppen 18-34, dvs 868 personer (enligt tabellen i Uppgift 1.1). Av dessa föredrog 354 personer iPhone (enligt tabellen i Uppgift 1.1). Den korrekta andelen berĂ€knas sĂ„ledes till \(354/868\approx 0.41\), eller ca 41%.
Terminologin Ă€r viktig att lĂ€ra sig. Vi sĂ€ger att vi har en betingad fördelning för smartphonepreferens betingat pĂ„ Ă„ldersgrupp (i det hĂ€r fallet Ă„ldersgruppen 18-34). Ofta sĂ€ger vi enbart âfördelning för smartphonepreferens betingat pĂ„ Ă„ldersgrupp (i det hĂ€r fallet Ă„ldersgruppen 18-34)â, och det Ă€r underförstĂ„tt att det Ă€r en betingad fördelning. Vi skulle ocksĂ„ kunna sĂ€ga âfördelningen för smartphonepreferens givet Ă„ldersgruppâ.
Vi skulle kunna kunna rÀkna ut den betingade fördelningen för de övriga utfallen för hand, Asus (\(25/868\)), Casio (\(25/868\)), HTC (\(26/868\)), osv. RÀkna ut för hand Àr viktigt pÄ tentan, sÄ se till att ni förstÄr hur man gör! I praktiken lÄter vi datorn göra rÀkningarna Ät oss. I R kan vi göra detta genom anvÀndning av formula-syntax tecknet |, som betyder betingat pÄ variabeln som följer efter tecknet.
tally(~ brand | age.group, data = SmartPhones, format = "proportion") age.group
brand 18-34 35-49 50+
Asus 0.028801843 0.025679758 0.008264463
Casio 0.028801843 0.016616314 0.006198347
HTC 0.029953917 0.016616314 0.010330579
iPhone 0.407834101 0.330815710 0.293388430
LG 0.100230415 0.125377644 0.171487603
Motorola 0.038018433 0.045317221 0.073347107
Others 0.006912442 0.009063444 0.045454545
Samsung 0.304147465 0.380664653 0.362603306
Sony 0.055299539 0.049848943 0.028925620Ovan ser vi tre olika betingade fördelningar, en för varje Äldersgrupp. En betingad fördelning Àr ocksÄ en fördelning: En fördelning summerar alltid till 1. Vi kan bekrÀfta detta för varje Äldersgrupp genom att anvÀnda margins = TRUE som argument.
tally(~ brand | age.group, data = SmartPhones, format = "proportion", margins = TRUE) age.group
brand 18-34 35-49 50+
Asus 0.028801843 0.025679758 0.008264463
Casio 0.028801843 0.016616314 0.006198347
HTC 0.029953917 0.016616314 0.010330579
iPhone 0.407834101 0.330815710 0.293388430
LG 0.100230415 0.125377644 0.171487603
Motorola 0.038018433 0.045317221 0.073347107
Others 0.006912442 0.009063444 0.045454545
Samsung 0.304147465 0.380664653 0.362603306
Sony 0.055299539 0.049848943 0.028925620
Total 1.000000000 1.000000000 1.000000000Notera att en Total för radsummorna inte finns i tabellen. Detta beror pÄ att R vet att vi betingar pÄ a~ age.group, och dÄ saknar radsummorna en betydelse och dÀrför skrivs de inte ut.
Ett alternativ sĂ€tt att skriva betingningen ~ brand | age.group Ă€r brand ~ age.group , utan att anvĂ€nda |-tecknet. R tolkar dĂ„ variabeln vĂ€nster om ~ som variabeln fördelningen ska berĂ€knas för, och variabeln till höger om ~ som variabeln det ska betingas pĂ„. Det hĂ€r alternativa sĂ€ttet att skriva formula-syntax Ă€r det som anvĂ€nds i mĂ„nga R funktioner vi kommer att stöta pĂ„ senare i kursen, till exempel lm() funktionen som anvĂ€nds i regression. Ăven plot() och boxplot() som vi anvĂ€nder senare i den hĂ€r labben anvĂ€nder den alternativa formula-syntaxen.
đȘ Uppgift 2.1
AnvÀnd den allternativa formula-syntaxen för att berÀkna den betingade fördelningen ovan. StÀm av att resultaten blir desamma.
# Write your code heređȘ Uppgift 2.2
Svara pÄ följande frÄgor:
Vilken Àr den populÀraste smartphonen i Äldersgruppen 35-49?
Hur stor andel av de i Äldersgruppen 18-34 föredrar en Sony?
Hur stor andel av de i Äldersgruppen 18-34 föredrar en Sony eller en Samsung?
Vad Ă€r problemet med pĂ„stĂ„endet âBland de som föredrar iPhone Ă€r det en större andel, ca 40%, som tillhör Ă„ldersgruppen 18-34 gentemot ca 29% som tillhör Ă„ldersgruppen 50+â?
đȘ Uppgift 2.3
Uttryck de betingade fördelningarna ovan i procentform istÀllet för andelar.
# Write your code here
Tips
Kom ihÄg att argumentet format bestÀmmer tabellens format.
I sista delfrĂ„gan pĂ„ Uppgift 2.2 frĂ„gades vad problemet med pĂ„stĂ„endet âBland de som föredrar iPhone Ă€r det en större andel, ca 40%, som tillhör Ă„ldersgruppen 18-34 gentemot ca 29% som tillhör Ă„ldersgruppen 50+â . Det finns tvĂ„ varningstecken hĂ€r. Det första Ă€r formuleringen âBland de som föredrar en iPhoneâŠâ, vilket betyder att vi bara betraktar de som har en iPhone â dvs vi mĂ„ste betinga pĂ„ iPhone! Andra varningstecknet Ă€r att vi uttalar oss om andelar för Ă„ldersgrupper â den betingade fördelningen vi skapade ovan Ă€r andelar över smartphonepreferens, inte över Ă„ldersgrupper! De andelar vi har angett i pĂ„stĂ„endet Ă€r sĂ„ledes felaktiga.
đȘ Uppgift 2.4
BerÀkna (med hjÀlp av R) den betingade fördelningen som hjÀlper oss att rÀtta till det felaktiga pÄstÄendet.
# Write your code here
Tips
Om du har förstÄtt innebörden av formula-syntaxen ~ brand | age.group (eller alternativt brand ~ age.group) borde denna uppgift vara enkel.
đȘ Uppgift 2.5
RÀtta det felaktiga pÄstÄendet i sista delfrÄgan pÄ Uppgift 2.2.
3. JÀmföra fördelningar för kategoriska variabler
NÀr vi jÀmför variablers fördelningar Àr det ofta via grafiska metoder. Vilken typ av plot som Àr lÀmpligast beror pÄ om variablerna vi jÀmför Àr numeriska eller kategoriska. Vi börjar med fallet dÄ bÄda variablerna Àr kategoriska, dvs samma scenario som i Avsnitt 2 ovan.
Visst skulle vi till exempel kunna kolla pÄ andelarna för smartphonepreferens betingat pÄ Äldersgrupp vi rÀknade i Avsnitt 2, dvs:
tally(~ brand | age.group, data = SmartPhones, format = "proportion") age.group
brand 18-34 35-49 50+
Asus 0.028801843 0.025679758 0.008264463
Casio 0.028801843 0.016616314 0.006198347
HTC 0.029953917 0.016616314 0.010330579
iPhone 0.407834101 0.330815710 0.293388430
LG 0.100230415 0.125377644 0.171487603
Motorola 0.038018433 0.045317221 0.073347107
Others 0.006912442 0.009063444 0.045454545
Samsung 0.304147465 0.380664653 0.362603306
Sony 0.055299539 0.049848943 0.028925620FrÄn dessa betingade fördelningar ser vi till exempel att iPhone Àr den populÀraste smartphonen i Äldersgruppen 18-34 och nÀst populÀraste i de andra tvÄ Äldersgrupperna. Vi ser ocksÄ att andelen som föredrog iPhone var ca 41% i Äldersgruppen 18-34, gentemot ca 33% bland Äldersgruppen 35-49 och ca 29% i Äldersgruppen 50+.
Informationen i tabellen ovan blir Ànnu tydligare i nÀr den presenteras med ett stapeldiagram för varje Äldersgrupp.
bargraph(~ brand | age.group, data = SmartPhones, type = "proportion")
đȘ Uppgift 3.1
Svara pÄ följande frÄgor:
Vad kan man sÀga om preferensen för Others-kategorin mellan Äldersgrupperna?
Vilken Àr den populÀraste smartphonen för respektive Äldersgrupp?
Ăr pĂ„stĂ„endet ââSumman av alla staplarna för Ă„ldersgruppen 50+ Ă€r 1ââ korrekt? Varför/varför inte?
Ăr pĂ„stĂ„endet ââSumman av iPhone staplarna över de tre Ă„ldersgrupperna Ă€r 1ââ korrekt? Varför/varför inte?
đȘ Uppgift 3.2
AnvÀnd bargraph() för att plotta de betingade fördelningarna för Äldersgrupperna betingat pÄ smartphonepreferens.
# Write your code heređȘ Uppgift 3.3
Svara pÄ följande frÄgor:
Hur ser Äldersfördelningen ut för smartphonepreferensen Others gentemot Casio?
Hur ser Äldersfördelningen ut för smartphonepreferensen HTC?
Ăr pĂ„stĂ„endet ââSumman av 50+ staplarna över de nio smartphonespreferensgrupperna Ă€r 1ââ korrekt? Varför/varför inte?
Vi kan ocksÄ göra ett staplat stapeldiagram (stacked barchart) för att kompakt visa resultaten i Uppgift 3.2 För detta anvÀnder vi barplot() funktionen i base R, dÄ mosaic-paketet inte har en sÄdan funktion (se nedan för en liknande plot, mosaicplotten). Vi behöver först skapa och spara den betingade fördelningen med hjÀlp av tally().
t <- tally(~ age.group | brand, data = SmartPhones, format = "prop")
barplot(t, col = c("aquamarine", "cornflowerblue", "darkblue"), xlab = "Brands", ylim = c(0,1.7), legend=c("18-34", "35-49", "50+"), las=2)
Det enda argumentet till funktionen ni inte har stött pÄ tidigare Àr las = 2, som ser till att namnen pÄ kategorierna hamnar vertikalt istÀllet för horisontellt (sÄ de fÄr plats i figuren).
đȘ Uppgift 3.4
Vilken av de tvÄ figurerna ovan tycker ni Àr enklast att jÀmföra proportionerna mellan grupperna i?
đȘ Uppgift 3.5
Om vi, för varje grupp vi betingar pÄ, istÀllet vill ha staplarna bredvid varandra (och inte ovanpÄ varandra) kan vi anvÀnda argumentet beside beside = TRUE. Gör en sÄdan plot.
# Write your code hereVi har gÄtt igenom mosaicplot pÄ FörelÀsning 4. Notera att namnet pÄ plotten inte har med mosaic-paketet att göra. En mosaicplot Àr en stacked barchart, men innehÄller ocksÄ information om marginalfördelningen pÄ variabeln man betingar pÄ. Vi kan anvÀnda funktionen mosaic() frÄn vcd-paketet. Viktigt att notera att variabeln vi betingar pÄ hamnar pÄ y-axeln, till skillnad frÄn x-axeln som mosaicplottarna i boken.
t <- tally(~ age.group | brand, data = SmartPhones, format = "prop")
vcd::mosaic(~ age.group | brand, data = SmartPhones)
Notera att vi anropar funktionen via vcd::mosaic() istÀllet för mosaic(), dÄ den senare refererar till mosaic-paketet och inte till funktionen mosaic() frÄn vcd-paketet som Àr den vi vill anvÀnda.
Variabeln smartphonepreferens har för mÄnga utfall (kategorier), vilket gör mosaicplotten olÀslig. TyvÀrr har inte funktionen argumentet las vi anvÀnde oss av för att rotera texten i bargraph(). För enkelhetens skull, lÄt oss definiera en ny dataframe dÀr vi endast behÄller repondenter som har svarat iPhone, Samsung, LG eller Sony och illustrera mosaicplotten med dessa nya data. HÀr Àr ett exempel dÀr tidyverse löser uppgiften pÄ ett mer elegant sÀtt Àn base R med funktionen filter().
SmartPhones_new <- filter(SmartPhones, brand == "iPhone" |brand == "Samsung"|brand == "LG" | brand == "Sony")Notera att i koden ovan sĂ„ har inte |-tecknet samma betydelse som i formula-syntax, dĂ€r den betyder betingat pĂ„. I koden ovan betyder |-tecknet âellerâ. filter anropet ovan behĂ„ller alltsĂ„ de observationer som har vĂ€rdet brand lika med âiPhoneâ eller âSamsungâ eller âLGâ eller âSonyâ. Tecknet == Ă€r en jĂ€mförelse, dvs filter() funktionen kollar om brand Ă€r lika med de vĂ€rdena vi har angett och om sĂ„ Ă€r fallet sĂ„ behĂ„lls den. Om brand till exempel skulle ha vĂ€rdet âCasioâ sĂ„ kommer den inte att behĂ„llas.
đȘ Uppgift 3.6
AnvÀnd tally() för att berÀkna den betingade fördelningen age.group givet brand samt marginalfördelningen för brand för den nya dataframen SmartPhones_new.
# Write your code heređȘ Uppgift 3.7
AnvÀnd vcd::mosaic() för att göra en mosaicplot för age.group givet brand för den nya dataframen SmartPhones_new. StÀmmer de betingade fördelningarna samt marginalfördelningen i plotten överens med de ni rÀknade fram i Uppgift 3.6?
# Write your code here.4. Samband mellan en numerisk och en kategorisk variabel
I Lab 2 presenterades histogrammet som ett verktyg att beskriva fördelningen för en numerisk variabel. Om vi vill studera hur sambandet mellan en numerisk och en kategorisk variabel ter sig Àr det naturligt att göra ett histogram av den numeriska variabeln för varje utfall av den kategoriska variabeln. Ett sÄdant histogram Àr en betingad fördelning, eftersom endast observationerna som tillhör utfallet vi betingar pÄ för den kategoriska variabeln anvÀnds vid konstruktionen av histogrammet.
LÄt oss illustrera koncepten med ett dataset över forcerad utandningsvolym, forced expiratory volume (FEV) pÄ engelska, hos 606 barn och ungdomar. De numeriska variblerna i datasetet Àr forcerad utandningsvolym (fev), lÀngd (height) och Älder (age). De kategoriska variablerna Àr rökare (smoking, ja kodat som 1 och nej kodat som 0), Äldersgrupp (age.group) och kön (gender). Data finns i FevChildren.RData.
load(file = url("https://github.com/StatisticsSU/SDA1/blob/main/datorlab/lab3/FevChildren.RData?raw=true"))
str(FevChildren)'data.frame': 606 obs. of 6 variables:
$ fev : num 1.71 1.72 1.72 1.56 1.9 ...
$ height : num 145 171 138 135 145 ...
$ smoking : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
$ gender : Factor w/ 2 levels "f","m": 1 1 1 2 2 1 1 1 1 1 ...
$ age : num 9 8 7 9 9 8 6 6 8 9 ...
$ age.group: Factor w/ 3 levels "6-9","10-14",..: 1 1 1 1 1 1 1 1 1 1 ...head(FevChildren) fev height smoking gender age age.group
1 1.708 144.78 0 f 9 6-9
2 1.724 171.45 0 f 8 6-9
3 1.720 138.43 0 f 7 6-9
4 1.558 134.62 0 m 9 6-9
5 1.895 144.78 0 m 9 6-9
6 2.336 154.94 0 f 8 6-9đȘ Uppgift 4.1
AnvÀnd histogram() funktionen frÄn mosaic-paketet för att plotta den marginella fördelningen för fev. BerÀkna de olika fördelningsmÄtten med hjÀlp av favstats() och beskriv fördelningen.
# Write your code hereFinns det nÄgot samband mellan den forcerade utandningsvolymen (fev) och Äldersgrupp (age.group)? LÄt oss göra ett histogram för den forcerade utandningsvolymen betingat pÄ de olika Äldersgrupp samt rÀkna ut fördelningsmÄtten för att besvara frÄgan.
histogram(~ fev | age.group, data = FevChildren)
favstats(~ fev | age.group, data = FevChildren) age.group min Q1 median Q3 max mean sd n missing
1 6-9 1.165 1.72100 2.0735 2.45175 3.842 2.113778 0.4813427 270 0
2 10-14 1.458 2.57575 3.0430 3.45950 5.224 3.089574 0.6965719 296 0
3 15-17 2.198 2.99825 3.6810 4.28775 5.793 3.707000 0.9137702 40 0Vi ser att de med de största forcerade utandningsvolymen Äterfinns i den Àldsta Äldersgruppen. Vi kan ocksÄ se att bÄde medelvÀrdet och medianen ökar med Äldersgrupp. Resultaten Àr förvÀntade dÄ Àldre barn och ungdomar har ocksÄ större lungor.
đȘ Uppgift 4.2
Hur ser sambandet mellan den forcerade utandningsvolymen (fev) och rökning (smoking) ut? Ăr resultaten förvĂ€ntade?
# Write your code here.Vi kan ocksÄ betinga pÄ mer Àn en variabel. Exempelvis kan vi plotta fördelningen för fev betingat pÄ smoking och Äldersgrupp (age.group).
histogram(~ fev | smoking + age.group, data = FevChildren)
favstats(~ fev | smoking + age.group, data = FevChildren) smoking.age.group min Q1 median Q3 max mean sd n
1 0.6-9 1.165 1.7200 2.076 2.455 3.842 2.114375 0.4821395 269
2 1.6-9 1.953 1.9530 1.953 1.953 1.953 1.953000 NA 1
3 0.10-14 1.458 2.5625 2.993 3.446 5.224 3.064416 0.6890278 255
4 1.10-14 1.694 2.9530 3.152 3.498 4.789 3.246049 0.7311560 41
5 0.15-17 2.635 3.5000 3.985 4.500 5.793 4.079952 0.8827685 21
6 1.15-17 2.198 2.7370 3.122 3.763 4.872 3.294789 0.7756381 19
missing
1 0
2 0
3 0
4 0
5 0
6 0Vi sÄg tidigare att de största forcerade utandningsvolymen Äterfanns i den Àldsta Äldersgruppen. Histogrammen i första raden ovan visar att de ocksÄ fanns bland icke-rökarna. Notera att ett av utfallen, Äldersgrupp 6-9 har endast 1 observation (en alltför ung rökare). Det Àr dÄ inte sÄ meningsfullt att vare sig plotta ett histogram eller rÀkna fördelningsmÄtt.
Bland Ă„ldersgrupp 15-17 ser vi att bĂ„de medianen och medelvĂ€rdet Ă€r högre bland de som inte röker jĂ€mfört med rökarna. För Ă„ldersgrupp 10-14 Ă€r medianen och medelvĂ€rdet högre bland de som röker jĂ€mfört med de som inte röker. Betyder det att det Ă€r bĂ€ttre för lugnkapaciteten att röka om man Ă€r i Ă„ldersgruppen 10-14? Notera först att skillnaden mellan medelvĂ€rdena (och medianerna) för Ă„ldersgrupp 15-17 mellan rökare och icke-rökare Ă€r betydligt större jĂ€mfört med de motsvarande skillnaderna för Ă„ldersgrupp 10-14. Skillnaderna för Ă„ldersgrupp 10-14 Ă€r ganska smĂ„ â kan det vara sĂ„ att det Ă€r slumpfaktorn som har resulterat i högre medelvĂ€rde (och median) för rökarna (jĂ€mfört med icke-rökarna) i Ă„ldersgrupp 10-14? Hade vi tagit ett nytt stickprov pĂ„ 606 nya barn och ungdomar hade vi kanske fĂ„tt det omvĂ€nda resultatet, dvs att medelvĂ€rdet (och medianen) för rökarna i Ă„ldersgrupp 10-14 hade varit lĂ€gre Ă€n icke-rökarna i samma Ă„ldersgrupp. Ett resultat som Ă€r orsakat av slumpen sĂ€gs vara icke signifikant (statistiskt insignifikant). För att statistiskt sĂ€kerstĂ€lla att skillnaderna inte beror pĂ„ slumpen behöver vi vĂ€nta till kursens andra del, men det Ă€r viktigt att ni redan nu förstĂ„r att smĂ„ avvikelser kan bero pĂ„ slumpen och inte dra för mycket slutsatser.
Vi kan ocksÄ notera att formula-syntax med mosaicpaketet som vi har anvÀnt ovan Àr kompakt och lÀtt att implementera. Vi kan göra samma sorts figurer i base R, men det krÀver betydligt mer kod som Àr dessutom krÄngligare att förstÄ (klicka i marginalen om du Àr intresserad).
đȘ Uppgift 4.3
AnvÀnd histogram() för att plotta fördelningen för fev betingat pÄ kön (gender) och age.group. BerÀkna Àven fördelningsmÄtten med favstats() och kommentera resultaten.
# Write your code hereEtt annat sÀtt att illustrera samband mellan en numerisk variabel och kategoriska variabler Àr att anvÀnda en boxplot dÀr man betingar pÄ de kategoriska variablerna. Fördelen med denna ansats Àr att vi kan summera allt i en figur, istÀllet för att anvÀnda mÄnga figurer som vi gjorde med histogram ansatsen. Det gÄr att plotta flera histogram i samma plot med tidyverse-syntax och ocksÄ med base R, men figuren kan bli ganska stökig (klicka i marginalen om du Àr intresserad).
Följande kod plottar samma fördelning som i Uppgift 4.3, men med boxplot() istÀllet för histogram(). Som vi noterade i Avsnitt 2 sÄ anvÀnder boxplot() den alternativa formula-syntaxen, dvs istÀllet för ~ fev | gender + age.group skriver vi fev ~ gender + age.group.
boxplot(fev ~ gender + age.group, data = FevChildren)
En fördel med boxplotten kontra histogrammet Àr att alla medianer visas, dessutom i samma figur, vilket förenklar jÀmförelser mellan grupperna utan att behöva anvÀnda favstats().
đȘ Uppgift 4.4
AnvÀnd boxplot() för att jÀmföra fördelningarna för fev betingat pÄ olika utfall av smoking och age.group. Ovan hade vi en diskussion om statistisk signifikans för samma exempel. Ser ni i figuren att skillnaderna för medianerna i Äldersgrupp 10-14 mellan rökare och icke rökare Àr liten jÀmfört med samma skillnad i Äldersgrupp 15-17? Vilken av dessa skillnader Àr mest trolig att vara statistiskt signifikant?
# Write your code here5. Tidsserier
Tidsserieanalys Àr ett stort omrÄde inom statistik. I den hÀr kursen gÄr vi inte in sÄ mycket pÄ tidsserier, mer Àn att lÀr oss att plotta dom.
Filen CAPM_data.RData innehÄller tidsserier över mÄnadsavkastningar för olika finansiella tillgÄngar samt makroekonomiska variabler sÄsom obligationsrÀnta (RKFREE) och konsumentprisindex (CPI). Observationerna Àr ordnade i tid, dÀr första observationen Àr frÄn januari 1978 och sista observationen Àr frÄn december 1987. Observationerna Àr pÄ mÄnadsfrekvens, sÄ vi har totalt 120 observationer (12 per Är under 10 Är). MÄnadsavkastningen \(r_t\) för en finansiell tillgÄng som Àr vÀrderar till \(S_t\) i period \(t\) definieras som \[r_t=\frac{S_t-S_{t-1}}{S_{t-1}}.\]
Exempelvis, antag att tillgÄngen var vÀrderad till 105 vid tidpunkten \(t-1\), dvs \(S_{t-1}=100\), och 110 vid tidpunkten \(t\), dvs \(S_{t}=110\). DÄ Àr mÄnadsavkastningen i perioden \(t\)
\[r_t = \frac{110 - 105}{105} \approx 0.048,\]
dvs en uppgÄng pÄ nÀstan 5% jÀmfört med föregÄende mÄnad.
Vi lÀser in data samt inspekterar variabeltyperna.
load(file = url("https://github.com/StatisticsSU/SDA1/blob/main/datorlab/lab3/CAPM_data.RData?raw=true"))
str(CAPM)'data.frame': 120 obs. of 22 variables:
$ DATE : Date, format: "1978-01-31" "1978-02-28" ...
$ MARKET: num -0.045 0.01 0.05 0.063 0.067 ...
$ RKFREE: num 0.00487 0.00494 0.00526 0.00491 0.00513 ...
$ BOISE : num -0.079 0.013 0.07 0.12 0.071 ...
$ CONED : num -0.079 -0.003 0.022 -0.005 -0.014 ...
$ CONTIL: num -0.129 0.037 0.003 0.18 0.061 ...
$ CITCRP: num -0.115 -0.019 0.059 0.127 0.005 ...
$ DATGEN: num -0.084 -0.097 0.063 0.179 0.052 ...
$ DEC : num -0.1 -0.063 0.01 0.165 0.038 ...
$ DELTA : num -0.028 -0.033 0.07 0.15 -0.031 ...
$ GENMIL: num -0.099 0.018 -0.023 0.046 0.063 ...
$ GERBER: num -0.048 0.16 -0.036 0.004 0.046 ...
$ IBM : num -0.029 -0.043 -0.063 0.13 -0.018 ...
$ MOBIL : num -0.046 -0.017 0.049 0.077 -0.011 ...
$ PANAM : num 0.025 -0.073 0.184 0.089 0.082 ...
$ PSNH : num -0.008 -0.025 0.026 -0.008 0.019 ...
$ TANDY : num -0.075 -0.004 0.124 0.055 0.176 ...
$ TEXACO: num -0.054 -0.01 0.015 0 -0.029 ...
$ WEYER : num -0.116 -0.135 0.084 0.144 -0.031 ...
$ POIL : num 7.9 7.87 7.79 7.86 7.89 ...
$ FRBIND: num 126 128 128 129 130 ...
$ CPI : num 167 167 168 168 169 ...Notera att dataframen som laddas in heter CAPM och inte CAPM_data. Det finns ingen koppling mellan filnamnet och namnet pÄ dataframen i R. Namnet pÄ dataframen bestÀms av vad den döptes till nÀr den sparades, i det hÀr fallet CAPM.
Vi noterar att alla variabler förutom DATE Àr numeriska. DATE Àr i sÄkallad Date-format, vilket gör att R förstÄr att det Àr ett datum, trots att den verkar lagras som strÀngar enligt utskriften ovan. R behandlar alltsÄ inte DATE som en kategorisk variabel.
Datasetet innehÄller variabeln MARKET, vilket Àr marknadsportföljensavkastning (dvs alla tillgÄngar i marknaden). Vi kan göra en tidsserieplot för marknadsportföljens avkastning.
plot(MARKET ~ DATE, data = CAPM, col = "darkblue")
Eftersom det vi plottar Ă€r en tidsserie sĂ„ Ă€r observationerna ordnade och dĂ„ Ă€r det vanligt att man binder ihop observationerna med linjer, dvs en linjeplot. Detta kan vi göra genom att anvĂ€nda argumentet type = "l" i funktionen plot(), dĂ€r âlâ stĂ„r för line.
plot(MARKET ~ DATE, data = CAPM, col = "darkblue", type = "l")
đȘ Uppgift 5.1
Om vi tycker att linjen Àr för tunn kan vi kontrollera dess tjocklek med plot() funktionens argument lwd som stÄr för linewidth. DefaultvÀrdet (dvs vÀrdet som anvÀnds om inget annat anges) Àr lwd = 1. Om vi vill ha en linje som Àr 50% tjockare anger vi lwd = 1.5. Gör om plotten ovan med en dubbelt sÄ tjock linje.
# Write your code hereFrÄn tidsserieplotten ser vi ser att, över en lÄng period, sÄ kan marknaden ha ganska kraftiga svÀngningar.
đȘ Uppgift 5.2
Black Monday Àr namnet pÄ en av de vÀrsta börskrascherna nÄgonsin. Black Monday intrÀffade oktober 1987. Kan ni se Black Monday hÀndelsen i tidsserien? Hur mycket föll marknadsportföljen den mÄnaden?
# Write your code here
Tips
Ibland kan det vara bra att kunna hÀmta ut en variabel frÄn en dataframe. Antag att vi har en dataframe my_data med en variabel vid namn x. Vi kan anvÀnda $-tecken för att hÀmta ut den variabeln frÄn dataframen. Vi skriver dÄ my_data$x.
đȘ Uppgift 5.3
Plotta tidsserien för tillgÄngen CITCRP. Hur drabbades CITCRP av Black Monday jÀmfört med marknadsportföljen? Skriv ut return för CITCRP och jÀmför med Uppgift 5.2.
# Write your code hereNotera att varje gÄng vi anropar plot() funktionen sÄ skapas en ny figur. Ibland vill man plotta flera tidsserierna i samma graf. Detta kan Ästadkommas genom att skapa första figuren med plot(), precis som vi gjort innan, och efterföljande figurer med lines(). lines() stödjer inte formula-syntax, sÄ det hÀr Àr ytterligare ett exempel pÄ varför det Àr viktigt att kunna hÀmta ut variabler frÄn en dataframe med hjÀlp av $-tecknet. Följande kod plottar tidsserierna för marknadsportföljen och CITCRP i samma graf. Koden lÀgger ocksÄ till etiketter via legend() funktionen och justerar y-axelns definitionsomrÄde via argumentet ylim, som ni ocksÄ anvÀnde i Lab 2.
plot(MARKET ~ DATE, data = CAPM, col = "darkblue", type = "l", lwd = 2, ylim = c(-0.4, 0.4))
lines(x = CAPM$DATE, y = CAPM$CITCRP, col = "lightblue", type = "l", lwd = 2)
legend(x = "topright", lty = 1 , lwd = 2, col = c("darkblue", "lightblue"), legend=c("Market", "CITCRP"))
đȘ Uppgift 5.4
Plotta de tre tidsserierna MARKET, CITCRP och DATGEN i samma figur. AnvÀnd lÀmpliga fÀrger samt etiketter. Verkar tidsserierna följas Ät? Vad tror ni att det beror pÄ?
# Write your code heređȘ Uppgift 5.5
Antag att ni skapar en portfölj som bestÄr av 30% IBM, 50% CITCRPoch 20% DATGEN. Plotta mÄnadsavkastningen för er portfölj.
Tips
Avkastningen för en portfölj Àr en viktad linjÀr kombination av avkastningarna i tillgÄngarna som portföljen innehÄller. I en portfölj med tre aktier vars vikter \(0<w_1<1, 0<w_2<1, 0<w_3<1\), med \(\sum_{k=1}^3w_k=1\), Àr portföljens avkastning \[r_t = w_1r^{(1)}_t + w_2r^{(2)}_t+w_3r^{(3)}_t,\]
dÀr \(r^{(k)}_t\) Àr avkastningen för tillgÄng \(k\) i portföljen.
Nu nÀr ni vet hur man kan anvÀnda $-tecken för att hÀmta ut variabler frÄn en dataframe, kan ni lista ut hur man skapar en ny variabel som viktar variablerna enligt ovan (30% IBM, 50% CITCRPoch 20% DATGEN).
6. Sammanfattning
A. Appendix
Lista över nÄgra vanliga argument i grafer
col = fÀrg, kan markeras med siffror eller med namnet pÄ fÀrgen, oftast med smÄ bokstÀver. Glöm ej att texten mÄste ligga inom citattecken.
main = rubrik pÄ plotten, detta sÀtts till en textstrÀng, alltsÄ en text inom citattecken.
xlab = rubrik pÄ x-axeln, detta sÀtts till en textstrÀng, alltsÄ en text inom citattecken.
ylab = rubrik pÄ y-axeln, detta sÀtts till en textstrÀng, alltsÄ en text inom citattecken.
xlim = definitionsomrÄde för x-axeln. Exempelvis betyder xlim = c(0, 14.3) att det lÀgsta vÀrdet som ritas kommer att vara 0 och det högsta vÀrdet kommer att vara 14.3 (pÄ x-axeln).
ylim = definitionsomrÄde för y-axeln. Exempelvis betyder ylim = c(-2, 20.7) att det lÀgsta vÀrdet som ritas kommer att vara -2 och det högsta vÀrdet kommer att vara 20.7 (pÄ y-axeln).
lwd = tjocklecken pÄ linjer (eller prickar) i ett linjediagram (spridningsdiagram), anges med ett nummer.
lty = typ av linje (rak, streckad, prickad etc) i ett linjediagram, anges med ett nummer.
pch = typ av prick (rund, fyrkantig + * etc), sÀtts till ett nummer.
breaks = antal staplar i ett histogram.
cex = justering av etiketternas storlek i en figur (exempelvis cex = 0.5 i en
legend()funktion minskar storleken som âlegend-textenâ tar upp i en figur med hĂ€lften). Finns Ă€ven exempelvis cex.axis, cex.main för att justera storleken av texter pĂ„ axlar respektive rubrik.bg = bakgrundsfĂ€rg i en figur.
col.main = rubrikens fÀrg.
col.lab = fÀrger för rubrikerna pÄ axlarna.
font = specificerar vilken typ av text man vill ha, exempelvis ger font = 3 kursiv text.