Instruktioner

Skriv din kode indenfor de forberedte kode-chunks og sikr undervejs, at dit dokument kan knit.
Det kan være nyttigt at bruge in-line kode i nogle af dine sætninger.
Husk at give dine plots passende etiketter som akseetiketter, en titel og et tema.
I er mere end velkommen til at hjælpe hinanden, men I skal aflevere hver jeres egen opgave.
Husk at aflevere dit HTML-dokument til FeedbackFruits, når du er færdig (se slides fra i dag for yderligere instruktioner om det).
- Efter afleveringsperioden er overstået, vil du modtage løsningerne, som du kan bruge til at give feedback til en anden person i FeedbackFruits.

Indlæs følgende pakker:

#notice warning=F and comment=F to stop 'noisy' package output upon knit
library(ggplot2)
library(ggrepel)
library(tidyverse)

FALSE ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
FALSE ✔ dplyr     1.1.4     ✔ readr     2.1.5
FALSE ✔ forcats   1.0.0     ✔ stringr   1.5.1
FALSE ✔ lubridate 1.9.4     ✔ tibble    3.2.1
FALSE ✔ purrr     1.0.4     ✔ tidyr     1.3.1
FALSE ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
FALSE ✖ dplyr::filter() masks stats::filter()
FALSE ✖ dplyr::lag()    masks stats::lag()
FALSE ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Opgave 1

Vi arbejder med datasættet msleep (mammals sleep) - se venligst ?msleep for yderlige informationer.

Kør følgende kode for at læse datasættet ind og fjerne nogle observationer med manglende værdier:

data(msleep)
?msleep #se beskrivelserne af kolonner i datasættet

## starting httpd help server ... done

msleep <- msleep[!is.na(msleep$vore),] #rengøre datasættet ved at fjerne manglende værdier i variablen `vore`
head(msleep) #se de første 6 rækker

## # A tibble: 6 × 11
##   name    genus vore  order conservation sleep_total sleep_rem sleep_cycle awake
##   <chr>   <chr> <chr> <chr> <chr>              <dbl>     <dbl>       <dbl> <dbl>
## 1 Cheetah Acin… carni Carn… lc                  12.1      NA        NA      11.9
## 2 Owl mo… Aotus omni  Prim… <NA>                17         1.8      NA       7  
## 3 Mounta… Aplo… herbi Rode… nt                  14.4       2.4      NA       9.6
## 4 Greate… Blar… omni  Sori… lc                  14.9       2.3       0.133   9.1
## 5 Cow     Bos   herbi Arti… domesticated         4         0.7       0.667  20  
## 6 Three-… Brad… herbi Pilo… <NA>                14.4       2.2       0.767   9.6
## # ℹ 2 more variables: brainwt <dbl>, bodywt <dbl>

1.1) Lav et barplot / stolpediagram for at visualisere antallet af dyr opdelt efter vore (brug fill=vore).

Ændr rækkefølgen på kolonnerne (hint: scale_x_discrete()), således at søjlen med det mindste antal dyr er til venstre og søjlen med det det største antal er til højre.
Du kan fjerne legenden, hvis du synes, at den ikke er nødvendig.

msleep$vore <- factor(msleep$vore, 
                      levels = names(sort(table(msleep$vore), decreasing = TRUE)))

ggplot(data = msleep , aes(x = vore, fill=vore)) +
  geom_bar(position = "dodge") +
  theme_minimal() + 
  scale_x_discrete()

1.2) Lav et scatterplot / punktdiagram af hjernestørrelse (brainwt) på x-aksen og total sleep (sleep_total) på y-aksen.

Transformer x-aksen til \(\log_{2}\)-skalen, så at man bedre kan visualisere datapunkterne på plotområdet.
Giv hvert niveau af variablen vore sin egen farve og punkt form.

OBS! Her får man en advarsel om 25 manglende værdier - se bare bort fra det!

ggplot(data = msleep , aes(x =  log2(brainwt), y = sleep_total, color =vore)) +
  geom_point() +
  xlab("brainwt") + 
  ylab("total sleep") + 
  theme_minimal()

## Warning: Removed 25 rows containing missing values or values outside the scale range
## (`geom_point()`).

1.3) For delmængden af datasættet med kun omnivores (dvs. “omni” fra variablen vore), lav samme scatterplot / punktdiagram som i 1.2) med hjernestørrelse (brainwt) på x-aksen og total sleep (sleep_total) på y-aksen.

Tilføj labels / etiketter direkte på plottet med de forskellige dyre navne (variablen name).
Brug annotate til at tilføje en orange rektangel omkring Chimpanzee og Human.
Skriv en kort sætning om hvordan Human sammenlignes med Chimpanzee med hensyn til hjernevægt og søvnbehov.

msleep_sub <- subset(msleep, vore == "omni")

