library(tibble)
library(dplyr)
library(ggplot2)

Referències

L. Cayuela & M. De la Cruz, Análisis de datos ecológicos en R, 2022.
A. García Pérez, Cuadernos de estadística aplicada: biología y ciencias ambientales, 2023.

Introducció

En els models lineals la variable de resposta sempre ha de ser quantitativa, mentre que les variables explicatives poden ser quantitatives o categòriques.

Regressió simple: una sola variable explicativa contínua.
Regressió múltiple: dues o més variables explicatives continues.
ANOVA monofactorial: una sola variable categòrica
ANOVA multifactorial: dues o més variables categòriques.
ANCOVA: variables categòriques i contínues.

Tots cinc son casos particulars de una única categoria de model que denominem models linials.

Tots tenen en comú dos elements:

La variable de resposta \(y\) és contínua.
Sobre el model es fan 4 supòsits: independència, linealitat, normalitat, homoscedasticitat.

En R tots aquests 5 models es poden ajustar amb una única funció, lm(), amb unes manera d’interpretar els models és molt semblant.

Regressió simple

La regressió simple en R es duu a terme amb lm() i l’argument principal és formula que relaciona \(y\) i \(x\). Els models linials utilitzen de mínims quadrats per estimar-ne els coeficients.

Exemple la Cladonia rangiferina

El liquen Cladonia rangiferina creix al sotabosc dels pinars. Investigarem quina és la influència té el gruix de la capa de humus en la cobertura (en %) de l’espècie, establint la hipòtesi nula \(H_0: \beta_1=0\).

Utilitzarem el dataset cladonia de la llibreria ADER

library(ADER)
data("cladonia")

plot(cladrang~Humdepth, data=cladonia,
     xlab= "Gruix capa de humus (cm)",
     ylab = "Cobertura de Cladonia (%)")

Ajustament del model amb `lm()`

lm.clad_H <- lm(cladrang~Humdepth, data=cladonia)

Després d’ajustar el model podem analitzar-lo des de vessants diferents. Per fer-ho tenim a disposició dues funcions anova() i summary().

Test d’hipòtesi i variació explicada pel model

La primera cosa que cal fer un cop ajustat el model és verificar si el model explica significativament les nostres dades. Això es pot fer amb la taula d’anàlisi de la variança que reumeix les fonts de variació de les dades que han quedat explicades (o no) pel model. Això es du a terme amb la funció anova().

anova(lm.clad_H)

## Analysis of Variance Table
## 
## Response: cladrang
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## Humdepth   1 2389.8 2389.84  19.305 0.0002306 ***
## Residuals 22 2723.5  123.79                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Estructura de la taula ANOVA

La taula presenta una sèrie de columnes que resumeixen la partició de la suma dels quadrats. Les files són una per a cada variable més una de la variança residual (Residuals). El sumatori de les sumes de quadrats (Sum Sq) representa la suma dels quadrats de la variable \(y\) (\(SS_Y\)) que correspon a \[ SS_Y= \sum_{i=1}^n (Y_i-\bar Y_i)^2 \] A la nostra taula seria la suma del la columna: 2389.8 + 2723.5 = 5113.32.

A partir d’aquesta dada podem calcular el coeficient de determinació \(R^2\), que sabem que s’obté de la suma residual dels quadrats \(RSS\) y de la suma dels quadrats de la variable \[ \small R^2=1-\frac{RSS}{SS_Y}=1-\frac{\sum_{i=1}^n(y_i-\hat y_i)^2}{\sum_{i=1}^n(y_i-\bar y_i)^2} \]

Considerant que la suma dels quadrats explicada pel model de regressió \(SS_{reg}\) és \(SS_{reg}=SS_Y-RSS\), \(R^2\) serà \[ R^2=\frac{SS_{reg}}{SS_Y} \] i, en el nostre cas, serà

R2 <- anova(lm.clad_H)$`Sum Sq`[1] / sum(anova(lm.clad_H)$`Sum Sq`)
R2

## [1] 0.4673767

