EXERCICIS

Exercici 1

Haureu de:
a) Buscar un resum estadístic de les variables del dataset Iris i Orange.
data("iris")
summary(iris) # Busquem el resum estadístic
##   Sepal.Length    Sepal.Width     Petal.Length    Petal.Width   
##  Min.   :4.300   Min.   :2.000   Min.   :1.000   Min.   :0.100  
##  1st Qu.:5.100   1st Qu.:2.800   1st Qu.:1.600   1st Qu.:0.300  
##  Median :5.800   Median :3.000   Median :4.350   Median :1.300  
##  Mean   :5.843   Mean   :3.057   Mean   :3.758   Mean   :1.199  
##  3rd Qu.:6.400   3rd Qu.:3.300   3rd Qu.:5.100   3rd Qu.:1.800  
##  Max.   :7.900   Max.   :4.400   Max.   :6.900   Max.   :2.500  
##        Species  
##  setosa    :50  
##  versicolor:50  
##  virginica :50  
##                 
##                 
## 
data("Orange")
summary(Orange) # Busquem el resum estadístic
##  Tree       age         circumference  
##  3:7   Min.   : 118.0   Min.   : 30.0  
##  1:7   1st Qu.: 484.0   1st Qu.: 65.5  
##  5:7   Median :1004.0   Median :115.0  
##  2:7   Mean   : 922.1   Mean   :115.9  
##  4:7   3rd Qu.:1372.0   3rd Qu.:161.5  
##        Max.   :1582.0   Max.   :214.0
b) Generar una taula de freqüències absolutes i una taula de freqüències relatives amb el dataset Iris. Totes les taules generades tenen sentit per a vosaltres?

En aquest cas, de les taules generades, només la de les espècies (variable categòrica) ens aporta informació ràpidament, indicant que de la mostra total, hi ha una repartició igual entre les tres espècies. Pel que fa a les altres quatre variants, no ens aporten una informació tant ràpid (ja que son variables numèriques contínues), tot i que es pot observar, principalment en Sepal.Length i .Width, que hi ha una distribució normal (campana de Gauss) en les mesures. Tot i això, en aquest cas seria molt més útil representar les dades en un histograma que no pas en una taula.

# Carreguem les dades
data("iris")

# Carreguem la llibreria knitr per usar la funció kable
library(knitr) # kable s'usa per millorar el format de la taula, doncs sense ella,
# els resultats es presenten de forma poc ordenada