ggplot(data = msleep_sub , aes(x =  log2(brainwt), y = sleep_total, color =vore)) +
  geom_point() +
  xlab("brainwt") + 
  ylab("total sleep") + 
  theme_minimal() + 
  geom_text_repel(aes(label = name)) + 
   annotate("rect", xmin = -2, xmax = 0.5, ymin = 7, ymax = 10, 
           color = "orange", fill = NA)

## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_text_repel()`).

Chimpanzens hjernevægt er den tættest på menneskets, og deres mængde af søvn er næsten den samme.

Opgave 2

Vi arbejder med LungCapData, som har følgende variabler:

LungCap - lungekapacitet
Age - i antal år
Height - i tommer
Smoke - ryger/ikke-ryger
Gender - male/female
Caesarean - fået kejsersnit no/yes

Datasættet findes i filen LungCapData.csv, som er tilgængelig på Absalon, og man kan bruge funktionen read.csv til at åbne filen. Bemærk, at det også er muligt (og nemmere!) at indlæse det samme datasæt ved hjælp af Dropbox-linket nedenfor.

#indlæs de data her, enten ved read.csv (data fil er på Absalon) eller med følgende link

#LungCapData <- read.csv(...)

LungCapData <- read.csv("https://www.dropbox.com/s/ke27fs5d37ks1hm/LungCapData.csv?dl=1")
head(LungCapData) #se variabler navne

##   LungCap Age Height Smoke Gender Caesarean
## 1   6.475   6   62.1    no   male        no
## 2  10.125  18   74.7   yes female        no
## 3   9.550  16   69.7    no female       yes
## 4  11.125  14   71.0    no   male        no
## 5   4.800   5   56.9    no   male        no
## 6   6.225  11   58.7    no female        no

Kør også følgende chunk for at tilføje en variabel, der hedder Age.Group, som vi bruger i analysen (selve koden er ikke vigtig men jeg benytter funktionen cut for at opdele den kontinuerlige variabel Age i fire aldersgrupper).

#Her laver jeg nogle alder grupper som kan bruges i analysen
LungCapData$Age.Group <- cut(LungCapData$Age,breaks=c(1,13,15,17,19),right=FALSE,include.lowest = TRUE)
levels(LungCapData$Age.Group) <- c("<13","13-14","15-16","17+")

2.1) Lav et barplot / stolpediagram for at undersøge om antallet af personer der ryger er forskellige for “males” og “females”.

LungCapData %>%
  filter(Smoke == "yes") %>% 
ggplot(aes(x = Smoke, fill= Gender)) +
  geom_bar(position = "dodge") +
  theme_minimal()

2.2) Lav et density plot / tæthedsplot for variablen LungCap og skriv en kort kommentar om fordelingen - ligner den tilnærmelsesvis en normalfordeling?

LungCapData %>%
ggplot(aes(x = LungCap, fill= Gender)) +
  xlab("LungCap") +
  theme_minimal() + 
  geom_density(alpha = 0.5)

Hvis vi antager at middelværdien er ca 7.5, med en spredning på 2.5 så ligner det meget en normalfordeling.

2.3) a) Lav et boxplot / boksplot for LungCap opdelt efter rygning-status.

LungCapData %>%
ggplot(aes(x = Smoke, y = LungCap, fill= Smoke)) +
  xlab("Smoker") +
  theme_minimal() + 
  geom_boxplot()

b) Er der en signifikant forskel i variablen LungCap mellem rygere og ikke-rygere? Hint: anvend t.test (?t.test) til at sammenligne LungCap for rygere med LungCap for ikke-rygere og skriv p-værdien ned.

t.test(LungCap ~ Smoke, data = LungCapData)

## 
##  Welch Two Sample t-test
## 
## data:  LungCap by Smoke
## t = -3.6498, df = 117.72, p-value = 0.0003927
## alternative hypothesis: true difference in means between group no and group yes is not equal to 0
## 95 percent confidence interval:
##  -1.3501778 -0.4003548
## sample estimates:
##  mean in group no mean in group yes 
##          7.770188          8.645455

p-value = 0.0003927

Der er signifikant forskel mellem rygere og ikke rygeres lungekapacitet i dette datasæt.

2.4) a) Lav det samme boxplot / boksplot af LungCap opdelt efter rygestatus, men kun for personer som er over 16 år gamle.

LungCapData %>%
  filter(Age > 16) %>%
  ggplot(aes(x = Smoke, y = LungCap, fill = Smoke)) +
  geom_boxplot() +
  theme_minimal()

b) Er der en signifikant forskel i LungCap mellem rygere og ikke-rygere blandt personer, der er over 16 år gamle? Får du det samme forhold som i 2.3)?

t.test(LungCap ~ Smoke, data = subset(LungCapData, Age > 16))

## 
##  Welch Two Sample t-test
## 
## data:  LungCap by Smoke
## t = 2.9447, df = 47.129, p-value = 0.005008
## alternative hypothesis: true difference in means between group no and group yes is not equal to 0
## 95 percent confidence interval:
##  0.2720388 1.4450112
## sample estimates:
##  mean in group no mean in group yes 
##          11.06237          10.20385