Dividint les sumes dels quadrats per els graus de llibertat (Df) obtenim les variances o sumes mitjanes dels quadrats (Mean Sq). Aquesta informació s’utilitza per fer el test d’hipòtesi en models linials, car podem obtenir l’estidístic \(F\) que és \[ \small F= \frac{SS_{reg}/1}{RSS/(n-2)} \] que justament és 2389.84/(123.79)=19.305 de la columna F value. La columna Pr(>F) ens dona el valor p, és a dir, la probabilitat que l’estadístic \(F\) caigui dins de la distribució nul·la de valors. Una probabilitat baixa (típicament <0.05) ens fa rebutjar la \(H_0\) que \(\beta_1=0\), com en aquest cas.

Interpretació del model

Llavors, si es rebutja la hipòtesi nul·la, el pas següent és explorar els coeficients estimats del model amb la funció summary().

summary(lm.clad_H)

## 
## Call:
## lm(formula = cladrang ~ Humdepth, data = cladonia)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -22.833  -6.973  -2.367   5.762  24.399 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   49.937      8.008   6.236 2.82e-06 ***
## Humdepth     -15.337      3.491  -4.394 0.000231 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.13 on 22 degrees of freedom
## Multiple R-squared:  0.4674, Adjusted R-squared:  0.4432 
## F-statistic:  19.3 on 1 and 22 DF,  p-value: 0.0002306

Fixar-se en la taula dels coeficients. Els valors es faciliten a la columna Estimate. Segons aquesta la relació entre la cobertura (\(y\)) de la Cladonia i el gruix de la capa de humus (\(x\)) és \[ y= 49.937-15.337x \]

Fixar-se també en l’estadístic t value que és igual a Estimate/Std. Error. El valor absolut d’aquest estadístic es compara amb una distribució t amb els graus de llibertat de la regressió (degrees of freedom). Per a obtenir un p valor d’un test de 2 cues a la columna Pr(>|t|).

A sota la taula de coeficients trobem informació addicional: L’error estàndard residual (Residual standard error), que no és més que l’arrel quadrada del valor de Mean Sq de la columna Residuals a la taula ANOVA: es tracta doncs de la desviació estàndard (\(\sigma\)) dels valors residuals, que en aquest cas és \(\sigma=11.2\), bastant alta, indicant dispersió de les dades al voltant de la recta. Els graus de llibertat són \(n-2\text{ coeficients}=24-2=22\). Finalment hi ha un estadístic \(F\) que testa la hipòtesi que tots els coeficients de la regressió son iguals a 0. En el cas de la regressió simple coincideix amb l’\(F\) de la taula ANOVA. A més a més, com que només hi ha un coeficient, aquí \(F=t^2=-4.394^2=19.3\).

In capítol a part mereix \(R^2\), que ja hem vist a la taula ANOVA, i \(R^2~\text{ajustat}\), que en regressió lineal simple és pràcticament igual a \(R^2\). Tanmateix, l’augment important del nombre de paràmetres (és a dir si el nobre de paràmetres s’aproxima al nombre d’observacions), provoca un augment de \(R^2\) que s’aproximarà a 1, encara que cap de les variables tingui efectes en la variable de resposta. En aquest cas, és molt més eficaç utilitzar \(R^2\text{-adjusted}\).

Representació gràfica del model

par(mfrow=c(1,2))

plot(cladonia$Humdepth,cladonia$cladrang,
     xlab = "",ylab = "")
abline(coef = lm.clad_H$coefficients) # Recta de regressió
title(main = "Regressió simple Cladonia",
      cex.main=1,
      x= "Gruix capa d'humus (cm)",
      y= "Cobertura Cladonia (%)")


# Per a generar els intervals de confiança
x.ic <- seq(0,4, 0.01) # genera una sequència de valors x per l'int. confiança
predict.clad <- predict(lm.clad_H, 
                        newdata = list(Humdepth=x.ic), 
                        interval = "confidence")