# Generem taules per cadascuna de les variables d'iris
kable(table(iris$Species))
Var1 Freq
setosa 50
versicolor 50
virginica 50 
kable(table(iris$Sepal.Length))
Var1 Freq
4.3 1
4.4 3
4.5 1
4.6 4
4.7 2
4.8 5
4.9 6
5 10
5.1 9
5.2 4
5.3 1
5.4 6
5.5 7
5.6 6
5.7 8
5.8 7
5.9 3
6 6
6.1 6
6.2 4
6.3 9
6.4 7
6.5 5
6.6 2
6.7 8
6.8 3
6.9 4
7 1
7.1 1
7.2 3
7.3 1
7.4 1
7.6 1
7.7 4
7.9 1 
kable(table(iris$Sepal.Width))
Var1 Freq
2 1
2.2 3
2.3 4
2.4 3
2.5 8
2.6 5
2.7 9
2.8 14
2.9 10
3 26
3.1 11
3.2 13
3.3 6
3.4 12
3.5 6
3.6 4
3.7 3
3.8 6
3.9 2
4 1
4.1 1
4.2 1
4.4 1 
kable(table(iris$Petal.Length))
Var1 Freq
1 1
1.1 1
1.2 2
1.3 7
1.4 13
1.5 13
1.6 7
1.7 4
1.9 2
3 1
3.3 2
3.5 2
3.6 1
3.7 1
3.8 1
3.9 3
4 5
4.1 3
4.2 4
4.3 2
4.4 4
4.5 8
4.6 3
4.7 5
4.8 4
4.9 5
5 4
5.1 8
5.2 2
5.3 2
5.4 2
5.5 3
5.6 6
5.7 3
5.8 3
5.9 2
6 2
6.1 3
6.3 1
6.4 1
6.6 1
6.7 2
6.9 1 
kable(table(iris$Petal.Width))
Var1 Freq
0.1 5
0.2 29
0.3 7
0.4 7
0.5 1
0.6 1
1 7
1.1 3
1.2 5
1.3 13
1.4 8
1.5 12
1.6 4
1.7 2
1.8 12
1.9 5
2 6
2.1 6
2.2 3
2.3 8
2.4 3
2.5 3 
# I per les mateixes variables generem la taula de freqüències relativa
kable(prop.table(table(iris$Species)))
Var1 Freq
setosa 0.3333333
versicolor 0.3333333
virginica 0.3333333 
kable(prop.table(table(iris$Sepal.Length)))
Var1 Freq
4.3 0.0066667
4.4 0.0200000
4.5 0.0066667
4.6 0.0266667
4.7 0.0133333
4.8 0.0333333
4.9 0.0400000
5 0.0666667
5.1 0.0600000
5.2 0.0266667
5.3 0.0066667
5.4 0.0400000
5.5 0.0466667
5.6 0.0400000
5.7 0.0533333
5.8 0.0466667
5.9 0.0200000
6 0.0400000
6.1 0.0400000
6.2 0.0266667
6.3 0.0600000
6.4 0.0466667
6.5 0.0333333
6.6 0.0133333
6.7 0.0533333
6.8 0.0200000
6.9 0.0266667
7 0.0066667
7.1 0.0066667
7.2 0.0200000
7.3 0.0066667
7.4 0.0066667
7.6 0.0066667
7.7 0.0266667
7.9 0.0066667 
kable(prop.table(table(iris$Sepal.Width)))
Var1 Freq
2 0.0066667
2.2 0.0200000
2.3 0.0266667
2.4 0.0200000
2.5 0.0533333
2.6 0.0333333
2.7 0.0600000
2.8 0.0933333
2.9 0.0666667
3 0.1733333
3.1 0.0733333
3.2 0.0866667
3.3 0.0400000
3.4 0.0800000
3.5 0.0400000
3.6 0.0266667
3.7 0.0200000
3.8 0.0400000
3.9 0.0133333
4 0.0066667
4.1 0.0066667
4.2 0.0066667
4.4 0.0066667 
kable(prop.table(table(iris$Petal.Length)))
Var1 Freq
1 0.0066667
1.1 0.0066667
1.2 0.0133333
1.3 0.0466667
1.4 0.0866667
1.5 0.0866667
1.6 0.0466667
1.7 0.0266667
1.9 0.0133333
3 0.0066667
3.3 0.0133333
3.5 0.0133333
3.6 0.0066667
3.7 0.0066667
3.8 0.0066667
3.9 0.0200000
4 0.0333333
4.1 0.0200000
4.2 0.0266667
4.3 0.0133333
4.4 0.0266667
4.5 0.0533333
4.6 0.0200000
4.7 0.0333333
4.8 0.0266667
4.9 0.0333333
5 0.0266667
5.1 0.0533333
5.2 0.0133333
5.3 0.0133333
5.4 0.0133333
5.5 0.0200000
5.6 0.0400000
5.7 0.0200000
5.8 0.0200000
5.9 0.0133333
6 0.0133333
6.1 0.0200000
6.3 0.0066667
6.4 0.0066667
6.6 0.0066667
6.7 0.0133333
6.9 0.0066667 
kable(prop.table(table(iris$Petal.Width)))
Var1 Freq
0.1 0.0333333
0.2 0.1933333
0.3 0.0466667
0.4 0.0466667
0.5 0.0066667
0.6 0.0066667
1 0.0466667
1.1 0.0200000
1.2 0.0333333
1.3 0.0866667
1.4 0.0533333
1.5 0.0800000
1.6 0.0266667
1.7 0.0133333
1.8 0.0800000
1.9 0.0333333
2 0.0400000
2.1 0.0400000
2.2 0.0200000
2.3 0.0533333
2.4 0.0200000
2.5 0.0200000
c) Generar una taula de freqüències absolutes amb cadascuna de les variables del conjunt de dades Orange. Totes les taules generades tenen sentit per a vosaltres?

No, no totes les taules tenen sentit. Si bé les resultants d’analitzar ‘Tree’ i ‘Age’ poden ser útils (Tree és una variable categòrica, que indica que tenim 5 arbres i cada un ha sigut mesurat 7 cops; i Age, que indica els 7 dies en que es van prendre mesures, i cada dia es van fer mesures a 5 arbres). La taula que representa els resultats de circumferència no aporta cap informació addicional. Al ser una variable numèrica i continua, cada arbre té un valor diferent i rarament es repeteixen, pel que per comptes d’agrupar els arbres entre si, els deixa tal i com estan.

data("Orange")