p-value = 0.005008

Ja der er stadig en signifikant forskel, men p værdien er anderledes, naturligvis.

2.5) Vi vil gerne undersøge vores observationer fra de to foregående spørgsmål nærmere.

a) Opret de samme boxplots / boksplotter af LungCap for alle personer, men brug facet_grid til at adskille efter Age.Group.

b) Tilføj geom_jitter() for at få punkterne til at ligge ovenpå boxplots / boksplotter.

c) Kan du skrive din fortolkning / idé ned for hvorfor vi så modsatte forhold i de to foregående plots?

ggplot(LungCapData, aes(x = Smoke, y = LungCap, fill = Smoke)) +
  geom_boxplot() +
  facet_grid(. ~ Age.Group) + 
  theme_minimal() +
  geom_jitter() +
  xlab("Smoker")

# Der er større standard afvigelse pga flere data punkter hos none smokers.

2.6) Vi vil gerne undersøge, om vores statistiske resultater stemmer overens med vores observationer i de tidligere spørgsmål. Kør følgende to lineære modeller - én model der ikke tager højde for variablen Age (til sammenligning med 2.3) og én model der inddrager Age (til sammenligning med 2.5).

# Effekten af rygning på lungkapacitet
mylm_Smoke <- lm(LungCap ~ Smoke, data = LungCapData)

# Effekten af både alder og rygning på lungkapacitet
mylm_Age_Smoke <- lm(LungCap ~ Age + Smoke, data = LungCapData)

a) Lav et summary() for begge modeller og angiv Estimate for Smokeyes (effekten af rygning på lungkapaciteten).

b) Stemmer fortegnet på dine estimeringer overens med dine observationer i 2.3) og 2.5)?

c) Tilføj Height til formlen i modellen mylm_Age_Smoke med +-tegnet. Hvor meget af variansen i LungCap forklarer de tre variable samlet set (kig på multiple R-squared værdien)?

summary(mylm_Smoke)

## 
## Call:
## lm(formula = LungCap ~ Smoke, data = LungCapData)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.2632 -1.7202  0.1048  1.9048  6.9048 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   7.7702     0.1041   74.64   <2e-16 ***
## Smokeyes      0.8753     0.3194    2.74   0.0063 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.65 on 723 degrees of freedom
## Multiple R-squared:  0.01028,    Adjusted R-squared:  0.008908 
## F-statistic: 7.507 on 1 and 723 DF,  p-value: 0.006297

#smoke yes estimate når man kun sammenligner Lungcap med smoke or no smoke,  Smokeyes = 0.8753

summary(mylm_Age_Smoke)

## 
## Call:
## lm(formula = LungCap ~ Age + Smoke, data = LungCapData)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.8559 -1.0289 -0.0363  1.0083  4.1995 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.08572    0.18299   5.933 4.61e-09 ***
## Age          0.55540    0.01438  38.628  < 2e-16 ***
## Smokeyes    -0.64859    0.18676  -3.473 0.000546 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.514 on 722 degrees of freedom
## Multiple R-squared:  0.6773, Adjusted R-squared:  0.6764 
## F-statistic: 757.5 on 2 and 722 DF,  p-value: < 2.2e-16

#Smoke yes estimate når Alder bliver taget i betragtning, Smokeyes -0.64859

#De stemmer ikke overens fordi vi ikke tager højde for variablen "alder", så når vi i 2.3 undersøgte om der var en sammenhæng mellem røgning og lungcap så er den højere for dem der ryger. Det er fordi dem som er unge på måske 13 år har en gennemsnitlig højere lungecap end dem som ikke gør, og der er mange færre data points i smokesyes så de punkter er meget mere statisically significant. Når alder kommer ind i billedet så bliver estimate -negativt for smokeyes.

mylm_Age_Smoke_Height <- lm(LungCap ~ Age + Smoke + Height, data = LungCapData)
summary(mylm_Age_Smoke_Height)

## 
## Call:
## lm(formula = LungCap ~ Age + Smoke + Height, data = LungCapData)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4840 -0.7135  0.0358  0.7088  3.0392 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -11.808195   0.469120 -25.171  < 2e-16 ***
## Age           0.136916   0.017677   7.745 3.24e-14 ***
## Smokeyes     -0.648583   0.128099  -5.063 5.24e-07 ***
## Height        0.278432   0.009761  28.525  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.039 on 721 degrees of freedom
## Multiple R-squared:  0.8484, Adjusted R-squared:  0.8477 
## F-statistic:  1345 on 3 and 721 DF,  p-value: < 2.2e-16

#Vores estimate fortæller os at height har størst indflydelse og Age er derefter, Smokeyes påvirker lungcap negativt.

Workshop 1 - ggplot2

Den 2. maj 2025

Opgave 1

Opgave 2