Anàlisi de dades amb R
Malu Calle
1 Introdució
En aquest tutorial revisarem els principals mètodes d’anàlisi de dades amb R. Per il·lustrar els procediments utilitzarem el conjunt de dades, “nadons.xlsx”, que podeu trobar aquí: https://malucalle.github.io/statistical-pills/nadons.xlsx
2 Preparar les dades
2.1 Estructurar bé les dades
Abans de començar una anàlisi de dades cal plantejar-se si les dades estan ben estructurades:
Quins són els “individus” o “casos” de l’estudi?
Quina informació recollim de cada individu? variables
Una fila no pot contenir informació de més d’un individu
Les dades han de tenir una estructura quadrada:
Files = individus = casos Columnes = variables
2.2 Especificar la carpeta de treball
Especificar la carpeta de treball directament des de RStudio:
Session -> Set working directory -> choose directory2.3 Importar les dades a R
Importar les dades a R directament des de RStudio
File -> Import Dataset
-> From Excel... -> Browse (carregueu el fitxer "nadons.xlsx")
Comprovar que les variables numèriques s'importen com a numèriques o double També es pot fer amb codi d’R amb la llibreria d’R readxl:
Instal·lar la llibreria. Només cal fer-ho una vegada. El nom ha d’anar entre cometes:
install.packages("readxl")Activar la llibreria sempre que es comença una nova sessió d’R. El nom no cal que estigui entre cometes:
library(readxl)Utilitzem la funció read_excel per carregar les dades:
nadons <- read_excel("nadons.xlsx")Mirar la taula de dades que hem carregat:
Clicar sobre el nom de la taula a la finestra “Environment” i comprovar que les dades s’han carregat bé. També es pot fer amb codi d’R:
View(nadons)Canviarem el nom de la taula a “data”:
data<-nadonsI ara eliminem (rm=remove) la taula “nadons”:
rm(nadons) Mirar l’estructura de les dades i el tipus de cada variable:
str(data)## tibble [200 x 11] (S3: tbl_df/tbl/data.frame)
## $ id : chr [1:200] "001" "002" "003" "004" ...
## $ nado_sexe : chr [1:200] "nen" "nena" "nen" "nena" ...
## $ nado_pes : num [1:200] 3.83 3.62 4.15 2.67 3.76 3.77 3.11 3.49 2.37 3.13 ...
## $ nado_baix_pes : num [1:200] 0 0 0 0 0 0 0 0 1 0 ...
## $ setmanes_gestacio: num [1:200] 40 38 41 39 39 40 40 37 40 39 ...
## $ prematur : chr [1:200] "no" "no" "no" "no" ...
## $ edat_mare : num [1:200] 1 2 2 2 2 1 2 1 2 2 ...
## $ nombre_visites : num [1:200] 2 4 2 2 3 3 3 4 10 1 ...
## $ mare_fumadora : chr [1:200] "no" "no" "no" "si" ...
## $ mare_pes_ini : num [1:200] 62 56 66 63 76 62 63 70 70 63 ...
## $ mare_pes_fi : num [1:200] 74 73 80 76 95 79 81 78 81 76 ...Veiem que totes les variables s’han importat com a numèriques (num) o caràcter (chr)
2.4 Resum numèric conjunt de totes les variables
summary(data)## id nado_sexe nado_pes nado_baix_pes
## Length:200 Length:200 Min. :1.510 Min. :0.000
## Class :character Class :character 1st Qu.:2.748 1st Qu.:0.000
## Mode :character Mode :character Median :3.230 Median :0.000
## Mean :3.201 Mean :0.185
## 3rd Qu.:3.612 3rd Qu.:0.000
## Max. :5.250 Max. :1.000
## setmanes_gestacio prematur edat_mare nombre_visites
## Min. :22.00 Length:200 Min. :1.000 Min. : 0.00
## 1st Qu.:36.00 Class :character 1st Qu.:1.000 1st Qu.: 2.00
## Median :38.00 Mode :character Median :2.000 Median : 3.00
## Mean :37.37 Mean :1.825 Mean : 3.72
## 3rd Qu.:40.00 3rd Qu.:2.000 3rd Qu.: 4.00
## Max. :44.00 Max. :3.000 Max. :20.00
## mare_fumadora mare_pes_ini mare_pes_fi
## Length:200 Min. :50 Min. :63.00
## Class :character 1st Qu.:62 1st Qu.:75.00
## Mode :character Median :65 Median :79.00
## Mean :65 Mean :78.64
## 3rd Qu.:69 3rd Qu.:82.00
## Max. :77 Max. :95.00Com es pot veure el resum de les dades categòriques no es adequat i cal preparar-les per l’anàlisi
2.5 Preparar les dades categòriques
Amb la funció factor() preparem les dades categòriques per l’anàlisi:
Variables caràcter (chr):
data$prematur<-factor(data$prematur)
data$nado_sexe<-factor(data$nado_sexe)
data$mare_fumadora<-factor(data$mare_fumadora)Variables numèriques (num) que en realitat són categòriques:
Cal especificar què representa cada valor:
data$edat_mare<-factor(data$edat_mare, levels = c(1,2,3), labels = c("<=30", "30-40", ">40"))
data$nado_baix_pes<-factor(data$nado_baix_pes, levels = c(0,1), labels = c("no", "si"))Tornem a fer un resum numèric conjunt de totes les variables:
summary(data)## id nado_sexe nado_pes nado_baix_pes setmanes_gestacio
## Length:200 nen : 94 Min. :1.510 no:163 Min. :22.00
## Class :character nena:106 1st Qu.:2.748 si: 37 1st Qu.:36.00
## Mode :character Median :3.230 Median :38.00
## Mean :3.201 Mean :37.37
## 3rd Qu.:3.612 3rd Qu.:40.00
## Max. :5.250 Max. :44.00
## prematur edat_mare nombre_visites mare_fumadora mare_pes_ini
## no:141 <=30 : 65 Min. : 0.00 no:160 Min. :50
## si: 59 30-40:105 1st Qu.: 2.00 si: 40 1st Qu.:62
## >40 : 30 Median : 3.00 Median :65
## Mean : 3.72 Mean :65
## 3rd Qu.: 4.00 3rd Qu.:69
## Max. :20.00 Max. :77
## mare_pes_fi
## Min. :63.00
## 1st Qu.:75.00
## Median :79.00
## Mean :78.64
## 3rd Qu.:82.00
## Max. :95.002.6 Gràfics amb R
Utilitzarem la llibreria ggpubr que permet fer gràfics més interessants que els que es poden fer amb les funcions bàsiques d’R.
Podeu veure exemples dels gràfics més importants en aquest tutorial: https://malucalle.github.io/statistical-pills/nice-plots-with-r.html
#install.packages("ggpubr")
library(ggpubr) # activar la llibreria. S'ha de fer cada cop que s'engega R2.7 Histogrames de les variables numèriques
Exemple: Setmanes de gestació:
gghistogram(data, x="setmanes_gestacio",
bins = 20,
color="aquamarine2", fill="aquamarine2")
Podeu representar totes les variables numèriques de cop:
p<-gghistogram(data,
x=c("setmanes_gestacio","nombre_visites","nado_pes","mare_pes_ini","mare_pes_fi"), bins = 20,
color = "coral1", fill = "coral1")
p## $setmanes_gestacio
##
## $nombre_visites
##
## $nado_pes
##
## $mare_pes_ini
##
## $mare_pes_fi
2.8 Boxplots de les variables numèriques
Exemple: Setmanes de gestació:
p<-ggboxplot(data, y="setmanes_gestacio", color = "blue")
p
2.9 Diagrames de barres de les variables categòriques
Exemple: Edat de la mare:
Per poder representar les variables categòriques primer hem de fer la taula de freqüències:
table(data$edat_mare)##
## <=30 30-40 >40
## 65 105 30i la taula de freqüències relatives:
prop.table(table(data$edat_mare))##
## <=30 30-40 >40
## 0.325 0.525 0.150Gràfic de barres de la taula de freqüències:
barplot(table(data$edat_mare))
Gràfic de barres de les freqüències relatives
barplot(prop.table(table(data$edat_mare)),
col=c("aquamarine2", "coral1","darkturquoise"))
2.10 Variables comptatge (recompte)
Mirar si alguna de les variables són comptatges (recomptes) perquè requeriran una atenció especial.
Les variables comptatges es poden analitzar com si fossin numèriques o com categòriques i mirar quin dels dos anàlisis ens aporta més informació sobre la variable:
La variable “nombre_visites” és un recompte. La podem analitzar segons una variable numèrica:
summary(data$nombre_visites)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 2.00 3.00 3.72 4.00 20.00gghistogram(data, x="nombre_visites",
bins = 20,
color = "coral1", fill = "coral1")
o com si fos categòrica:
table(data$nombre_visites)##
## 0 1 2 3 4 5 8 9 10 11 15 16 17 20
## 4 14 36 52 70 10 1 2 4 2 1 2 1 1prop.table(table(data$nombre_visites))##
## 0 1 2 3 4 5 8 9 10 11 15 16 17
## 0.020 0.070 0.180 0.260 0.350 0.050 0.005 0.010 0.020 0.010 0.005 0.010 0.005
## 20
## 0.005barplot(table(data$nombre_visites))
i també la agrupar en noves categories:
Podem subdividir el nombre de visites en 2 categories: “<= 4” i “superiors a 4”
data$visites_cat2<-cut(data$nombre_visites, c(0,4,20))
table(data$visites_cat2)##
## (0,4] (4,20]
## 172 24prop.table(table(data$visites_cat2))##
## (0,4] (4,20]
## 0.877551 0.122449barplot(table(data$visites_cat2),
col=c("aquamarine2", "coral1"))
També podem subdividir el nombre de visites en 3 categories: “<=2”, “>2 i <=4” i “>4”:
data$visites_cat3<-cut(data$nombre_visites, c(0,2,4,20))
table(data$visites_cat3)##
## (0,2] (2,4] (4,20]
## 50 122 24prop.table(table(data$visites_cat3))##
## (0,2] (2,4] (4,20]
## 0.255102 0.622449 0.122449barplot(table(data$visites_cat3),
col=c("aquamarine2", "coral1","darkturquoise"))
2.11 Informe descripció inicial de les dades
Resum dels resultats obtinguts fins al moment:
És molt important tenir una bona descripció inicial de les dades que analitzarem. És recomanable fer un informe d’aquesta descripció inicial, tot i que aquesta part segurament no caldrà incloure-la a l’informe final (es pot afegir com a annex).
En aquesta fase es poden detectar possibles valors anòmals i decidir què fer amb ells. En general, eliminar valors no és una opció acceptable, llevat que es tracti de valors extremadament diferents de la majoria i que distorsionen molt els resultats. Si es decideix eliminar algun valor cal fer-ho de manera transparent i, en tot cas, informar d’aquest fet.
3 Analisi de les dades
Concretar i prioritzar els objectius:
Definir l’objectiu principal de l’estudi
- Identificar els factors associats al baix pes dels nens
Identificar la variable principal de l’estudi (outcome)
- Variables principals: “nado_pes”, “nado_baix_pes”
Definir objectius que involucrin dues variables i mirar de quin tipus són:
Quin és l’efecte de l’edat de la mare sobre el pes del nadó?
"edat_mare" (categòrica) i "nado_pes" (numèrica) "edat_mare" (categòrica) i "nado_baix_pes" (categòrica)Hi ha relació entre el fet que la mare fumi i el pes del nadó?
"mare_fumadora" (categòrica) i "nado_pes" (numèrica) "mare_fumadora" (categòrica) i "nado_baix_pes" (categòrica)Quin és l’efecte del nombre de setmanes de gestació sobre el pes del nadó?
"setmanes_gestacio" (numèrica) i "nado_pes" (numèrica) "setmanes_gestacio" (numèrica) i "nado_baix_pes" (categòrica)
Definir objectius que involucrin més de dues variables.
Quin és l’efecte de l’edat de la mare sobre el pes del nadó en funció de si fuma o no fuma?
"edat_mare" (categòrica), "mare_fumadora" (categòrica) i "nado_pes" (numèrica) "edat_mare" (categòrica), "mare_fumadora" (categòrica) i "nado_baix_pes" (categòrica)
4 Anàlisi relació dues var. numèriques
Exemple: Analitzar la possible relació entre el nombre de setmanes de gestació i el pes dels nadons
“setmanes_gestacio” (numèrica) i “nado_pes” (numèrica)
Fem un diagrama de punts (scatter plot):
ggscatter(data, x = "setmanes_gestacio", y = "nado_pes",
add = "reg.line", # regression line
conf.int = TRUE, # confidence interval
add.params = list(color = "coral1", fill = "bisque") # color
)+
stat_cor(method = "pearson") # correlation coefficient
5 Anàlisi relació var. numèrica i var.categòria
5.1 Histograma de la var. numèrica en funció de la categòrica
Exemple: Analitzar la possible relació entre el nombre de setmanes de gestació i que el nadó tingui baix pes o no
“setmanes_gestacio” (numèrica) i “nado_baix_pes” (categòrica)
Fem l’histograma posant la variable numèrica a la “x” i la categòrica a “color”:
gghistogram(data,
x="setmanes_gestacio",
color="nado_baix_pes",
add = "mean",
bins = 20)
5.2 Boxplot de la var. numèrica en funció de la categòrica
Exemple: Analitzar la possible relació entre el nombre de setmanes de gestació i que el nadó tingui baix pes o no
“setmanes_gestacio” (numèrica) i “nado_baix_pes” (categòrica)
Fem el gràfic de caixes múltiple posant la variable numèrica a la “y” i la categòrica a la “x”:
boxp<-ggboxplot(data,
y="setmanes_gestacio",
x="nado_baix_pes",
color = "nado_baix_pes",
)
boxp
5.3 Resum de la variable numèrica en funció de les categories
Exemple: Analitzar la possible relació entre el nombre de setmanes de gestació i que el nadó tingui baix pes o no
“setmanes_gestacio” (numèrica) i “nado_baix_pes” (categòrica)
El resum numèric es pot fer de dues maneres, amb la funció tapply():
tapply(data$setmanes_gestacio, data$nado_baix_pes, summary)## $no
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 22.0 37.0 39.0 38.3 40.0 44.0
##
## $si
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 30.00 31.00 33.00 33.24 34.00 40.00o amb la funció desc_statby() de la llibreria ggpubr:
desc_statby(data,"setmanes_gestacio", "nado_baix_pes") # amb la funció desc_statby de ggpubr ## nado_baix_pes length min max median mean iqr mad sd se
## 1 no 163 22 44 39 38.30061 3 1.4826 3.233954 0.2533028
## 2 si 37 30 40 33 33.24324 3 2.9652 2.361796 0.3882768
## ci range cv var
## 1 0.5002010 22 0.08443608 10.458456
## 2 0.7874618 10 0.07104588 5.5780786 Proves d’hipòtesis d’igualtat de mitjanes
Si les dades en cada categoria són normals (p-valor shapiro test >= alpha) apliquem:
t.test (2 categories)
anova (> 2 categories)
Si les dades en cada categoria no són normals (p-valor shapiro test < alpha) apliquem:
wilcoxon test (2 categories)
Kruskal-wallis (> 2 categories)
Exemple1: Analitzar la possible relació entre el nombre de setmanes de gestació i que el nadó tingui baix pes o no
“setmanes_gestacio” (numèrica) i “nado_baix_pes” (categòrica)
Apliquem el test de normalitat:
tapply(data$setmanes_gestacio, data$nado_baix_pes, function(x) shapiro.test(x))## $no
##
## Shapiro-Wilk normality test
##
## data: x
## W = 0.8354, p-value = 2.893e-12
##
##
## $si
##
## Shapiro-Wilk normality test
##
## data: x
## W = 0.92247, p-value = 0.01319Com que les dades no segueixen una distribució normal apliquem el test de Wilcoxon d’igualtat de dues mitjanes:
wilcox.test(data$setmanes_gestacio~data$nado_baix_pes)##
## Wilcoxon rank sum test with continuity correction
##
## data: data$setmanes_gestacio by data$nado_baix_pes
## W = 5457.5, p-value = 1.066e-14
## alternative hypothesis: true location shift is not equal to 0Exemple2: Analitzar la possible relació entre el pes inicial de la mare i que el nadó tingui baix pes o no
“mare_pes_ini” (numèrica) i “nado_baix_pes” (categòrica):
Apliquem el test de normalitat:
tapply(data$mare_pes_ini, data$nado_baix_pes, function(x) shapiro.test(x))## $no
##
## Shapiro-Wilk normality test
##
## data: x
## W = 0.98949, p-value = 0.2677
##
##
## $si
##
## Shapiro-Wilk normality test
##
## data: x
## W = 0.97341, p-value = 0.5082Com que les dades són normals, podrem aplicar el t-test per comparar les mitjanes. Però prèviament necessitem verificar si les variàncies són iguals o no:
var.test(data$mare_pes_ini~data$nado_baix_pes)##
## F test to compare two variances
##
## data: data$mare_pes_ini by data$nado_baix_pes
## F = 1.6908, num df = 162, denom df = 36, p-value = 0.06545
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.9658489 2.7172032
## sample estimates:
## ratio of variances
## 1.690776Apliquem el t-test amb variàncies iguals:
t.test(data$mare_pes_ini~data$nado_baix_pes, var.equal=T)##
## Two Sample t-test
##
## data: data$mare_pes_ini by data$nado_baix_pes
## t = -1.5509, df = 198, p-value = 0.1225
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -3.1777358 0.3798582
## sample estimates:
## mean in group no mean in group si
## 64.73620 66.13514Exemple 3: Analitzar la possible relació entre l’edat de la mare i el pes del nadó
“edat_mare” (categòrica) i “nado_pes” (numèrica)
Comprovem si les dades són normals:
tapply(data$nado_pes, data$edat_mare, function(x) shapiro.test(x))## $`<=30`
##
## Shapiro-Wilk normality test
##
## data: x
## W = 0.96744, p-value = 0.08446
##
##
## $`30-40`
##
## Shapiro-Wilk normality test
##
## data: x
## W = 0.9862, p-value = 0.3518
##
##
## $`>40`
##
## Shapiro-Wilk normality test
##
## data: x
## W = 0.95714, p-value = 0.2614Abans d’aplicar l’ANOVA cal comprovar si les variàncies dels diferents grups es poden considerar iguals:
#install.packages("lmtest")
library(lmtest)
bptest(lm(data$nado_pes ~data$edat_mare),studentize = FALSE)##
## Breusch-Pagan test
##
## data: lm(data$nado_pes ~ data$edat_mare)
## BP = 0.90261, df = 2, p-value = 0.6368Podem aplicar l’ANOVA per comparar les mitjanes de pes dels 3 grups d’edat:
summary(aov(data$nado_pes~data$edat_mare))## Df Sum Sq Mean Sq F value Pr(>F)
## data$edat_mare 2 0.17 0.0848 0.188 0.828
## Residuals 197 88.67 0.4501La conclusió és que no hi ha diferències. Si n’hi haguéssin aplicaríem un test per determinar quins grups tenen mitjanes diferents:
TukeyHSD(aov(data$nado_pes~data$edat_mare))## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = data$nado_pes ~ data$edat_mare)
##
## $`data$edat_mare`
## diff lwr upr p adj
## 30-40-<=30 0.006241758 -0.2438059 0.2562894 0.9980859
## >40-<=30 0.085051282 -0.2646475 0.4347501 0.8339475
## >40-30-40 0.078809524 -0.2491809 0.4068000 0.8375801Exemple 4: Analitzar la possible relació entre l’edat de la mare i el nombre de setmanes de gestació
“edat_mare” (categòrica) i “setmanes_gestacio” (numèrica)
Comprovem si les dades són normals:
tapply(data$setmanes_gestacio, data$edat_mare, function(x) shapiro.test(x))## $`<=30`
##
## Shapiro-Wilk normality test
##
## data: x
## W = 0.89778, p-value = 5.835e-05
##
##
## $`30-40`
##
## Shapiro-Wilk normality test
##
## data: x
## W = 0.92059, p-value = 9.623e-06
##
##
## $`>40`
##
## Shapiro-Wilk normality test
##
## data: x
## W = 0.90979, p-value = 0.0147Com que no es verifica la condició de normalitat apliquem el test de Kruskal-Wallis per comparar les mitjanes dels 3 grups:
kruskal.test(data$setmanes_gestacio~data$edat_mare)##
## Kruskal-Wallis rank sum test
##
## data: data$setmanes_gestacio by data$edat_mare
## Kruskal-Wallis chi-squared = 0.97933, df = 2, p-value = 0.61287 Anàlisi relació dues var. categòriques
Exemple: Analitzar la possible relació entre el fet que la mare fumi i que el nadó tingui baix pes o no
“mare_fumadora” (categòrica) i “nado_baix_pes” (categòrica)
Comencem fent la taula de freqüències i la taula de freqüències relatives per files:
freq<-table(data$mare_fumadora, data$nado_baix_pes)
freqbyrow<-prop.table(table(data$mare_fumadora, data$nado_baix_pes),1)Convertim la taula de freqüències per files en un data frame i posem bé el nom de les variables:
freqbyrow<-as.data.frame(freqbyrow) # transform the table into a data frame
freqbyrow## Var1 Var2 Freq
## 1 no no 0.89375
## 2 si no 0.50000
## 3 no si 0.10625
## 4 si si 0.50000colnames(freqbyrow)<-c("mare_fumadora","nado_baix_pes","Freq")
freqbyrow## mare_fumadora nado_baix_pes Freq
## 1 no no 0.89375
## 2 si no 0.50000
## 3 no si 0.10625
## 4 si si 0.50000Fem el gràfic de barres de la taula de freqüències per files:
ggbarplot(freqbyrow, x="mare_fumadora", y="Freq", # x = variable files de la taula
color = "nado_baix_pes",
label = TRUE,
lab.nb.digits = 2) # nombre de dígits de les proporcions
Apliquem la prova xi-quadrat d’independència entre dos factors:
chisq.test(freq)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: freq
## X-squared = 30.345, df = 1, p-value = 3.616e-08Quan la taula de freqüències és 2x2 podem aplicar el test de Fisher:
fisher.test(freq)##
## Fisher's Exact Test for Count Data
##
## data: freq
## p-value = 1.864e-07
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 3.49944 20.11074
## sample estimates:
## odds ratio
## 8.280577De vegades pot ser interessant mesurar el grau d’associació entre dos variables categòriques. Una possible mesura d’associació és la V de Cramer:
library(rcompanion)
cramerV(data$mare_fumadora, data$nado_baix_pes, bias.correct = T)## Cramer V
## 0.40048 Models de regressió
Els models de regressió són útils per estudiar la relació d’una o més variables sobre una variable d’interès (variable resposta).
8.1 Regressió lineal
El model de regressió lineal s’utilitza quan la variable resposta és numèrica (contínua):
Exemple: Quin és l’efecte de l’edat i el pes inicial de la mare sobre el pes del nadó en funció de si fuma o no fuma?
y = “nado_pes” (numèrica)
x1 = “edat_mare” (categòrica), x2 = “mare_fumadora” (categòrica), x3 = “mare_pes_ini”
model1<-lm(nado_pes~mare_fumadora+edat_mare+mare_pes_ini, data = data)
summary(model1)##
## Call:
## lm(formula = nado_pes ~ mare_fumadora + edat_mare + mare_pes_ini,
## data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.75729 -0.41119 -0.03987 0.42235 1.91645
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.783093 0.587203 6.443 8.98e-10 ***
## mare_fumadorasi -0.587786 0.111838 -5.256 3.85e-07 ***
## edat_mare30-40 -0.007793 0.100558 -0.078 0.938
## edat_mare>40 0.118196 0.139082 0.850 0.396
## mare_pes_ini -0.007362 0.009093 -0.810 0.419
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6294 on 195 degrees of freedom
## Multiple R-squared: 0.1304, Adjusted R-squared: 0.1126
## F-statistic: 7.313 on 4 and 195 DF, p-value: 1.655e-05Cal comprovar que es compleixen els requisits de la regressió lineal:
- Normalitat dels residus:
shapiro.test(residuals(model1))##
## Shapiro-Wilk normality test
##
## data: residuals(model1)
## W = 0.99039, p-value = 0.2037- Els residus estan centrats al voltant del zero:
plot(residuals(model1))
abline(a=0, b=0)
8.2 Regressió logística
El model de regressió logística s’utilitza quan la variable resposta és binària (categòrica, amb dues categories):
Exemple: Quin és l’efecte de l’edat de la mare en el risc que el nen tingui baix pes en funció de si la mare fuma o no fuma?
y = “nado_baix_pes” (binària)
x1=“edat_mare” (categòrica), x2=“mare_fumadora” (categòrica)
model2<-glm(nado_baix_pes~mare_fumadora+edat_mare, data = data, family = "binomial")
summary(model2)##
## Call:
## glm(formula = nado_baix_pes ~ mare_fumadora + edat_mare, family = "binomial",
## data = data)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.2062 -0.4861 -0.4614 -0.4540 2.1558
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.18641 0.39274 -5.567 2.59e-08 ***
## mare_fumadorasi 2.14358 0.41060 5.221 1.78e-07 ***
## edat_mare30-40 0.11039 0.44787 0.246 0.805
## edat_mare>40 -0.03424 0.61464 -0.056 0.956
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 191.56 on 199 degrees of freedom
## Residual deviance: 163.71 on 196 degrees of freedom
## AIC: 171.71
##
## Number of Fisher Scoring iterations: 49 Selecció de variables: Step-wise regression
La selecció de variables s’utilitza per identificar quines variables estan més relacionades amb la variable d’interès.
Exemple: Quines variables estan relacionades amb el pes del nadó?
L’objectiu és identificar quines variables estan realment relacionades amb el pes del nadó i quines es poden eliminar del model. Afegirem en l’anàlisi l’augment de pes de la mare durant l’embaràç:
data$augment_pes<-data$mare_pes_fi-data$mare_pes_ini
model3<-lm(nado_pes~., data = data[, c(2,3,7,8,9,14)])
summary(model3)##
## Call:
## lm(formula = nado_pes ~ ., data = data[, c(2, 3, 7, 8, 9, 14)])
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.68535 -0.39363 -0.04404 0.37876 2.00964
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.437544 0.202720 12.024 < 2e-16 ***
## nado_sexenena -0.174646 0.084276 -2.072 0.0396 *
## edat_mare30-40 0.049287 0.093836 0.525 0.6000
## edat_mare>40 0.106804 0.129932 0.822 0.4121
## nombre_visites 0.007366 0.015220 0.484 0.6289
## mare_fumadorasi -0.508147 0.105874 -4.800 3.18e-06 ***
## augment_pes 0.065103 0.012075 5.391 2.03e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5877 on 193 degrees of freedom
## Multiple R-squared: 0.2497, Adjusted R-squared: 0.2264
## F-statistic: 10.71 on 6 and 193 DF, p-value: 2.889e-10Apliquem el mètode de selecció de variables (step-wise regression):
step(model3)## Start: AIC=-205.77
## nado_pes ~ nado_sexe + edat_mare + nombre_visites + mare_fumadora +
## augment_pes
##
## Df Sum of Sq RSS AIC
## - edat_mare 2 0.2457 66.896 -209.03
## - nombre_visites 1 0.0809 66.731 -207.53
## <none> 66.651 -205.77
## - nado_sexe 1 1.4830 68.134 -203.37
## - mare_fumadora 1 7.9551 74.606 -185.22
## - augment_pes 1 10.0384 76.689 -179.71
##
## Step: AIC=-209.03
## nado_pes ~ nado_sexe + nombre_visites + mare_fumadora + augment_pes
##
## Df Sum of Sq RSS AIC
## - nombre_visites 1 0.0866 66.983 -210.78
## <none> 66.896 -209.03
## - nado_sexe 1 1.5094 68.406 -206.57
## - mare_fumadora 1 7.9952 74.891 -188.46
## - augment_pes 1 10.1911 77.087 -182.68
##
## Step: AIC=-210.78
## nado_pes ~ nado_sexe + mare_fumadora + augment_pes
##
## Df Sum of Sq RSS AIC
## <none> 66.983 -210.78
## - nado_sexe 1 1.4496 68.432 -208.49
## - mare_fumadora 1 7.9650 74.948 -190.31
## - augment_pes 1 10.1165 77.099 -184.65##
## Call:
## lm(formula = nado_pes ~ nado_sexe + mare_fumadora + augment_pes,
## data = data[, c(2, 3, 7, 8, 9, 14)])
##
## Coefficients:
## (Intercept) nado_sexenena mare_fumadorasi augment_pes
## 2.51355 -0.17157 -0.50559 0.06445Veiem que el procediment de selecció de variables ha seleccionat el sexe del nadó, si la mare fuma i l’augment de pes de la mare com a variables rellevants en relació al pes del nadó. Fem un resum d’aquest model:
summary(step(model3))## Start: AIC=-205.77
## nado_pes ~ nado_sexe + edat_mare + nombre_visites + mare_fumadora +
## augment_pes
##
## Df Sum of Sq RSS AIC
## - edat_mare 2 0.2457 66.896 -209.03
## - nombre_visites 1 0.0809 66.731 -207.53
## <none> 66.651 -205.77
## - nado_sexe 1 1.4830 68.134 -203.37
## - mare_fumadora 1 7.9551 74.606 -185.22
## - augment_pes 1 10.0384 76.689 -179.71
##
## Step: AIC=-209.03
## nado_pes ~ nado_sexe + nombre_visites + mare_fumadora + augment_pes
##
## Df Sum of Sq RSS AIC
## - nombre_visites 1 0.0866 66.983 -210.78
## <none> 66.896 -209.03
## - nado_sexe 1 1.5094 68.406 -206.57
## - mare_fumadora 1 7.9952 74.891 -188.46
## - augment_pes 1 10.1911 77.087 -182.68
##
## Step: AIC=-210.78
## nado_pes ~ nado_sexe + mare_fumadora + augment_pes
##
## Df Sum of Sq RSS AIC
## <none> 66.983 -210.78
## - nado_sexe 1 1.4496 68.432 -208.49
## - mare_fumadora 1 7.9650 74.948 -190.31
## - augment_pes 1 10.1165 77.099 -184.65##
## Call:
## lm(formula = nado_pes ~ nado_sexe + mare_fumadora + augment_pes,
## data = data[, c(2, 3, 7, 8, 9, 14)])
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.66988 -0.39415 -0.06879 0.39318 2.05793
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.51355 0.17321 14.511 < 2e-16 ***
## nado_sexenena -0.17157 0.08330 -2.060 0.0408 *
## mare_fumadorasi -0.50559 0.10473 -4.828 2.78e-06 ***
## augment_pes 0.06445 0.01185 5.441 1.57e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5846 on 196 degrees of freedom
## Multiple R-squared: 0.246, Adjusted R-squared: 0.2345
## F-statistic: 21.32 on 3 and 196 DF, p-value: 5.417e-12