kable(table(Orange$Tree))
Var1 Freq
3 7
1 7
5 7
2 7
4 7 
kable(table(Orange$age))
Var1 Freq
118 5
484 5
664 5
1004 5
1231 5
1372 5
1582 5 
kable(table(Orange$circumference))
Var1 Freq
30 3
32 1
33 1
49 1
51 1
58 1
62 1
69 1
75 1
81 1
87 1
108 1
111 1
112 1
115 2
120 1
125 1
139 1
140 1
142 2
145 1
156 1
167 1
172 1
174 1
177 1
179 1
203 2
209 1
214 1
d) Generar una taula de doble entrada entre les variables Tree i Age d’Orange.
data("Orange")

kable(table(Orange$Tree, Orange$age))
118 484 664 1004 1231 1372 1582
3 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1
5 1 1 1 1 1 1 1
2 1 1 1 1 1 1 1
4 1 1 1 1 1 1 1

Exercici 2

vect1 <- c(1,2,1,2,1,2,1,2,1,2,1,1,1,1,2,2,1,1,2,1)

vect2 <- c(1,1,2,2,2,1,2,1,1,2,1,2,1,1,1,2,1,1,1,1)

a) Usant el vect1 creem un nou vector anomenat Baix_pes que sigui un factor amb dos nivells. Les etiquetes corresponen a 1 = Baix pes i 2 = Pes normal.
b) Usant el vect2 creem un nou vector anomenat Fumador que sigui un factor amb dos nivells. Les etiquetes corresponen a 1 = Fuma i 2 = No fuma.
c) Creem una taula de contingència amb les dues variables anteriors amb el nom Taula.
d) Mirem la relació de les variables anteriors amb la prova de khi quadrat.
e) Mirem també com resulta el test de Fisher.
# a)
vect1 <- c(1,2,1,2,1,2,1,2,1,2,1,1,1,1,2,2,1,1,2,1)
vect2 <- c(1,1,2,2,2,1,2,1,1,2,1,2,1,1,1,2,1,1,1,1)

Baix_pes <- factor(vect1,levels=c("1","2"),labels=c("Baix Pes", "Pes Normal"))

# b)
Fumador <- factor(vect2,levels=c("1","2"),labels=c("Fuma","No Fuma"))

# c)
Taula <- table(Baix_pes,Fumador)
Taula
##             Fumador
## Baix_pes     Fuma No Fuma
##   Baix Pes      8       4
##   Pes Normal    5       3
# d)
chisq.test(Taula)
## Warning in chisq.test(Taula): Chi-squared approximation may be incorrect
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  Taula
## X-squared = 0, df = 1, p-value = 1
# e)
fisher.test(Taula)
## 
##  Fisher's Exact Test for Count Data
## 
## data:  Taula
## p-value = 1
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##   0.1200513 10.9278345
## sample estimates:
## odds ratio 
##   1.189031

Havent obtingut els resultats, es pot determinar que no hi ha motius suficients per demostrar una relació entre el fet de fumar i una disminució de pes. Un p-value = 1 ens indica que no podem rebutjar la hipòtesi nul·la, i per tant, per la mostra actual diriem que el pes i el tabaquisme son variables independents.

Exercici 3

a) Activeu el paquet de dades airquality del paquet datasets i genereu els gràfics següents:
• Un gràfic de dispersió de la variable Ozone de color blau i amb coixinets (#) en comptes de punts.
• Un gràfic de caixa de color vermell amb la variable Temp amb el títol Temperatura (en graus Farenheit).
b) Activeu el paquet de dades airmiles del paquet datasets i genereu els gràfics següents:
• Un gràfic de línies de la sèrie de dades airmiles amb el títol Dades de passatgers en vols comercials (en milers) i de color blau (cadetblue2) i l’etiqueta de l’eix x amb Milers de passatgers.
• Un histograma de la sèrie de dades airmiles de color marró (xocolata2).
c) Representeu els quatre gràfics en una única imatge on els vegem junts.
# c) Per a generar els 4 gràfics junts hem de posar aquest fragment de codi a l'inici com a un canvas que s'anirà omplint amb els diferents gràfics que anem creant
par(mfrow=c(2,2)) # Informació agafada de: https://nathanieldphillips-yarrr.share.connect.posit.cloud/arranging-plots-with-parmfrow-and-layout.html

# a)
data("airquality")

# Generem el primer punt, de la variable Ozone amb els valors representats mitjançant coixinets (#)
plot(airquality$Ozone, col="cadetblue2", pch="#")

# Generem el segon punt, el gràfic de caixa de la variable Temp
boxplot(airquality$Temp, col="red", main="Temperatura (en graus Fahrenheit)")


# b) 
data("airmiles")

# Generem el gràfic de línies per airmiles. En l'eix de les x deixem el títol que ens indiquen, però sembla que representa els anys. 
plot(airmiles, col="cadetblue2", main="Dades de passatgers en vols comercials (en milers)", xlab="Milers de passatgers", lwd=3) # La línia era molt prima, així que seguint els apunts he fet servir la variable lwd, per augmentar el gruix