plot(cladonia$Humdepth,cladonia$cladrang,xlab = "",ylab = "")
title(main = "Regressió simple Cladonia amb \nintervals de confiança",
      cex.main=1,
      x= "Gruix capa d'humus (cm)",
      y= "Cobertura Cladonia (%)")
lines(x.ic, predict.clad[,1]) # Recta de regressió
lines(x.ic, predict.clad[,3], col="red", lty=2) # Int. de confiança sup.
lines(x.ic, predict.clad[,2], col="darkgreen", lty=2) # Int. de confiança inf.

Revisió dels supòsits del model

Pel que fa al supòsit d’independència, es valora en base a l’avaluació del model experimental, és a dir, en base a com es generen les dades.

Gràfics de diagnòstic de residus

Els supòsits de linealitat, normalitat i homoscedasticitat es poden comprobar mitjant els gràfics de diagnòstic de residus:

par(mfrow= c(2,2))
plot(lm.clad_H)

`Residual vs Fitted`

El gràfic Residual vs Fitted mostra els valors residuals contra els valors esperats. S’utilitza per comprovar els supòsits de linealitat, però també per comprovar si els residus son normals i homoscedàstics.

Si tot està bé, s’hauria d’observar un núvol de punts sense cap patró aparent, com a l’exemple aquí sota.

alcohol <- read.table("datosr/alcohol.txt", header = T)

lm.alc <- lm(strength~alcohol, data=alcohol)

plot(lm.alc,1)

Comprovem que aquest no és el cas, amb una pauta amb forma d’embut que es va obrint per la dreta.

`Normal Q-Q`

El gràfic Normal Q-Q representa els residus estandarditzats contra quantils teòrics. Si els residus están normalment distribuïts els punts seguiran la diagonal dibuixada al gràfic. Si se-n’allunyen hi ha violació del supòsit de normalitat. Tot i així, és relativament freqüent observar una mica d’allunyament, com en aquest cas.

`Scale-location`

Aquesta gràfica presenta la relació de l’arrel quadrada dels valors absoluts dels residus estandardidtats contra els valors esperats. Aquest gràfic permet comprovar específicament el supòsit d’homoscedasticitat: la situació desitjable és un núvol de punts sense pautes abservables, com el gràfic aquí sota

plot(lm.alc,3, add.smooth = T)

i no és el cas del nostre model on s’observa una clara tendència a augmentar la variança residual augmentant la mitjana.

`Residuals vs Leverage`

Aquest gràfic s’utilitza per identificar observacions amb potencial de influència sobre la pendent de la recta de regressió: són els punts palanca (no uns simples outliers) que es troben als extrems de les x i de les y. Són dades que ens han de preocupar perquè poden influir de manera significativa en els resultats de l’anàlisi.

La configuració ideal pot ser la que s’ensenya aquí sota.

plot(lm.alc,5, add.smooth = T)

Es consideren dades molt influents, amb capacitat d’influir en l’anàlisi quan superen la línia de 1 de la distància de Cook. En el nostre model, cap dada supera aquest marge, tot i que n’hi ha una, la 5, que supera 0.5 de la distància.

Tests addicionals

Test de normalitat de Shapiro-Wilk

Aquest comprova la normalitat dels residus. L’hipòtesi nul·la és que les dades són distribuïdes nromalment normalment.

