¡Descarga Entregable 4 GLM rcloud rstudio y más Ejercicios en PDF de Análisis de Datos Avanzado solo en Docsity!
ENTREGABLE_
Estela Torres Serrano
Cal que feu una anàlisi per veure si el fet de trobar o no una espècie invasora es pot relacionar amb els diferents blocs d’un índex que mesura l’estat de conservació d’una llacuna.
En total hi ha 5 blocs: bloc1 (fa referència a la morfologia de la llacuna); bloc2 (activitat humana, valors alts indiquen poca activitat); bloc3 (aspecte de l’aigua); bloc4 (vegetació emergent); bloc5 (vegetació aquàtica). Quan més alt és el valor més ben conservat està aquell aspecte a la llacuna. Les dades les trobareu a l’arxiu: (entregable_4.xlsx).
heu de tenir en compte que si no podeu fer una simplificació automàtica del model, i tot i així cal simplificar per arribar al model òptim, l’haureu de fer manualment. Ho podeu fer, si és necessari, seguint el procediment descrit per la regressió multiple a les sessions d’R"
Especifiqueu en les seccions corresponents les comandes utilitzades en cada cas, afegiu els comentaris pertinents i responeu les preguntes quan calgui.
- Quina variable creus que és la variable resposta i quin tipus de distribució té?
RESPOSTA: En aquest cas la variable resposta és la “sp_invasora”. Pren valors de 0 i 1, així que té una distribució binomial.
- Quines variables són les explicatives?
RESPOSTA: Tenim 5 variables explicatives: “bloc1”, “bloc2”, “bloc3”, “bloc4” i “bloc5”.
Abans d’entrar les dades en R, pensa si cal organitzar de manera diferent l’arxiu d’excel. Un cop tinguis les dades organitzades correctament entra les dades a l’R i guardar l’arxiu en format R.
IMPORTAR DADES load("entregable_4.RData")
GUARDAR DADES save(entregable_4,file="entregable_4.RData")
Ara, per tal de verificar que s’han entrat correctament demana un petit resum de les dades i la seva estructura.
RESUM DADES summary(entregable_4)
sp_invasora bloc1 bloc2 bloc3 bloc
Min. :0.0000 Min. : 0 Min. : 8.00 Min. :10 Min. :20.
1st Qu.:0.0000 1st Qu.:20 1st Qu.:14.00 1st Qu.:10 1st Qu.:25.
Median :1.0000 Median :20 Median :18.00 Median :10 Median :30.
Mean :0.5333 Mean :19 Mean :16.76 Mean :10 Mean :27.
3rd Qu.:1.0000 3rd Qu.:20 3rd Qu.:20.00 3rd Qu.:10 3rd Qu.:30.
Max. :1.0000 Max. :20 Max. :20.00 Max. :10 Max. :30.
bloc
Min. : 3.
1st Qu.:20.
Median :20.
Mean :19.
3rd Qu.:20.
Max. :20.
ESTRUCTURA DE LES DADES str(entregable_4)
Classes 'tbl_df', 'tbl' and 'data.frame': 45 obs. of 6 variables:
$ sp_invasora: num 1 1 1 1 1 1 1 1 1 1 ...
$ bloc1 : num 20 20 20 20 20 20 5 20 20 20 ...
$ bloc2 : num 14 11 15 15 20 16 13 20 20 13 ...
$ bloc3 : num 10 10 10 10 10 10 10 10 10 10 ...
$ bloc4 : num 30 30 30 30 30 30 30 30 30 30 ...
$ bloc5 : num 20 20 20 20 20 20 20 20 20 20 ...
- Creus que l’estructura i el tipus de variables del teu arxiu és l’adient? Justifica la resposta
RESPOSTA: No, la variable “sp_invasora” conta com si fos una variable numèrica hi hauria de ser un factor ja que només indica presència i absència però amb els números 1 i 0. Per tal de que les dades siguin adients, haurem de transformar-la i guardar-la al document perquè l’Rstudio la pugui llegir bé: sp_invasora_factor <- factor(entregable_4$sp_invasora) entregable_4$sp_invasora <- factor(entregable_4$sp_invasora)
Finalment, comprovem que ara sí es llegeix com un factor: summary(entregable_4)
sp_invasora bloc1 bloc2 bloc3 bloc
0:21 Min. : 0 Min. : 8.00 Min. :10 Min. :20.
1:24 1st Qu.:20 1st Qu.:14.00 1st Qu.:10 1st Qu.:25.
Median :20 Median :18.00 Median :10 Median :30.
Mean :19 Mean :16.76 Mean :10 Mean :27.
3rd Qu.:20 3rd Qu.:20.00 3rd Qu.:10 3rd Qu.:30.
Max. :20 Max. :20.00 Max. :10 Max. :30.
bloc
Min. : 3.
1st Qu.:20.
Median :20.
Mean :19.
3rd Qu.:20.
Max. :20.
Sí, ja no calcula res i només ens diu la quantitat de 0 i 1 que té la variable “sp_invasora”.
Finalment necessitem un gràfic que mostri la distribució dels valors de totes les variables i que et pugui donar una idea de possibles relacions entre elles.
GRÀFIC pairs(~entregable_4$sp_invasora+entregable_4$bloc1+entregable_4$bloc2+entregable_4$bloc3+entregable_4$bl
bloc2 -0.33376 0.13149 -2.538 0.0111 *
bloc3 NA NA NA NA
bloc4 0.17796 0.14157 1.257 0.
bloc5 0.02023 0.17303 0.117 0.
---
Signif. codes: 0 '' 0.001 '' 0.01 '' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 62.183 on 44 degrees of freedom
Residual deviance: 52.140 on 40 degrees of freedom
AIC: 62.
Number of Fisher Scoring iterations: 4
Només les variables amb un p-valor < 0.05 són significatives. En aquest cas, en primera instància, només trobariem la variable “bloc2” (0.0111). Aixó pot ser degut a que les variables han de ser transformades. Probarem a transformar-les amb logaritmes i elevant-les al quadrat en el PAS 2. Però també podem fer un gràfic del model per veure si compleix o no els supòsits.
plot(sp_invasora_glm)
Predicted values
Residuals
glm(sp_invasora_factor ~ bloc1 + bloc2 + bloc3 + bloc4 + bloc5)
Residuals vs Fitted
22
19
14
Theoretical Quantiles
Std. Pearson resid.
glm(sp_invasora_factor ~ bloc1 + bloc2 + bloc3 + bloc4 + bloc5)
Normal Q−Q
22
19
14
Predicted values
Std. Pearson resid.
glm(sp_invasora_factor ~ bloc1 + bloc2 + bloc3 + bloc4 + bloc5)
Scale−Location
22
19 14
Histogram of entregable_4$bloc
bloc
Frequency
hist(log(entregable_4$bloc1), xlab='log.bloc1')
Histogram of log(entregable_4$bloc1)
log.bloc
Frequency
hist(sqrt(entregable_4$bloc1), xlab='sqrt.bloc1')
Histogram of sqrt(entregable_4$bloc1)
sqrt.bloc
Frequency
A la vari- able “bloc1” no hi ha canvis significatius així que la deixarem sense transformar.
hist(entregable_4$bloc2, xlab='bloc2')
Histogram of entregable_4$bloc
bloc
Frequency
hist(entregable_4$bloc3, xlab='bloc3')
Histogram of entregable_4$bloc
bloc
Frequency
hist(log(entregable_4$bloc3), xlab='log.bloc3')
Histogram of log(entregable_4$bloc3)
log.bloc
Frequency
hist(sqrt(entregable_4$bloc3), xlab='sqrt.bloc3')
Histogram of sqrt(entregable_4$bloc3)
sqrt.bloc
Frequency
A la vari- able “bloc3” no hi ha canvis significatius així que la deixarem sense transformar.
hist(entregable_4$bloc4, xlab='bloc4')
Histogram of entregable_4$bloc
bloc
Frequency
hist(entregable_4$bloc5, xlab='bloc5')
Histogram of entregable_4$bloc
bloc
Frequency
hist(log(entregable_4$bloc5), xlab='log.bloc5')
Histogram of log(entregable_4$bloc5)
log.bloc
Frequency
hist(sqrt(entregable_4$bloc5), xlab='sqrt.bloc5')
Histogram of sqrt(entregable_4$bloc5)
sqrt.bloc
Frequency
A la vari- able “bloc5” no hi ha canvis significatius així que la deixarem sense transformar. Finalment veiem que no val la pena transformar cap variable i les deixarem tal i com les teníem al principi. PAS 3 Passem a mirar ara si tenim problemes de col·linealitat ja que quan es fa una regressió múltiple sempre s’ha de vigilar de no tenir col·linealitat. Per fer-ho utilitzarem el variance inflation factor. Com que la variable “bloc3” és constant, no sortirà.
(AQUÍ VAIG POSAR LA FUNCIÓ “vif(sp_invasora_glm)” I EM VA SORTIR LA TAULA VIF, PERÒ AL PASSAR-HO A PDF EM SORTIA “Error in vif(sp_invasora_glm) : could not find function”vif" Calls:... handle -> withCallingHandlers -> withVisible -> eval -> eval Execution halted AIXÍ QUE US POSSARÉ ELS VALORS QUE M’HA DONAT LA TAULA VIF I HO BORRARÉ PER PODER VER EL PDF): bloc1: 1.244701 bloc2: 1.256663 bloc4: 1.150236 bloc5: 1.009136 No hi ha cap variable que tingui el vif per sobre de 5, així que podem assumir que no tenim problemes de col·linealitat. PAS 4 Ara em de veure si tenim overdispresion (o underdispersion). summary(sp_invasora_glm)
Call:
glm(formula = sp_invasora_factor ~ bloc1 + bloc2 + bloc3 + bloc4 +
bloc5, family = binomial, data = entregable_4)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.9271 -1.0023 0.3828 0.8590 1.
Coefficients: (1 not defined because of singularities)
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.32381 5.64715 -0.057 0.
summary(sp_invasora_glm)
Call:
glm(formula = sp_invasora_factor ~ bloc1 + bloc2 + bloc3 + bloc4 +
bloc5, family = binomial, data = entregable_4)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.9271 -1.0023 0.3828 0.8590 1.
Coefficients: (1 not defined because of singularities)
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.32381 5.64715 -0.057 0.
bloc1 0.04144 0.10836 0.382 0.
bloc2 -0.33376 0.13149 -2.538 0.0111 *
bloc3 NA NA NA NA
bloc4 0.17796 0.14157 1.257 0.
bloc5 0.02023 0.17303 0.117 0.
---
Signif. codes: 0 '' 0.001 '' 0.01 '' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 62.183 on 44 degrees of freedom
Residual deviance: 52.140 on 40 degrees of freedom
AIC: 62.
Number of Fisher Scoring iterations: 4
Com hem vist al principi, només la variable “bloc2” és la que té un p-valor menor de 0.05, o sigui que és la única variable significativa. Això ens mostra que podem significar el model. Ho farem a partir del mètode del model únic.
PAS 8 Per fer-ho utilitzarem la simplificació automàtica via la funció “stepAIC”.
step(sp_invasora_glm, method="bacward")
Start: AIC=62.
sp_invasora_factor ~ bloc1 + bloc2 + bloc3 + bloc4 + bloc
Step: AIC=62.
sp_invasora_factor ~ bloc1 + bloc2 + bloc4 + bloc
Df Deviance AIC
- bloc5 1 52.154 60.
- bloc1 1 52.275 60.
- bloc4 1 53.827 61.
52.140 62.
- bloc2 1 60.528 68.
Step: AIC=60.
sp_invasora_factor ~ bloc1 + bloc2 + bloc
Df Deviance AIC
- bloc1 1 52.288 58.
- bloc4 1 53.967 59.
52.154 60.
- bloc2 1 60.837 66.
Step: AIC=58.
sp_invasora_factor ~ bloc2 + bloc
Df Deviance AIC
- bloc4 1 53.968 57.
52.288 58.
- bloc2 1 61.329 65.
Step: AIC=57.
sp_invasora_factor ~ bloc
Df Deviance AIC
53.968 57.
- bloc2 1 62.183 64.
Call: glm(formula = sp_invasora_factor ~ bloc2, family = binomial,
data = entregable_4)
Coefficients:
(Intercept) bloc
5.1960 -0.
Degrees of Freedom: 44 Total (i.e. Null); 43 Residual
Null Deviance: 62.
Residual Deviance: 53.97 AIC: 57.
Observem com passem d’un AIC de 62.14 a un de 57.97 (com més petit aquest valor, millor és el model) progressivament mentres va eliminant les variables que no són significatives. Doncs, el model final ajustat només tindrà la variable “bloc2” com a variable explicativa de la nostra variable resposta “sp_invasora_factor”.
PAS 9 Fem el model final.
sp_invasora_glm_final<-glm(sp_invasora_factor~bloc2,family="binomial",data=entregable_4) summary(sp_invasora_glm_final)
Call:
glm(formula = sp_invasora_factor ~ bloc2, family = "binomial",
data = entregable_4)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.0285 -0.8762 0.5229 1.0004 1.
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 5.1960 2.0511 2.533 0.0113 *
bloc2 -0.2978 0.1169 -2.547 0.0109 *
---
petits) trobem menys presència de la espècie invasora a la llacuna estudiada.