# Generem l'histograma
hist(airmiles, col="chocolate2")

Exercici 4

Intenteu reproduir l’exemple 7 d’aquest LAB amb unes dades simulades per vosaltres. No fa falta que siguin semblants, però cal que pugueu fer diferents gràfics estadístics amb l’ordre plot().
set.seed(909)

# Simulem les variables 
x1 <- rnorm (400)
y1 <- ((4.32 * x1) / 10) +rnorm (400)

# Generem el gràfic de dispersió que conté una recta de regressió
plot(x1, y1, main="Gràfic de dispersió")
abline(lm(y1~x1), col="chocolate4")

# Generem el gràfic de densitats per x1 i y1
plot(density(x1), main = " Densitat d'X1")
lines(density(y1), col="blue")

Exercici 5

a) Amb les dades airquality del paquet datasets creeu un gràfic de dispersió simple, però amb la recta de regressió ajustada de color vermell (darkred).
b) Amb les dades airquality del paquet datasets creeu un gràfic de dispersió amb les variables Solar.R i Temp amb les seves respectives etiquetes d’eix. Els colors i la forma dels punts han de ser en funció de la variable Month convertida a factor.
c) Amb les dades de Birthwt del paquet MASS mostreu la distribució de la variable age que tingui la vora de color negre i blau a l’interior de les barres.
d) Amb les dades de Birthwt del paquet MASS mostreu la distribució de la variable age, però separada en dos gràfics (fumador/no fumador) per la variable smoke.
# Carreguem la llibreria tidyverse i activem les dades
library(tidyverse) 
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.0     ✔ readr     2.2.0
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.2     ✔ tibble    3.3.1
## ✔ lubridate 1.9.5     ✔ tidyr     1.3.2
## ✔ purrr     1.2.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
data("airquality")
 
 # a) 
 ggplot(data = airquality, aes(x = Ozone, y = Temp)) + geom_point(size=2) + geom_smooth(method = "lm", color = "darkred") + labs(title = "Relació entre l'Ozó i la Temperatura") # La funció labs ens serveix per posar text a xlab, ylab o main(title)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 37 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 37 rows containing missing values or values outside the scale range
## (`geom_point()`).

 # b) 
 ggplot(data = airquality, aes(x = Solar.R, y = Temp, col=factor(Month),shape=factor(Month))) + geom_point() + labs(title="Relació entre Solar.R i Temperatura", x = "Solar.R", y = "Temperatura (en Fahrenheit)")
## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_point()`).

 # c) 
 # Carreguem i activem les dades
 library("MASS")
## 
## Adjuntando el paquete: 'MASS'
## 
## The following object is masked from 'package:dplyr':
## 
##     select
 data("birthwt")
 
 ggplot(data = birthwt, aes(x = age)) + geom_histogram(col = "black", fill = "blue")
## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.

 # d)
 ggplot(data = birthwt, aes(x = age)) + geom_histogram(col = "black", fill = "blue") + facet_wrap(~smoke)
## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.

Exercici 6

Creeu un gràfic amb unes dades extretes del paquet Datasets d’R i deseu-lo com a imatge (.jpg) amb el nom migrafic1 i, també, com a document en PDF amb el nom migrafic2. Feu una captura de pantalla del fitxer generat.
data("faithful")

jpeg("migrafic1.jpg")
pdf("migrafic2.pdf")
plot(faithful, main = "Correlació entre els minuts d'erupció del gèiser i espera fins la següent", xlab = "Durada de l'erupció (en minuts)", ylab= "Espera (en minuts)", col = blues9)
dev.off()
## png 
##   2
jpeg("migrafic1.jpg")
plot(faithful, main = "Correlació entre els minuts d'erupció del gèiser i espera fins la següent", xlab = "Durada de l'erupció (en minuts)", ylab= "Espera (en minuts)", col = blues9)
dev.off()
## png 
##   2

Els documents s’han generat correctament, tal com es mostra en la captura de pantalla:

Captura
Captura

Exercici 7

Per practicar la regressió lineal simple usarem el conjunt de dades Orange que es troba a la llibreria tidyverse i que té informació sobre tres variables (arbre, edat en dies des que es va sembrar l’arbre i circumferència del tronc en centímetres) de 35 tarongers.
• Volem saber quin valor de circumferència tindrà un arbre 600 dies després de plantar-lo. Creeu per a això un model lineal i practiqueu la regressió lineal pas a pas amb els passos que heu vist en aquest laboratori.
# Pas 1. Carreguem i activem les dades
library(tidyverse)
data("Orange")