shapiro.test(lm.clad_H$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  lm.clad_H$residuals
## W = 0.96809, p-value = 0.6201

En el cas de la nostra regressió el p valor tan alt ens confirma la normalitat dels valors residuals.

Test de linealitat RESET

Per al test de linealitat, està disponible el test RESET (Ramsey, 1969). La hipòtesi nul·la d’aquest test és que la relació entre x i y és lineal. Per aquest test podem fer servir la funció resettest() de paquet lmtest.

lmtest::resettest(lm.clad_H)

## 
##  RESET test
## 
## data:  lm.clad_H
## RESET = 1.1693, df1 = 2, df2 = 20, p-value = 0.3309

En aquest cas el valor alt del p valor recomana acceptar la \(H_0\) de linealitat de la relació.

Test d’homoroscedasticitat de Beusch-Pagan

Podem testar el supòsit de la constança de la variança amb el test de Beusch-Pagan. La hipòtesi nul·la del test és que la variança és constant. El podem dur a terme amb la funció ncvTest del paquet car.

car::ncvTest(lm.clad_H)

## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 5.388179, Df = 1, p = 0.020274

Tal com ja havíem observat amb el gráfic “Scale-Location”, es detecta un problema d’heteroscedasticitat en els valors residuals de la nostra recta ajustada.

Comentari

És important ser conscients que la violació d’algun o varis supòsits de del model és més la tònica que l’excepció. La solució es pot trobar en algun tipus de transformació.

Què fer quan es violen els supòsits?

Sovint ens podem trobar amb violacions d’un o més supòsits, per exemple, quan la relació entre x i y no és linial. Un altre cas bastant comú és la variana no constant.En ambdós casos la simple transformació de la variable resposta pot resoldre el problema. La transormació logarítmica sol ser una pràctica comú, però també se’n fan servir altres, com l’arrel quadrada, la inversa, o la transformació de Box-Cox.

En cualsevol cas convé no aplicar transformacions molt complicades a què no siguem capaços de trobar després un significat biològic molt clar.

De vegades pot ser convenient transformar també la variable explicativa: la transformació logarítmica tan en la variblae de resposta com l’explicativa es coneix com transformació log-log.

Amb el següent exemple veiem com, utilitzant una transfomació log-log, podem linealitar i les dimensions de l’illa d’un arxipèlag i alhora controlada la variana residual.

library(faraway)

data("gala")

head(gala)

##              Species Endemics  Area Elevation Nearest Scruz Adjacent
## Baltra            58       23 25.09       346     0.6   0.6     1.84
## Bartolome         31       21  1.24       109     0.6  26.3   572.33
## Caldwell           3        3  0.21       114     2.8  58.7     0.78
## Champion          25        9  0.10        46     1.9  47.4     0.18
## Coamano            2        1  0.05        77     1.9   1.9   903.82
## Daphne.Major      18       11  0.34       119     8.0   8.0     1.84

lm.gala <- lm(Species~Area, data = gala)
summary(lm.gala)

## 
## Call:
## lm(formula = Species ~ Area, data = gala)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -99.495 -53.431 -29.045   3.423 306.137 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 63.78286   17.52442   3.640 0.001094 ** 
## Area         0.08196    0.01971   4.158 0.000275 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 91.73 on 28 degrees of freedom
## Multiple R-squared:  0.3817, Adjusted R-squared:  0.3596 
## F-statistic: 17.29 on 1 and 28 DF,  p-value: 0.0002748

Mirant la taula de summary, el model sembla correcte, amb un \(F\) molt gran, amb uns coeficients amb valor p infinitesimal que també en confirmen la validesa.

Tanmateix, observant la recta del model i el gràfic Residuals vs Fitted, veiem que el model presenta problemes de linialitat i de heteroscedasticitat.

par(mfrow=c(1,2))
plot(gala$Area,gala$Species)
abline(lm.gala)

plot(lm.gala,1, add.smooth = F)

Si executem el test RESET de linialitat, es confirma la nostra impressió visual sobre la falta de linialitat, atès el valor p infinitesimal.

lmtest::resettest(lm.gala)

## 
##  RESET test
## 
## data:  lm.gala
## RESET = 19.783, df1 = 2, df2 = 26, p-value = 5.995e-06

També quedaría rebutjada la hipòtesi de homogeneïtat de la variança amb el test de Beusch-Pagan

car::ncvTest(lm.gala)

## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 3.929863, Df = 1, p = 0.047436

Provarem per tant una transformació log-log, trasformant logarítmicament tant la variable de resposta (Species), com la variable explicativa (Area)

lm.log <- lm(log(Species)~log(Area), data = gala)

summary(lm.log)

## 
## Call:
## lm(formula = log(Species) ~ log(Area), data = gala)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.5442 -0.4001  0.0941  0.5449  1.3752 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.9037     0.1571  18.484  < 2e-16 ***
## log(Area)     0.3886     0.0416   9.342 4.23e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7842 on 28 degrees of freedom
## Multiple R-squared:  0.7571, Adjusted R-squared:  0.7484 
## F-statistic: 87.27 on 1 and 28 DF,  p-value: 4.23e-10

Tal com podem veure, els paràmetres del model presenten una validesa igual o fins i tot superior (con en el cas de \(R^2\)), per exemple.

Observant els dos gràfics que hem utilitzat pel model anterior

par(mfrow=c(1,2))
plot(log(gala$Area), log(gala$Species))
abline(lm.log)

plot(lm.log,1, add.smooth = F)

Podem observar en el gràfic de l’esquerra una linealitat més clara, i una homogeneïtat més evident en la variança dels valors residuals.

Cosa que queda confirmada pels test de linealitat i d’homogeneïtat de la variança (encara que amb algun dubte aquest últim).

lmtest::resettest(lm.log)

## 
##  RESET test
## 
## data:  lm.log
## RESET = 0.067046, df1 = 2, df2 = 26, p-value = 0.9353

car::ncvTest(lm.log)

## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 3.467564, Df = 1, p = 0.062583

Finalment, si les transformacions no donen efecte, caldrà canviar tipologia de model.

En el cas d’heteroscedasticitat, es podria intentar definir una estructura d’errors no normals amb algun model linial generalitzat.
Alguns casos de falta de linealitat es poden tractar amb regressions polinòmiques o models additius generalitzats que permeten modelar relacions no linials entre \(x\) i \(y\).
Finalment, quan hi hagi estructures de no independència podem incorporarles amb la inclusió de factors aleatoris en el context de models mixtes.

Regressió múltiple

La regressió múltiple és un cas particular de regressió lineal on tenim diverses variables explicatives contínues. \[ \small y_i=\beta_0+\beta_1x_1+...+\beta_kx_k+\epsilon_i~~~\epsilon_i \sim N(0,\sigma^2) \]

Ajustament del model

L’ajustament del model de regressió múltiple és similar al del model regressió simple. Continuant amb l’exemple sobre la Cladonia rangiferina, i com la seva cobertura es veu afectada per diversos factors ambientals, si volquessim a més a més saber com impacten la concentració de \(N\) i de \(P\), tindríem, en principi, tres hipòtesis nul·les i tres alternatives per comprovar: si la pendent de la recta de cada una d’aquestes variables és o no igual a 0.

Un error comú és plantejar una regressió lineal simple per a cada variable: no s’ha de fer, excepte per fer unes exploracions gràfiques (per exemple, una matriu exploratòria de diagrames de dispersió), perquè els efectes d’una variable de sobre la variable de resposta poden estar condicionats pels efectes d’una altra variable. D’altra banda, realitzant múltiples tests augmenta considerablement l’error de tipus I.

Anàlisi visual previ

Utilitzarem la funció ggpairs() de GGally que és molt más immediada de fer servir respecte a pairs().

GGally::ggpairs(cladonia[,c(1:3,14)], 
        lower = list(continuous = GGally::wrap(
          "smooth_lm",
          color = "blue"))) +
  theme_bw()

Encara que no es veuen unes relacions molt clares, podem dir que l’ajustament de rectes de regressió linear simple mostra una relació més o menys lineal, positiva amb el nitrogen i negativa amb les altres dues variables.

Ajustament

El següent pas és ajustar a les dades un model de regressió múltiple. També es fa servir la funció lm(), amb la diferència que, en la fórmula s’inclouen vàries variables explicatives separedes pel signe \(+\).

lm.clad <- lm( cladrang ~ N + P + Humdepth, data = cladonia)

Test d’hipòtesi i variació explicada pel model

Seguint la metodologia ja marcada, un cop ajustat el model a les dades, mitjançant la funció anova(), examinem si les variables explicatives tenen un efecte significatiu en la variable de resposta.

anova(lm.clad)

## Analysis of Variance Table
## 
## Response: cladrang
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## N          1  776.75  776.75  11.869  0.002560 ** 
## P          1  843.17  843.17  12.883  0.001833 ** 
## Humdepth   1 2184.49 2184.49  33.379 1.184e-05 ***
## Residuals 20 1308.91   65.45                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sumes de quadrats de tipus I

Quan l’anàlisi es duu a terme seguint la seqüència:

el primer test comprova l’efecte de N, ignorant P i Humdepth.
el segon test contrasta l’efecte de P depsrés d’aher tingut en compte N i ignorant els efecte d’Hum depth,
el tercer contrasta els efectes d’Humdepth d’haver ajustat els efectes de Ny P,

aquesta manera de calcular la contribució relativa de cada variable per explicar la variabilitat de la resposta es denomina Sumes de quadrats de tipus I, que és la que R executa per defecte.

Les sumes dels quadrats de tipus I s’obtenen calculant la reducció en la suma residual dels quadrats a mesura que afegim termes al model de manera seqüencial.

Exemple

Les variables \(A\) i \(B\) expliquen per elles soles la variabilitat \(SS_A\) i \(SS_B\). Però, a més a més, a causa d’una certa colinealitat existent entre \(A\) i \(B\), part de la variabilitat d’\(y\) és explicada conjuntament per totes dues \(SS_{AB}\). Si la primera variable que s’introdueix al model és \(A\), aquesta explicará \(SS_{A}\) que inclourà també la variabilitat compartida amb \(B\), \(SS_{AB}\). Per part seva, la variable \(B\) explicarà la variabilitat de \(y\), \(SS_{B|A}\), és a dir la part no explicada conjuntament per \(A\) y \(B\). El mateix passarà si primer introduïm \(B\) (és a dir la primera que entra es queda la variabilitat compartida).

Això vol dir que en aquest cas l’ordre importa, ja que afecta els test di significació que estan basats en les sumes de quadrats. En el cas d’independència de les variables obtindrem resultats molts similars independentment de quines variables entrin primer al model.

Sumes de quadrats de tipus III

Les sumes de de quadrats de tipus III es denominen sumes de quadrats marginals.

Aquest tipus de sumes de quadrats es calcula per diferència: apartant cada vegada una variable, la diferència entre la suma dels quadrats del model complet i la suma dels quadrats del model sense aquesta variable. Això te per conseqüència que la variabilitat compartida \(SS_{AB}\) queda exclosa del càlcul de les sumes dels quadrats explicades per cada variable. És a dir que cada variable contribueix a les sumes dels quadrats només amb la parti marginal de la variança.

R executa per defecte les sumes de tipus I. Si volem utilitzar les sumes de de tipus III, hem de fer servir la funció Anova() de la liberia car.

car::Anova(lm.clad, type= "III")

## Anova Table (Type III tests)
## 
## Response: cladrang
##              Sum Sq Df F value    Pr(>F)    
## (Intercept)  936.96  1 14.3167  0.001166 ** 
## N            631.17  1  9.6443  0.005575 ** 
## P            411.21  1  6.2832  0.020935 *  
## Humdepth    2184.49  1 33.3787 1.184e-05 ***
## Residuals   1308.91 20                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

En qualsevol cas, cal fer l’incís que cal tenir en compte que els coeficients no canviaran utilitzand diferents tipus de sumes de quadrats. Canvien les hipòtesis testades.

Quan utilitzar un tipus o l’altre

En principi, si el nostre model només conté variables continues el tipus de sumes de quadrats que utilitzaem no és rellevant sempre i quan no existeixi colinealitat entre variables explicatives.

Si existeix colinealitat, llavors ens hem de preguntar si hi ha jerarquització dels efectes sobre la variable de resposta. Per exemple, si volem comprovar l’efecte de la temperatura i de l’ús antròpic sobre la cobertura d’una espècies, podem considerar la temperatura un factor més general y global i l’ús antròpic com más localitzat. En aquest cas utilitzarem el Tipus I definint l’ordre de les variables segons aquest criteri. En cas que considerem que no existeixi cap tipus de jerarquització farem servir el Tipus III o el Tipus II (molt similar).

Interpretació i representació gràfica del model

La funció summary() ens proporciona els coeficients i els seus errors estàndard.

summary(lm.clad)

## 
## Call:
## lm(formula = cladrang ~ N + P + Humdepth, data = cladonia)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.1670  -5.4829  -0.2535   6.5184  14.5110 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  40.3875    10.6740   3.784  0.00117 ** 
## N             0.9861     0.3175   3.106  0.00557 ** 
## P            -0.2970     0.1185  -2.507  0.02093 *  
## Humdepth    -14.9431     2.5865  -5.777 1.18e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.09 on 20 degrees of freedom
## Multiple R-squared:  0.744,  Adjusted R-squared:  0.7056 
## F-statistic: 19.38 on 3 and 20 DF,  p-value: 3.915e-06

\[ \tt cladrang= \beta_0+ \beta_1N+\beta_2P+\beta_3Humdepth\\ \tt cladrang= 40.3875+ 0.9861N -0.297P-14.9431Humdepth \]

Els coeficients sempre s’han d’interpretar en el context de les altres variables: per això se’ls denomina coeficients parcials.

Quan es facin prediccions cal sempre tenir present que, en el cas de regressió múltiple, no s’ajusta una recta, sinó un hiperplà. En aquests casos, quan es representi els efectes de cada variable explicativa sobre la de resposta, serà bona pràctica mantenir constants la resta de variables explicatives, fent servir, per exemple el valor de la seva mitjana. I així, generar una taula de dades.

# Crear una taula de valors de les var. explicatives on la var. on les altres var. repeteixen el seu valor mitjà.

n <- 100

dades.noves <- data.frame(N=seq(min(cladonia$N),max(cladonia$N), length.out=n),
           P= rep(mean(cladonia$P)),
           Humdepth= rep(mean(cladonia$Humdepth)))

# Predicció de la variable de resposta als valors de N
resposta.N <- predict(lm.clad, newdata = dades.noves)

# Visualització amb plot()
plot(cladonia$N,cladonia$cladrang, 
     main = "Efecte parcial de la concentració de N \nen la cobertura de Cladonia",
     ylab = "Cobertura de Cladonia (%)",
     xlab = "Concentració de N")
lines(dades.noves$N,resposta.N)

Per poder examinar ràpidament i alhora els efectes parcials de les diferents covariables, es pot fer servir la funció allEffects del paquet effects

library(effects)
plot(effects::allEffects(lm.clad), row=1, col=3)

Coeficients de regressió estandarditzats

Quan, en regressió linial múltiple, a més d’avaluar la significància es vol comparar la magnitud de l’efecte de cada variable, cal fer servir els coeficients restandarditzats, altrament dits coeficients beta. per a això hem d’estandarditzar la resposta i les variables explicatives per a que tinguin mitjana \(\bar x=0\) i desviació estàndar \(s=1\). Tot això s’aconsegueix amb la funció scale

lm.clad.st <- lm( scale(cladrang) ~ scale(N) + scale(P) + scale(Humdepth), 
                  data = cladonia)
summary(lm.clad.st)

## 
## Call:
## lm(formula = scale(cladrang) ~ scale(N) + scale(P) + scale(Humdepth), 
##     data = cladonia)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.7489 -0.3677 -0.0170  0.4372  0.9732 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -1.906e-16  1.108e-01   0.000  1.00000    
## scale(N)         3.656e-01  1.177e-01   3.106  0.00557 ** 
## scale(P)        -2.977e-01  1.188e-01  -2.507  0.02093 *  
## scale(Humdepth) -6.661e-01  1.153e-01  -5.777 1.18e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5426 on 20 degrees of freedom
## Multiple R-squared:  0.744,  Adjusted R-squared:  0.7056 
## F-statistic: 19.38 on 3 and 20 DF,  p-value: 3.915e-06

De manera que el test de Wald y el de la \(F\) no canvien, com tampoc l’\(R^2\). Canvien els coeficients: \(\beta_0\) es posa pràcticament a 0 i els altres coeficients són comparables: per exemple, es pot veure que l’efecte gruix de la capa d’humus seria quasi el doble en termes absoluts que el del Nitrogen.

Multicolinealitat

Anàlisi amb visualització

GGally::ggpairs(cladonia[,c(2:7,14:15)], 
        lower = list(continuous = GGally::wrap(
          "smooth_lm",
          color = "blue"))) +
  theme_bw()

En aquest correlograma es poden observar visualment els alts nivells de colinialitat entre algunes variables.

Anàlisi amb el Factor d’inflació de la variança (Variance Inflation Factor - VIF)

El VIF estableix com augmenta d’un coeficient estimat del model quan les varibles expicatives estan correlacionades entre si.

\[ \small \hat V(\beta_j) = \frac{S^2}{(n-1)S_j^2} \frac{1}{1-R_j^2} \] Com es veu aquest índex és una funció del coeficient de determinació \(R_j^2\). L’arrel quadrada de VIF indica quant s’incrementa l’interval de confiança de cada coeficient estimat \(\beta_j\) amb relació a les dades correlacionades. Un valor de 4 (que és el que es considera el límit d’acceptabilitat) indicaria doncs una duplicació de la variança.

Càlcul de VIF amb R

Amb R, VIF es pot calcular amb la funció vif() del paquet car:

lm.clad.complet <- lm(cladrang ~ ., data = cladonia)

car::vif(lm.clad.complet)

##         N         P         K        Ca        Mg         S        Al        Fe 
##  1.884049  6.184742 11.746873  9.825650  9.595203 18.460717 20.331905  8.849137 
##        Mn        Zn        Mo  Baresoil  Humdepth        pH 
##  5.168278  7.755110  4.575108  2.213443  5.638807  6.910118

Veiem doncs que excepte N i la proporció de sòl nu (Baresoil), la resta de variables presenten un nivell més que significatiu de colinealitat amb una alta inflació de la variança.

Què es pot fer en aquests casos?

La primera cosa, senzilla, que podem fer és seleccionar de tot el model les varibles que tinguin rellevància biològica i de pas, poc correlacionades. Veiem, per exemple el cas de les variables que hem utilitzat inicialment: N, P i Humdepth:

car::vif(lm.clad)

##        N        P Humdepth 
## 1.082762 1.102065 1.038539

Veiem que el valor són, d’aquesta manera, molt propers a 1 i per tant amb poquísima correlació.

Unes alternatives poden ser l’anàlisi de components principals i la regressió parcial.

Introducció als models lineals en ecologia

Notes d’estadística aplicada a les ciències ambientals: Nota Núm. 1

Stefano Sanfilippo

Abril 2025

Referències

Introducció

Regressió simple

Exemple la Cladonia rangiferina

Ajustament del model amb lm()

Test d’hipòtesi i variació explicada pel model

Estructura de la taula ANOVA

Interpretació del model

Representació gràfica del model

Revisió dels supòsits del model

Gràfics de diagnòstic de residus

Residual vs Fitted

Normal Q-Q

Scale-location

Residuals vs Leverage

Tests addicionals

Test de normalitat de Shapiro-Wilk

Test de linealitat RESET

Test d’homoroscedasticitat de Beusch-Pagan

Comentari

Què fer quan es violen els supòsits?

Regressió múltiple

Ajustament del model

Anàlisi visual previ

Ajustament

Test d’hipòtesi i variació explicada pel model

Sumes de quadrats de tipus I

Exemple

Sumes de quadrats de tipus III

Quan utilitzar un tipus o l’altre

Interpretació i representació gràfica del model

Coeficients de regressió estandarditzats

Multicolinealitat

Anàlisi amb visualització

Anàlisi amb el Factor d’inflació de la variança (Variance Inflation Factor - VIF)

Càlcul de VIF amb R

Què es pot fer en aquests casos?

Ajustament del model amb `lm()`

`Residual vs Fitted`

`Normal Q-Q`

`Scale-location`

`Residuals vs Leverage`