# Pas 2. Fem un resum estadístic
summary(Orange)
##  Tree       age         circumference  
##  3:7   Min.   : 118.0   Min.   : 30.0  
##  1:7   1st Qu.: 484.0   1st Qu.: 65.5  
##  5:7   Median :1004.0   Median :115.0  
##  2:7   Mean   : 922.1   Mean   :115.9  
##  4:7   3rd Qu.:1372.0   3rd Qu.:161.5  
##        Max.   :1582.0   Max.   :214.0
# Pas 3. Generem un diagrama de dispersió
pairs(Orange) # Observant els gràfics generats, sembla que hi ha una relació bastant clara entre l'edat i la circumferència. 

# Pas 4. Mirem el grau de relació entre les dues variables
kable(cor(Orange$age, Orange$circumference))
x
0.9135189 
# Hem obtingut una relació de 0.9135189, que al ser molt propera a 1, ens indica que molt segurament hi ha una relació entre les dues variables

# Pas 5. Generem el nou model
nou_model <- lm(circumference ~ age, data = Orange)
summary(nou_model)
## 
## Call:
## lm(formula = circumference ~ age, data = Orange)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -46.310 -14.946  -0.076  19.697  45.111 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 17.399650   8.622660   2.018   0.0518 .  
## age          0.106770   0.008277  12.900 1.93e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 23.74 on 33 degrees of freedom
## Multiple R-squared:  0.8345, Adjusted R-squared:  0.8295 
## F-statistic: 166.4 on 1 and 33 DF,  p-value: 1.931e-14
# A partir d'aquest Summary() hem obtingut l'equació de la recta, que serà:
# y = 17.399650 + 0.106770x

# Pas 6. Generem el gràfic de dispersió
plot(Orange$age, Orange$circumference, xlab = "Age", ylab = "Circumferència del tronc")
abline(nou_model)

# Pas 7. Fem la validació del model
residus <- rstandard(nou_model)
valors.ajustats <- fitted(nou_model)
plot(valors.ajustats, residus)

# En aquest cas si que podem observar una mena de patró en la validació, que ens pot indicar que seguir un model linear pot no ser la millor opció, tot i tenir un R^2 elevat. Per a que la validació fós efectiva, hauriem de veure un núvol de punts que semblen posats a l'atzar.

# Pas 8. Fem test de normalitat
qqnorm(residus)
qqline(residus)

# Pas 9. Finalment fem la predicció per a 600 dies
prediccio <- data.frame(age = seq(600))
resultats <- predict(nou_model, prediccio)
tail(resultats, 10) # Per tal d'evitar tenir la llista de 600 mesures, fem que només mostri els últims 10
##      591      592      593      594      595      596      597      598 
## 80.50091 80.60768 80.71445 80.82122 80.92799 81.03476 81.14153 81.24830 
##      599      600 
## 81.35507 81.46185
# La predicció als 600 dies és que faci 81.46185 mm

# Pas 10. Fem un interval de confiança pels resultats que ens ha donat
confint(nou_model)
##                   2.5 %     97.5 %
## (Intercept) -0.14328303 34.9425835
## age          0.08993141  0.1236092

Exercici 8

Repetiu els gràfics (aquells que pugueu) de la regressió lineal simple de l’exercici anterior però ara amb ggplot2.
# Carreguem i activem les dades
library(tidyverse)
data("Orange")

nou_model <- lm(circumference ~ age, data = Orange)

# Generem el gràfic de dispersió del Pas 6. 
ggplot(Orange, aes(x = age, y = circumference)) + geom_point(size = 2) + geom_smooth(method = "lm", color = "darkred") +  labs(title = "Regressió Lineal Ajustada",
       subtitle = "y = 17.40 + 0.1068x",
       x = "Edat (dies)", y = "Circumferència (mm)")
## `geom_smooth()` using formula = 'y ~ x'

# Per al gràfic de la validació del model hem de crear un dataframe. Pas 7. 
df_residus <- data.frame(
  residus <- rstandard(nou_model),
  valors.ajustats <- fitted(nou_model)
  )
ggplot(df_residus, aes(x= valors.ajustats, y= residus)) + geom_point(size = 2, col= "blue")

Exercici 9

Activeu el conjunt de dades PlantGrowth del paquet datasets d’R. Aquest arxiu té els resultats d’un experiment per comparar els rendiments mitjans pel pes sec de les plantes (weight) obtinguts sota un control i dues condicions de tractament diferents (group factor).
• Creieu que hi ha diferències entre tractaments?

Com hem pogut observar en el boxplot, hi ha diferències significatives entre els tractaments. Tenim que el grup que ha rebut el tractament 1, ha tingut un augment de pes sec inferior al grup control, mentre que el grup que ha rebut el tractament 2 ha tingut un augment de pes superior al grup control

• Es compleixen les condicions per poder aplicar una ANOVA. Quines proves us plantegeu?

Per aplicar ANOVA s’han de complir tres condicions:

- Independència. Se suposa que per com s’ha plantejat l’experiment, hi ha independència entre les mostres

- Normalitat. Analitzarem els residus del model lineal pel Test de Shapiro-Wilk

- Homoscedasticitat. Comprovem si les variàncies entre els grups són iguals pel Test de Bartlett

Si les tres condicions es compleixen, podrem aplicar ANOVA

library(ggplot2)

# Activem les dades
data("PlantGrowth")


# Fem un resum de les dades
summary(PlantGrowth) # Veiem que tenim tres grups, ctrl, trt1 i trt2. 
##      weight       group   
##  Min.   :3.590   ctrl:10  
##  1st Qu.:4.550   trt1:10  
##  Median :5.155   trt2:10  
##  Mean   :5.073            
##  3rd Qu.:5.530            
##  Max.   :6.310
# Fem un gràfic per establir les diferències entre els tractaments visualment
ggplot(PlantGrowth, aes(x = group, y = weight)) + geom_boxplot() + labs(title = "Distribució del pes per grup de tractament",
       x = "Grup", y = "Pes sec (weight)")

# Pel segon apartat, començarem determinant la normalitat, amb el test de Shapiro-Wilk
by(PlantGrowth$weight, PlantGrowth$group, shapiro.test)
## PlantGrowth$group: ctrl
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.95668, p-value = 0.7475
## 
## ------------------------------------------------------------ 
## PlantGrowth$group: trt1
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.93041, p-value = 0.4519
## 
## ------------------------------------------------------------ 
## PlantGrowth$group: trt2
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.94101, p-value = 0.5643
# I per la heterocedasticitat fem el test de Bartlett
bartlett.test(PlantGrowth$weight, PlantGrowth$group)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  PlantGrowth$weight and PlantGrowth$group
## Bartlett's K-squared = 2.8786, df = 2, p-value = 0.2371
# En tots els casos obtenim que p > 0.05, pel que no hauriem de tenir problemes d'heterocedasticitat ni de normalitat

# Finalment, apliquem ANOVA
anv <- aov(weight ~ group, data = PlantGrowth)
summary(anv)
##             Df Sum Sq Mean Sq F value Pr(>F)  
## group        2  3.766  1.8832   4.846 0.0159 *
## Residuals   27 10.492  0.3886                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Podem fer Tukey per determinar entre quins grups hi ha diferències ( p-value < 0.05 en ANOVA ens indica que hi ha diferències significatives entre almenys dos grups)
TukeyHSD(anv)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = weight ~ group, data = PlantGrowth)
## 
## $group
##             diff        lwr       upr     p adj
## trt1-ctrl -0.371 -1.0622161 0.3202161 0.3908711
## trt2-ctrl  0.494 -0.1972161 1.1852161 0.1979960
## trt2-trt1  0.865  0.1737839 1.5562161 0.0120064
# Veiem que la diferència principal és entre el tractament 1 i tractament 2 (p-value = 0.012)

Exercici 10

Busqueu informació del paquet Plotly per a la creació de gràfics interactius i genereu un histograma o un gràfic de barres interactiu. Expliqueu què es pot fer amb aquest gràfic.

El paquet Plotly permet crear gràfics on qualsevol usuari pot moure el ratolí per sobre dels punts per veure quins valors tenen, o pots fer zoom en els diferents elements del gràfic i filtrar dades

library(plotly)
## 
## Adjuntando el paquete: 'plotly'
## The following object is masked from 'package:MASS':
## 
##     select
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(tidyverse)
data("iris")

histograma <- plot_ly(data = iris, x= iris$Sepal.Length, type = "histogram")

histograma

En aquest gràfic creat per Plotly podem:

  • Etiquetes: Si posem el ratolí sobre una barra, apareix un requadre amb informació dels valors que representa a l’eix de les x, aixi com al de les y.

  • Zoom: Mantenint el clic esquerre i dibuixant un rectangle podem fer zoom sobre les dades, si estiguessin molt juntes.

  • Exportar: A la barra superior hi tenim un icone d’una càmera, que permet exportar el gràfic a PDF.

  • Box Select. Ens serveix per seleccionar subconjunts, ressaltar visualment alguns d’aquests o filtrar dinàmicament

CAS PRÀCTIC

Seleccionaré les dades “Plasma_Retinol” de la biblioteca StatLib datasets archive, Carnegie-Mellon University.

library(tidyverse)
library(readr)

# Importem les dades corresponents a Plasma_Retinol
data_retinol <- read_table("C:/Users/carle/Documents/r_files/Plasma_Retinol.txt", skip = 34, col_names = FALSE, col_types = cols(.default = "d")) # Importem el database sense les primeres 34 linies, que corresponen a text, i no a les dades com a tal. Faig col_types perque no em detectava les primeres 4 columnes com a números, i no em feia el resum estadístic, amb aquesta variable li indico que totes les columnes son double (numèriques)
## Warning: 19 parsing failures.
## row col   expected    actual                                                  file
## 312  X1 a double   Therese   'C:/Users/carle/Documents/r_files/Plasma_Retinol.txt'
## 312  X2 a double   Stukel    'C:/Users/carle/Documents/r_files/Plasma_Retinol.txt'
## 312  -- 14 columns 2 columns 'C:/Users/carle/Documents/r_files/Plasma_Retinol.txt'
## 313  X1 a double   Dartmouth 'C:/Users/carle/Documents/r_files/Plasma_Retinol.txt'
## 313  X2 a double   Hitchcock 'C:/Users/carle/Documents/r_files/Plasma_Retinol.txt'
## ... ... .......... ......... .....................................................
## See problems(...) for more details.
# Definim les diferents columnes i la seva categoria
colnames(data_retinol) <- c(
  "Edat", # AGE: Age (years)
    "Sexe", # SEX: Sex (1=Male, 2=Female).
    "Tabac", # SMOKSTAT: Smoking status (1=Never, 2=Former, 3=Current Smoker)
    "IMC", # QUETELET: Quetelet (weight/(height^2))
    "Vitamines", # VITUSE: Vitamin Use (1=Yes, fairly often, 2=Yes, not often, 3=No)
    "Calories", # CALORIES: Number of calories consumed per day.
    "Greix", # FAT: Grams of fat consumed per day.
    "Fibra", # FIBER: Grams of fiber consumed per day.
    "Alcohol", # ALCOHOL: Number of alcoholic drinks consumed per week.
    "Colesterol", # CHOLESTEROL: Cholesterol consumed (mg per day).
    "Beta_dieta", # BETADIET: Dietary beta-carotene consumed (mcg per day).
    "Retinol", # RETDIET: Dietary retinol consumed (mcg per day)
    "Beta_plasma", # BETAPLASMA: Plasma beta-carotene (ng/ml)
    "Retinol_plasma"# RETPLASMA: Plasma Retinol (ng/ml)
  )

# Netegem files buides
data_retinol <- na.omit(data_retinol)

# Fem el resum estadístic de tot el dataset
summary(data_retinol)
##       Edat            Sexe           Tabac            IMC       
##  Min.   :19.00   Min.   :1.000   Min.   :1.000   Min.   :16.33  
##  1st Qu.:39.00   1st Qu.:2.000   1st Qu.:1.000   1st Qu.:21.82  
##  Median :48.00   Median :2.000   Median :1.000   Median :24.74  
##  Mean   :50.09   Mean   :1.865   Mean   :1.637   Mean   :26.20  
##  3rd Qu.:62.00   3rd Qu.:2.000   3rd Qu.:2.000   3rd Qu.:28.98  
##  Max.   :83.00   Max.   :2.000   Max.   :3.000   Max.   :50.40  
##    Vitamines        Calories          Greix            Fibra      
##  Min.   :1.000   Min.   : 445.2   Min.   : 14.40   Min.   : 3.10  
##  1st Qu.:1.000   1st Qu.:1354.4   1st Qu.: 53.95   1st Qu.: 9.15  
##  Median :2.000   Median :1666.8   Median : 72.90   Median :12.10  
##  Mean   :1.968   Mean   :1796.8   Mean   : 77.10   Mean   :12.74  
##  3rd Qu.:3.000   3rd Qu.:2090.6   3rd Qu.: 95.25   3rd Qu.:15.35  
##  Max.   :3.000   Max.   :6662.2   Max.   :235.90   Max.   :36.80  
##     Alcohol          Colesterol      Beta_dieta      Retinol      
##  Min.   :  0.000   Min.   : 37.7   Min.   : 214   Min.   :  30.0  
##  1st Qu.:  0.000   1st Qu.:155.0   1st Qu.:1116   1st Qu.: 480.0  
##  Median :  0.300   Median :206.3   Median :1788   Median : 707.0  
##  Mean   :  3.275   Mean   :242.9   Mean   :2175   Mean   : 834.2  
##  3rd Qu.:  3.200   3rd Qu.:308.9   3rd Qu.:2836   3rd Qu.:1047.5  
##  Max.   :203.000   Max.   :900.7   Max.   :9642   Max.   :6901.0  
##   Beta_plasma     Retinol_plasma  
##  Min.   :   0.0   Min.   : 179.0  
##  1st Qu.:  89.0   1st Qu.: 466.0  
##  Median : 138.0   Median : 564.0  
##  Mean   : 189.7   Mean   : 601.0  
##  3rd Qu.: 230.0   3rd Qu.: 709.5  
##  Max.   :1415.0   Max.   :1727.0
# El resultat d'aquest ens indica, per exemple, que tenim més dones que homes en les mostres (Sexe té una mitjana més propera a 2), que més de la meitat de les persones no son fumadores (Median = 1, per tant, el 50% dels valors son 1), i a més ens indiquen el rang de valors per totes les diferents variables, a més de les seves mitjanes, medianes i quartils. 

# Procedim a fer 5 gràfics bàsics i els guardem en format imatge
jpeg("01_hist.jpg")
hist(data_retinol$Edat)
dev.off()
## png 
##   2
jpeg("02_boxplot.jpg")
boxplot(data_retinol$Sexe, data_retinol$Greix)
dev.off()
## png 
##   2
jpeg("03_plot.jpg")
plot(data_retinol$Beta_dieta, data_retinol$Retinol_plasma)
dev.off()
## png 
##   2
jpeg("04_barplot.jpg")
barplot(table(data_retinol$Tabac))
dev.off()
## png 
##   2
jpeg("05_hist.jpg")
hist(data_retinol$Fibra)
dev.off()
## png 
##   2
# Generem dos gràfics amb ggplot2:
# Analitzem els nivells d'ingesta de greix segons l'edat
jpeg("05_ggplot1.jpg")
ggplot(data_retinol, aes(x = Edat, y = Greix)) + geom_point(size = 2) + geom_smooth(method = "lm", color = "darkred") +  labs(title = "Regressió Lineal Ajustada",
       x = "Edat (anys)", y = "Greix consumit per dia")
## `geom_smooth()` using formula = 'y ~ x'
dev.off()
## png 
##   2
# Analitzem els nivells de retinol en plasma segons els tres tipus de suplementació amb vitamines. 
jpeg("06_ggplot2.jpg")
ggplot(data_retinol, aes(x = factor(Vitamines), y = Retinol_plasma, fill = factor(Vitamines))) + geom_violin(alpha=0.7) + labs(title = "Nivells de Retinol segons Vitamines", x = "Ús de Vitamines", y = "Retinol (ng/ml)")
dev.off()
## png 
##   2
# Finalment, fem una regressió lineal entre les variables Edat i Retinol en plasma
kable(cor(data_retinol$Retinol_plasma, data_retinol$Edat))
x
0.2080947 
modelreg <- lm(Retinol_plasma~Edat, data = data_retinol)
summary(modelreg)
## 
## Call:
## lm(formula = Retinol_plasma ~ Edat, data = data_retinol)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -458.60 -126.77  -32.79  101.70 1105.35 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 451.0212    41.7456   10.80  < 2e-16 ***
## Edat          2.9935     0.8004    3.74 0.000219 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 205 on 309 degrees of freedom
## Multiple R-squared:  0.0433, Adjusted R-squared:  0.04021 
## F-statistic: 13.99 on 1 and 309 DF,  p-value: 0.0002194
plot(data_retinol$Edat, data_retinol$Retinol_plasma, xlab = "Edat", ylab = "Nivells de retinol en plasma")
abline(modelreg)

residus <- rstandard(modelreg)
valors.ajustats <- fitted(modelreg)
plot(valors.ajustats, residus)

qqnorm(residus)
qqline(residus)

# Volem predir quins seran els nivells de Retinol en plasma en una persona de 75-80 anys
pred_edat <- data.frame(Edat=seq(75, 80))
predict(modelreg, pred_edat)
##        1        2        3        4        5        6 
## 675.5369 678.5304 681.5240 684.5175 687.5111 690.5046
# Calculem l'interval de confiaça al 95%
confint(modelreg)
##                  2.5 %     97.5 %
## (Intercept) 368.879578 533.162783
## Edat          1.418527   4.568559
# Amb l'anàlisi dut a terme, veient que R^2 (0.04) i la correlació (0.2) son tant baixes, podiem determinar que no hi ha relació (o aquesta és insuficient) entre les dues variables.

Confirmació de que els 5 gràfics bàsics s’han generat i guardat correctament: Captura