Mostra la font per a Manual de R (Estadística)/Relació entre dos variables

= Introducció =

Es presentaran els mètodes estadístics bàsics per avaluar l'associació entre dues variables. En la taula següent es resumeixen les proves estadístiques bàsiques per les tres combinacions possibles, segons si les dues variables són qualitatives, una qualitativa i l'altre quantitativa o les dues quantitatives (<ref name="Schwartz1994">Schwartz D. Méthodes statistiques à l'usage des médecins et des biologistes 4 ed. Paris: Médecine-Sciences Flammarion; 1994.</ref> pàg. 271):

{| class="wikitable"
|-
! Tipus de variables !! Prova estadística !! Condicions d'aplicació !! Conclusions possibles
|-
| 2 qualitatives || Khi-quadrat || Efectius teòrics suficientment grans || Existència d'associació entre les dues variables.
|-
| 1 qualitativa 1 quantitativa || Anàlisi de la variància (o ''t''-Student si la variable qualitativa només te dos valors possibles).  || Pels diferents valors de la variable qualitativa, la quantitativa te que presentar una distribució normal de la mateixa variància. || Existència d'associació entre les dues variables. La diferència entre les mitjanes informa sobre la importància de l'associació.
|-
| 2 quantitatives || Recta de regressió o correlació. || La variable dependent ha de tenir una distribució normal amb una variància constant pels diferents valors de la variable independent. || Existència d'associació lineal entre les dues variables. Si s'utilitza un model de regressió lineal, el coeficient de determinació (R<math>^2</math>) informa sobre el grau d'associació.
|-
| colspan="4"| Font: D. Schwartz 1994 <ref name="Schwartz1994" />  pàg. 271.
|}

A continuació, i en el mateix ordre que el presentat a la taula, es comentarà com es realitzaran aquestes proves amb el '''R'''.

= Relació entre dues variables qualitatives =
== Dades independents: Prova de la khi quadrat ==

El fitxer de dades Aids2 del paquet MASS conté les dades de la mortalitat dels malalts diagnosticats de sida a Austràlia abans de 1991. <ref name=Aids2MASS>[http://stat.ethz.ch/R-manual/R-patched/library/MASS/html/Aids2.html Aids2, Australian AIDS Survival Data]</ref> Si es vol determinar si les dones s'han mort més que les homes, s'han de comparar la proporció de defuncions de les dones amb la dels homes <ref>S'assumeix que de tots els malalts es coneix si ha mort o no.</ref>. Les dades són:

{| class="wikitable" align=center 
|-
!Sexe !! Vius !! Morts
|-
| Dones || 36 || 53
|-
| Homes || 1046 || 1708
|}

Per avaluar l'associació estadística entre dues variables qualitatives amb una taula de contingència de ''l'' línies i ''c'' columnes mitjançant la prova de la khi quadrat, primer es calcula per cada casella els efectius esperats sota la hipòtesis de que les dues variables són independents. Els efectius esperats és calculen multiplicant el total de la fila de la casella pel total de la columna, dividit per total general de la taula. Un cop calculats els esperats de cada casella, l'estadístic χ<sup>2</sup> es calcula amb el sumatori:

:<math> \chi^{2} = \sum_{i=1}^{n} {(O_i - E_i)^2 \over E_i} ,</math>

on:

:<math>\chi^{2}</math> = és l'estadístic que s'apropa asimptòticament a una distribució ''&chi;''<sup>2</sup>;
:<math>O_i</math> = la freqüència observada;
:<math>E_i</math> = la freqüència esperada;
:<math>n </math> = el nombre de cel·les.

L'estadístic permet estimar el valor de ''p'' sota la hipòtesi de que no existeix associació entre les dues variables (la hipòtesi nul·la és certa):

* Si el valor de ''p'' és ≥ 0,05, no es pot rebutjar la hipòtesis nul·la i es conclou que no existeix associació estadística entre les dues variables.
* Si el valor de ''p'' és < 0,05, es rebutja la hipòtesi nul·la i es conclou que existeix una associació estadística entre les dues variables (pot ser real o no).

Aquest mètode només és vàlid si:

* Tots els efectius esperats calculats són com a mínim de 5. Si són més petits, es pot utilitzar la prova de Fisher.
* Les dades són independents, és a dir, els dos grups que es comparen (en l'exemple, el sexe) no estan relacionats i per cada individu d'un dels dos grups a comparar, no hi ha cap individuo sistemàticament similar en l'altre grup.

No existeix en '''R''' una funció que sigui similar al ''proc freq'' del SAS o del ''crosstabs'' del SPSS. Però el Grup de Bioestadística i Biomatemàtica de la Universitat de Lleida <ref>[http://web.udl.es/Biomath/Bioestadistica/ Grup de Bioestadística i Biomatemàtica de la Universitat de Lleida]</ref> ha modificat una funció prèvia que funciona de forma molt similar als procediments del SAS i del SPSS <ref>La funció està en un fitxer que conté tot una sèrie de funcions i que es pot trobar [http://web.udl.es/Biomath/Bioestadistica/R/Instalacio/FuncionsAuxiliars.r aquí]</ref>. S'utilitzarà per analitzar les dades anteriors. La descripció de les dades es pot fer:

<source lang="rsplus">
## ES carrega el fitxer on està la funció CrossTabs.
source("http://web.udl.es/Biomath/Bioestadistica/R/Instalacio/FuncionsAuxiliars.r", local = F, echo=TRUE, encoding = "unknown")

## Es carrega la biblioteca on estan les dades que s'analitzaran
library(MASS)

## Es visualitza la taula
CrossTabs(Aids2$sex, Aids2$status)


   Cell Contents
|-----------------|
|           Count |
|-----------------|

Total Observations in Table:  2843 

             | Aids2$status 
   Aids2$sex |        A  |        D  | Row Total | 
-------------|-----------|-----------|-----------|
           F |       36  |       53  |       89  | 
-------------|-----------|-----------|-----------|
           M |     1046  |     1708  |     2754  | 
-------------|-----------|-----------|-----------|
Column Total |     1082  |     1761  |     2843  | 
-------------|-----------|-----------|-----------|
</source>

Per realitzar la prova de la khi quadrat es pot executar el codi:

<source lang="rsplus">
CrossTabs(Aids2$sex, Aids2$status, digits = 2, chisq = T)



   Cell Contents
|-----------------|
|           Count |
|-----------------|

Total Observations in Table:  2843 

             | Aids2$status 
   Aids2$sex |        A  |        D  | Row Total | 
-------------|-----------|-----------|-----------|
           F |       36  |       53  |       89  | 
-------------|-----------|-----------|-----------|
           M |     1046  |     1708  |     2754  | 
-------------|-----------|-----------|-----------|
Column Total |     1082  |     1761  |     2843  | 
-------------|-----------|-----------|-----------|

============================================================

Statistics for All Table Factors

Pearson's Chi-squared test 
------------------------------------------------------------
Chi^2 =  0.2228163     d.f. =  1     p =  0.6369023 

Pearson's Chi-squared test with Yates' continuity correction 
------------------------------------------------------------
Chi^2 =  0.1304118     d.f. =  1     p =  0.7180055 

       Minimum expected frequency: 33.87197
</source>

El valor de ''p'' és de 0,64, major de 0,05 i, per tant, no es pot rebujar la hipotesi nul·la de que no existeix associació entre el sexe i la defunció. Es concluirà que la mortalitat és la mateixa per les dones i els homes.

Per obtenir els percentatges de morts en homes i dones, es pot executar l'ordre anterior amb l'opció de que mostri la proporció de les files:

<source lang="rsplus">
## Mostrar els percentatges de fila:
CrossTabs(Aids2$sex, Aids2$status, digits = 2, chisq = T, row =T)

   Cell Contents
|-----------------|
|           Count |
|     Row Percent |
|-----------------|

Total Observations in Table:  2843 

             | Aids2$status 
   Aids2$sex |        A  |        D  | Row Total | 
-------------|-----------|-----------|-----------|
           F |       36  |       53  |       89  | 
             |    40.45% |    59.55% |     3.13% | 
-------------|-----------|-----------|-----------|
           M |     1046  |     1708  |     2754  | 
             |    37.98% |    62.02% |    96.87% | 
-------------|-----------|-----------|-----------|
Column Total |     1082  |     1761  |     2843  | 
-------------|-----------|-----------|-----------|

============================================================

Statistics for All Table Factors

Pearson's Chi-squared test 
------------------------------------------------------------
Chi^2 =  0.2228163     d.f. =  1     p =  0.6369023 

Pearson's Chi-squared test with Yates' continuity correction 
------------------------------------------------------------
Chi^2 =  0.1304118     d.f. =  1     p =  0.7180055 

       Minimum expected frequency: 33.87197 
</source>

El percentatge de defuncions en les dones és de 60% i el dels homes de 62%. Aquesta diferència no és estadísticament signigicativa. Moltes vegades és útil conexer els efectius esperats de cada casella sota la hipòtesi nul·la. Es poden obtenir amb l'ordre:


<source lang="rsplus">
## Mostrar els esperats
CrossTabs(Aids2$sex, Aids2$status, digits = 2, chisq = T, row = T, expected = T)

   Cell Contents
|-----------------|
|           Count |
| Expected Values |
|     Row Percent |
|-----------------|

Total Observations in Table:  2843 

             | Aids2$status 
   Aids2$sex |        A  |        D  | Row Total | 
-------------|-----------|-----------|-----------|
           F |       36  |       53  |       89  | 
             |    33.87  |    55.13  |           | 
             |    40.45% |    59.55% |     3.13% | 
-------------|-----------|-----------|-----------|
           M |     1046  |     1708  |     2754  | 
             |  1048.13  |  1705.87  |           | 
             |    37.98% |    62.02% |    96.87% | 
-------------|-----------|-----------|-----------|
Column Total |     1082  |     1761  |     2843  | 
-------------|-----------|-----------|-----------|

============================================================

Statistics for All Table Factors

Pearson's Chi-squared test 
------------------------------------------------------------
Chi^2 =  0.2228163     d.f. =  1     p =  0.6369023 

Pearson's Chi-squared test with Yates' continuity correction 
------------------------------------------------------------
Chi^2 =  0.1304118     d.f. =  1     p =  0.7180055 

       Minimum expected frequency: 33.87197 
</source>

Si la hipòtesi nul·la fos certa, s'espererien 55,13 defuncions en les dones i s'observen 53. Aquesta diferència és molt petita.

Si es volguessin conèixer els percentatges de les columnes, es pot executar:

<source lang="rsplus">
## Mostrar els percentatges de columna:
CrossTabs(Aids2$sex, Aids2$status, digits = 2, chisq = T, column = T)
</source>

Si es desitga la prova de Fisher, s'obté amb l'ordre:

<source lang="rsplus">
## Prova de Fisher
CrossTabs(Aids2$sex, Aids2$status, digits = 2, chisq = T, fisher = T)
</source>

== Mesures d'associació: tau de Kendall ==

XXX http://www.r-tutor.com/gpu-computing/correlation/kendall-rank-coefficient XXX

= Relació entre una variable qualitativa i una quantitativa =

Quan es parla de la comparació d'una variable quantitativa (p. ex., l'edat) i una variable qualitativa (p. ex., el sexe), en general es vol determinar si la distribució dels valors de la variable quantitativa és la mateixa pels grups definits pels valors de la variable qualitativa (p. ex., dones i homes). (<ref>La primera versió d'aquest apartat es va treure de Barnier J. [http://alea.fr.eu.org/pages/intro-R Introduction à R]. 20 avril 2012</ref>)

En el paquet '''Epicalc''' <ref>[http://cran.r-project.org/web/packages/epicalc/index.html epicalc: Epidemiological calculator]</ref>, hi han les dades d'una enquesta (transversal) comunitària sobre la pressió arterial i els seus determinants (en aquest cas, el consum de sal). La taula de dades s'anomena BP. És un estudi transversal amb dades independents. Aquest darrer detall és important ja que determina el tipus d'anàlisi a realitzar. Les dades són aquestes:

<source lang="rsplus">
> data(BP, package="epicalc")     # Carregar les dades en memòria
> head(BP)                        # Llistar el primer registres
  id    sex sbp dbp saltadd  birthdate
1  1   male 110  80     yes 1943-01-06
2  2 female  85  55      no 1969-01-03
3  3   male 167 112     yes 1933-06-10
4  4 female 145 110     yes 1946-11-23
5  5 female 180 120      no 1941-01-03
6  6   male 112  78      no 1942-04-16
</source>

La pregunta és si les persones que solen afegir sal al que mengen (<code>saltadd = yes</code>) tenen la pressió arterial més alta.

== Estadística descriptiva ==

El primer pas seria visualitzar la pressió arterial segons el nivell de consum de sal. Es pot fer un [[Manual_de_R_(Estadística)/Gràfics#Diagrama_de_caixa | Diagrama de caixa]]:

<source lang="rsplus">
boxplot(dbp ~ saltadd, data=BP)
</source>

En el diagrama es pot veure clarament que la pressió arterial diastòlica (la "mínima") és més elevada en els que prenen sal.

Amb la funció <code>tapply</code> es poden visualitzar els estadístiques bàsiques de cada un dels dos grups: nombre de persones, mitjana, mediana o variància. Els arguments d'aquesta funció són la variable quantitativa, la qualitativa i la funció resum que es desitja: <code>tapply(variable quantitativa, variable qualitativa, funció resum)</code>:

<source lang="rsplus">
> tapply(BP$dbp, BP$saltadd, length)  # Nombre de persones
 no yes 
 37  43 

> tapply(BP$dbp, BP$saltadd, mean)    # Mitjana
       no       yes 
 88.91892 104.86047 

> tapply(BP$dbp, BP$saltadd, median)  # Mediana
 no yes 
 86 105 

> tapply(BP$dbp, BP$saltadd, var)     # Variància
      no      yes 
269.8544 638.1705 

> tapply(BP$dbp, BP$saltadd, summary)
$no
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  55,00   80,00   86,00   88,92   97,00  124,00 

$yes
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   62,0    81,0   105,0   104,9   121,0   158,0 
</source>

També amb la funció <code>tapply</code> es poden visualitzar els estadístics amb una sola ordre (<ref name="Barnier2012">Barnier J. [http://alea.fr.eu.org/pages/intro-R Introduction à R]. 20 avril 2012</ref> pàg. 81, <ref>[http://stats.stackexchange.com/questions/8225/how-to-summarize-data-by-group-in-r How to summarize data by group in R?]</ref>):

<source lang="rsplus">
> tapply(BP$dbp, BP$saltadd, function(X) c(length(X), mean(X, na.rm=TRUE), median(X, na.rm=TRUE)))
$no
[1] 37.00000 88.91892 86.00000

$yes
[1]  43.0000 104.8605 105.0000
</source>

El codi es pot millorar amb etiquetes dels estadístics:

<source lang="rsplus">
> tapply(BP$dbp, BP$saltadd, function(X){c(mitjana=mean(X, na.rm=TRUE), d.t.=sqrt(var(X, na.rm=TRUE)), n=length(X))})
$no
 mitjana     d.t.        n 
88,91892 16,42724 37,00000 

$yes
  mitjana      d.t.         n 
104,86047  25,26204  43,00000 
</source>

Ni la ordre és intuïtiva ni el el format de sortida és molt clar (es pot millorar programant). Tanmateix, es veu, per cada grup, el nombre d'observacions (37 i 43), la mitjana (89 i 105) i la mediana (86 i 105).

La diferència, tant de la mitjana com de la mediana, és molt important (entre 15 i 20 punts). És aquesta diferència estadísticament significativa? Dit d'un altre manera, si a la població d'on s'extret aquesta mostra, la pressió arterial diastòlica és la mateixa es prengui o no sal, és probable que en una mostra d'aquesta mida la  pressió arterial diastòlica dels que prenent sal sigui tant diferent de la que no en prenen? Per contestar a aquesta pregunta cal fer una prova d'hipòtesi de comparació de les dues mitjanes. Per poder determinar quina prova estadística cal (o es pot) utilitzar, primer s'ha de determinar si es versemblant suposar que la pressió arterial diastòlica en la població segueix una llei normal.

Primer es pot fer un histograma amb la corba normal superposada par estimar la forma de la distribució: <ref>[http://www.statmethods.net/graphs/density.html Robert I. Kabacoff. Histograms and Density Plots ]</ref>

<source lang="rsplus">
# Histograma amb corba normal

# Es crea una variable "x"(vector) amb els valors de la pressió pels que NO PRENENT sal
#  i s'eliminen els altres o els valors desconeguts:
x = na.omit(BP$dbp[BP$saltadd == "no"])
H = hist(x, 
         col = "red",
         xlab = "Pressió arterial diastòlica",
         ylab = "Freqüència",
         main = "Histograma amb la corba normal")              # Dibuixa l'histograma
xfit = seq(min(x),max(x),length=40)
yfit = dnorm(xfit,mean=mean(x),sd=sd(x))
yfit = yfit*diff(H$mids[1:2])*length(x)
lines(xfit, yfit, col="blue", lwd=2)                           # Dibuixa la corba normal

# Es crea una variable "x" (vector) amb els valors de la pressió pels que PRENENT sal
#  i s'eliminen els altres o els valors desconeguts:
x = na.omit(BP$dbp[BP$saltadd == "yes"])
H = hist(x, 
         col  = "red",
         xlab = "Pressió arterial diastòlica",
         ylab = "Freqüència",
         main = "Histograma amb la corba normal")              # Dibuixa l'histograma
xfit = seq(min(x),max(x),length=40)
yfit = dnorm(xfit,mean=mean(x),sd=sd(x))
yfit = yfit*diff(H$mids[1:2])*length(x)
lines(xfit, yfit, col="blue", lwd=2)                           # Dibuixa la corba normal
</source>

El codi és complicat a fi d'evitar que la corba de normal quedi tallada. L'explicació es pot trobar a llibre de Dalgaard (<ref name="Dalgaard2008">Dalgaard P. Introductory Statistics with R 2ed. New York: Springer; 2008</ref> pàg. 42-43). Si no és volgués la línia normal, seria suficient les ordres:

<source lang="rsplus">
hist(BP$dbp[BP$saltadd == "no"])
hist(BP$dbp[BP$saltadd == "yes"])
</source>

La distribució de la pressió arterial diastòlica és aproximadament normal pels entrevistats que no afegien sal, però és més dubtós pels altres.

Tanmateix, si es consulta l'ajuda del '''R''' de l'ordre <code>hist</code> (<code>help(hist)</code>), es pot llegir:

:"Comparing data with a model distribution should be done with qqplot()!"

Els diagrames q-q representen en l’eix d’ordenades els quantils de la variable estudiada i en l’eix d’abscisses els quantils teòrics o de la distribució normal. Si la variable té una distribució normal o gairebé normal, la línia de punts del diagrama q-q es troba sobre la diagonal. En canvi, en distribucions asimètriques, els punts que s’aparten de la diagonal. El diagrama de comparació quartil-quartil normal es pot fer amb les ordres:

<source lang="rsplus">
# QQPlot
qqnorm(BP$dbp[BP$saltadd == "no"], main = "Grup no sal")   # Grup no sal
qqline(BP$dbp[BP$saltadd == "no"])                         # Línia a 45 graus

qqnorm(BP$dbp[BP$saltadd == "yes"], main = "Grup no sal")  # Grup sal
qqline(BP$dbp[BP$saltadd == "yes"])                        # Línia a 45 graus
</source>

Si les dades presentessin una distribució normal, els punts estarien sobre la recta de 45°. No queda clar si aquestes dades són compatibles amb una distribució normal de la pressió a la població. Es pot fer una prova de normalitat amb l'ordre <code> shapiro.test</code>:

<source lang="rsplus">
> # Proves de normalitat
> shapiro.test(BP$dbp[BP$saltadd == "no"])      # Grup no sal

	Shapiro-Wilk normality test

data:  BP$dbp[BP$saltadd == "no"] 
W = 0.9621, p-value = 0.2349

> shapiro.test(BP$dbp[BP$saltadd == "yes"])     # Grup sal

	Shapiro-Wilk normality test

data:  BP$dbp[BP$saltadd == "yes"] 
W = 0.962, p-value = 0.1639
</source>

En cap d'elles el valor de ''p'' és menor a 0,05. Per tant no es rebutja la hipòtesi de que a la població la pressió presenta una distribució normal i, a la pràctica, s'accepta que la distribució és normal. Per tant, es podrà fer servir proves d'hipòtesi paramètriques, concretament la prova ''t'' de ''Student''.

== Prova ''t'' de Student per dades independents ==

La comparació de dues mitjanes ''m''<sub>1</sub> i ''m''<sub>2</sub> observades en dues mostres independents de mida ''n''<sub>1</sub> i ''n''<sub>2</sub> es basa en el valor de:

:<math>t ~ = ~ \frac{{\displaystyle m_{1}-m_{2}}}{{\displaystyle \sqrt{\frac{{\displaystyle s^{2}}}{n_{1}}+\frac{{\displaystyle s^{2}}}{n_{2}}}}}</math>

on ''s''<sup>2</sup> és l'estimació de la variància comuna i que s'estima:
 
:<math> s^2 = \frac{(n_1-1)s_{1}^2+(n_2-1)s_{2}^2}{n_1+n_2-2}.</math>

on ''s''<sub>1</sub><sup>2</sup> i ''s''<sub>2</sub><sup>2</sup> són els variàncies de la mostra ''1'' i ''2'', respectivament.

Si |''t''| es superior a un determinat valor crític per un risc del 5%, la diferència de les mitjanes és estadísticament significativa.

La prova només es pot utilitzar si la variable quantitativa estudiada es distribueix, en les dues poblacions d'on provenen les mostres, segons una llei normal de variància igual. Per tant, per poder-la utilitzar cal que:

* Les mostres siguin independents.
* La distribució de la variable quantitativa sigui normal pels dues grups.
* La variància sigui la mateixa. Si no ho són, existeix una opció en l'ordre de la prova ''t'' de Student per poder fer la prova amb variàncies diferents. 

Per realitzar la prova ''t'' de Student per dades independents s'utilitza l'ordre <code>t.test</code>:

<source lang="rsplus">
> t.test(<Variable quantitativa> ~ <Variable qualitativa>,
         data=<DADES>,           # Nom de la taula de dades
         paired = FALSE,         # Si FALSE: dades independents
         var.equal = TRUE,       # Si TRUE assumeix igualtat variàncies
         subset =                # Si es vol analitzar només una part de les dades           
         )
</source>

Les dades comentades en l'apartat anterior sobre el consum de sal i hipertensió, compleixin els dos primers requisits (els dos grups a comparar són independents i la distribució de la pressió arterial ja s'ha verificat que era normal mitjançant la ordre <code>shapiro.test<code>).

Com es verifica que la variància de la pressió arterial diastòlica és la mateixa pels que prenen sal i pels que no en prenen? Pes estimar la variància segons els consum de sal s'utilitzarà la ordre <code>tapply</code> i la comparació de les dues es realitza amb l'ordre <code>var.test</code>:

<source lang="rsplus">
> data(BP, package="epicalc")       # Carregen les dades en memòria
> tapply(BP$dbp, BP$saltadd, var)   # Variància de cada grup
      no      yes 
269.8544 638.1705                   # Variàncies de cada grup
> var.test(dbp ~ saltadd, data=BP)  # Prova comparació variàncies

	F test to compare two variances

data:  dbp by saltadd 
F = 0.4229, num df = 36, denom df = 42, p-value = 0.009659
alternative hypothesis: true ratio of variances is not equal to 1 
95 percent confidence interval:
 0.2249793 0.8074667 
sample estimates:
ratio of variances 
         0.4228562 
</source>

La raó de les dues variàncies està molt allunyada de 1 (0,42) i el valor de ''p'' (0,0097) és inferior a 0,05. Per tant és rebutja la hipòtesi (nul·la) de que la variància de pressió arterial diastòlica és la mateixa pels que prenen salt i que no en prenen. És molt poc probable que si les variàncies són iguals, en unes mostres d'aquesta mida la variància del grup que no pren sal sigui un 58% (0,42 - 1 = 0,58) més petita que la del grup que en pren.

Per tant la prova ''t'' de Student s'ha de realitzar amb l'opció <code>var.equal = FALSE</code>:

<source lang="rsplus">
> t.test(dbp ~ saltadd, 
+        data=BP,
+        var.equal = FALSE)	

	Welch Two Sample t-test

data:  dbp by saltadd 
t = -3.3884, df = 72.887, p-value = 0.001137
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -25.318316  -6.564777 
sample estimates:
 mean in group no mean in group yes 
         88.91892         104.86047
</source>

El valor de ''p'' és inferior a 0,05 (és 0.0011). Per tant és rebutja la hipòtesi de que la mitjana de la pressió arterial és la mateixa es prengui o no sal.

Cal tenir clar que amb un anàlisi d'aquest tipus no és possible afirmar que el consum de sal augmenta el risc d'hipertensió. Aquesta associació estadística podria ser deguda a algun biaix o a la presència d'algun altre factor de risc d'hipertensió associat al consum de sal. Al ser un estudi transversal, amb aquestes dades únicament no és possible determinar si la sal incrementa la pressió o el fet de tenir la pressió alta indueix a consumir més sal. Associació estadística no vol dir que l'associació sigui real.

=== Anàlisi d'un subgrup ===

Si es volgués només avaluar els homes, es pot fer amb l'opció <code>subset=(sex=="male")</code>:

<source lang="rsplus">
> # Seleccionar només les home
> library(epicalc)             # Es carrega la biblioteca
> tab1(BP$sex)                 # Nombre d'homes i dones
BP$sex :  
        Frequency Percent Cum. percent
male           45      45           45
female         55      55          100
  Total       100     100          100

> t.test(dbp ~ saltadd,
+        data=BP,
+        var.equal = FALSE,
+        subset=(sex=="male")
+        )

	Welch Two Sample t-test

data:  dbp by saltadd 
t = -1.8191, df = 19.484, p-value = 0.0843
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -25.927274   1.793941 
sample estimates:
 mean in group no mean in group yes 
         90.33333         102.40000
</source>

La diferència, pels homes, no és estadísticament significativa (''p'' = 0.0843), probablement per manca de potència (la mida de la mostra s'ha reduït més de la meitat).

També es poden seleccionar els homes amb l'[[Manual_de_R_(Estadística)/Modificar_les_dades#Per_indexació | indexació]]:

<source lang="rsplus">
> # Selecció del homes amb indexació
> t.test(BP$dbp[BP$sex=="male"] ~ BP$saltadd[BP$sex=="male"],
+        var.equal = FALSE
+        )

	Welch Two Sample t-test

data:  BP$dbp[BP$sex == "male"] by BP$saltadd[BP$sex == "male"] 
t = -1.8191, df = 19.484, p-value = 0.0843
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -25.927274   1.793941 
sample estimates:
 mean in group no mean in group yes 
         90.33333         102.40000
</source>

=== Publicació del resultats ===

El resultat de l'ordre <code>t.test</code> s'exporta sense cap problema major a un fitxer HTML amb el R2HTML. El codi seria:

<source lang="rsplus">
library(R2HTML)
library(epicalc)                            
library(Epi)
HTMLStart(outdir=getwd(),                   # Deixa el fitxer al directori de treball
          file="t-test",                    # Desa els resultats en el fitxer "t-test"
          extension="html",
          echo=FALSE,                       # No salvis les ordres al fitxer HTML
          HTMLframe=F)                      # No crear un fitxer HTML amb marcs
HTML.title("Relació entre consum de sal i hipertensió", HR=1)  # Títol del informe
HTML.title("Descriptiva dels dos grups (prenen i no prenen salt)", HR=3)
stat.table(index=list(saltadd), 
           contents=list(count(),percent (saltadd),mean(dbp)),
           data=BP)
HTMLhr()                                    # Crea un línia horitzontal
HTML.title("Prova t de Student de la comparació de les mitjanes", HR=3)  
t.test(dbp ~ saltadd, 
       data=BP,
       var.equal = FALSE)
HTMLStop()                                  # Desa el fitxer HTML
</source>

Es crea el fitxer t-test.html en el directori de treball que es pot compartir. Amb l'ordre <code>stat.table</code> es crea una taula amb la mitjana de la pressió arterial sistòlica per cada grup. Aquesta taula exportada a HTML queda amb un format entenedor.

Un petit problema, és que en el fitxer HTML apareixent quatre línies amb el la paraula "NULL". Sembla que és un error del paquet <ref>[https://r-forge.r-project.org/tracker/index.php?func=detail&aid=1628&group_id=405&atid=1643 Functions with no value returned add NULL to report]</ref>.

== Prova ''t'' de Student per dades aparellades ==

XXX Per fer XXX

<source lang="rsplus">
> DADES = data.frame(FARMAC  = c(5,3,5,6,4,4,7,4,3),
+                    PLACEBO = c(2,3,2,4,2,2,3,4,2))
> DADES
  FARMAC PLACEBO
1      5       2
2      3       3
3      5       2
4      6       4
5      4       2
6      4       2
7      7       3
8      4       4
9      3       2
> t.test(DADES$FARMAC, DADES$PLACEBO, paired=TRUE, alt="two.sided")

	Paired t-test

data:  DADES$FARMAC and DADES$PLACEBO 
t = 4.1538, df = 8, p-value = 0.003192
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 0.8402524 2.9375254 
sample estimates:
mean of the differences 
               1.888889 
</source>

== Prova ''U'' de Mann Whitney ==

És una prova no paramètrica, és a dir, no te cap requeriment sobre la llei de probabilitat de les variables considerades. Mentre que la prova ''t'' de Student exigeix que la llei de probabilitat de la variable quantitativa sigui normal, amb la prova ''U'' de Mann Whitney aquest requeriment desapareix.

Aquesta prova és estrictament equivalent, malgrat calcular-se de modo diferent, a la ''W'' de Wilcoxon, el que pot portar a certa confusió. Per evitar-la, en ocasions s'anomena ''prova de Mann-Whitney-Wilcoxon (<ref name="Schwartz1994" /> pàg. 257).

Per realitzar la prova ''U'' de Mann-Whitney s'utilitza l'ordre <code>wilcox.test</code>:

<source lang="rsplus">
> wilcox.test(<Variable quantitativa> ~ <Variable qualitativa>,
              data=<DADES>,           # Nom de la taula de dades
              paired = FALSE,         # Si FALSE: dades independents
              exact  = FALSE,         # Si s'ha de calcular el valor exacte de ''p''
              subset =                # Si es vol analitzar només una part de les dades           
              )
</source>

Els valors mancants els tracta segon com estigui definit a "<code>na.action<code>". Per determinar-ho es pot executar l'ordre <code>getOption("na.action")</code>. Per les dades anteriors de la relació entre consum de sal i hipertensió, l'ordre seria:

<source lang="rsplus">
> data(BP, package="epicalc")     # Es carregen les dades en memòria
> getOption("na.action")          # Opció pels valors mancants que utilitza la prova
[1] "na.omit"
> wilcox.test(dbp ~ saltadd,      # Prova
+             data = BP
+             )

	Wilcoxon rank sum test with continuity correction

data:  dbp by saltadd 
W = 523, p-value = 0.008631
alternative hypothesis: true location shift is not equal to 0 

Warning message:
In wilcox.test.default(x = c(55L, 120L, 78L, 70L, 72L, 90L, 118L,  :
  cannot compute exact p-value with ties
</source>

La prova és estadísticament significativa (valor de ''p'' =  0,0086). Per tant, es rebutja la hipòtesis nul·la de que la pressió arterial diastòlica és independent del consum de sal.

En la prova, els valors desconeguts no s'han utilitzat ja que <code>"na.action"</code> té el valor <code>"na.omit"</code>.

El '''R''' adverteix de que la prova exacte no es pot realitzar, ja que hi ha almenys dos individus amb el mateix valor de la pressió arterial diastòlica. Per evitar l'advertència es pot utilitzar l'opció "<code>exact = F </code>". Encara que no fa falta, ja que és el valor per omissió, també es pot informar l'opció <code>paired = F</code>" a fi de que quedi clar que les dades són independents:

<source lang="rsplus">
> wilcox.test(dbp ~ saltadd, 
+             data   = BP,
+             paired = F,
+             exact  = F
+             )

	Wilcoxon rank sum test with continuity correction

data:  dbp by saltadd 
W = 523, p-value = 0.008631
alternative hypothesis: true location shift is not equal to 0
</source>

El valor de ''p'' és el  mateix que abans (0,0086), però ja no apareix l'advertència sobre els empats.

=== Anàlisi d'un subgrup ===

Si es volgués només avaluar les dones, es pot fer amb l'opció <code>subset=(sex=="female")</code>:

<source lang="rsplus">
> # Seleccionar només les dones
> wilcox.test(dbp ~ saltadd, 
+             data   = BP,
+             paired = F,
+             exact  = F,
+             subset = (sex=="female")
+             )

	Wilcoxon rank sum test with continuity correction

data:  dbp by saltadd 
W = 104.5, p-value = 0.03081
alternative hypothesis: true location shift is not equal to 0 
</source>

La prova és estadísticament significativa (''p'' = 0,031).

=== Publicació del resultats ===

El resultat de l'ordre <code>wilcox.test</code> s'exporta sense cap problema major a un fitxer HTML amb el R2HTML. El codi seria:

<source lang="rsplus">
library(R2HTML)
HTMLStart(outdir=getwd(),                   # Deixa el fitxer al directori de treball
          file="u-mann-whitney",            # Desa els resultats en el fitxer "u-mann-whitney"
          extension="html",
          echo=FALSE,                       # No salvis les ordres al fitxer HTML
          HTMLframe=F)                      # No crear un fitxer HTML amb marcs
HTML.title("Relació entre consum de sal i hipertensió", HR=1)  # Títol del informe
wilcox.test(dbp ~ saltadd, 
            data   = BP,
            paired = F,
            exact  = F
            )
HTMLStop()                                  # Desa el fitxer HTML
</source>

== Comparació de diverses mitjanes: anàlisi de la variància ==

XXX Consultar:
* [http://www.uv.es/mperea/R-Anovas.html R-Anovas]
* [http://www.statmethods.net/stats/anovaAssumptions.html Assessing Classical Test Assumptions]
XXX

= Relació entre dues variables quantitatives. Regressió i correlació =

Quan s'ha d'analitzar dues variables quantitatives, la relació entre les dues pot ser:

* El valor d'una d'elles, la  variable dependent, pot dependre de l'altre, la variable independent. L'interès és estimar quin és valor mitjà de la variable dependent en funció dels valors de la variable independent. P. ex., quin és el valor mitjà de la pressió arterial diastòlica en funció de l'edat? Aquest tipus de relació s'avalua amb la ''regressió''. 
* La relació entre les dues és simètrica i les dues variables es poden considerar dependents. P. ex., la pressió arterial diastòlica de dos bessons estan relacionades? Aquest tipus de problemes s'avalua amb la correlació.

Tanmateix, en la vida real, la distinció pot no ser tan clara.

== La regressió lineal simple ==
[[File:Regressió lineal simple.svg|thumb|Relació entre la concentració d'hemoglobina en la sang i l'edat en homes adults (dades simulades de sis homes adults).]]

La relació entre dues variables quantitatives,  ''X'' i ''Y'', una de les quals (la ''Y'') pot dependre de l'altre (''X''), és independent si per qualsevol valor de la variable ''X'', la ''Y'', en mitjana, presenta el mateix valor. Per exemple, el valor de l'hemoglobina en sang és independent de l'edat en els homes adults, si en la població, sigui quina sigui l'edat, la mitjana de l'hemoglobina és la mateixa.

L'anàlisi de la relació entre una variable dependent i una independent, ambdues quantitatives, s'inicia amb la visualització  de les dades amb un [[Manual_de_R_(Estadística)/Gràfics#Diagrama_de_dispersió | diagrama de dispersió]].

En la figura del costat es mostra de la relació entre l'edat i la concentració d'hemoglobina en la sang en homes adults. Si no existeix relació entre les dues , vol dir que la concentració de l'hemoglobina (representada en les ordenades) ha de ser, en mitjana, la mateixa per qualsevol edat (representada en la abscissa). La corba que representa la variació de la mitjana de la concentració de l'hemoglobina en funció de l'edat (anomenada recta de regressió de l'hemoglobina sobre l'edat) és una recta horitzontal (és a dir, la pendent és zero).

[[File:Regressió lineal simple relació positiva.svg|thumb|Relació entre l'edat i la pressió arterial diastòlica (dades simulades de 100 homes adults).]]

Si dues variables quantitatives estan relacionades (de forma lineal), la recta de regressió de la variable dependent sobre la independent no és horitzontal i presenta una pendent diferent de zero. En la figura adjunta de la relació entre l'edat i la pressió arterial diastòlica, al augmentar el valor de l'edat, la pressió arterial diastòlica també augmenta. És una relació positiva i la pendent de la recta és més gran de zero.

Per tant, si existeix una relació (lineal) entre dues variables quantitatives, en la població la pendent de la recta de regressió de la variable dependent sobre la independent és diferent de zero. En cas contrari és zero.

[[File:Regressió lineal simple del pes del nadó i PCB.svg|thumb|Relació entre el pes del nadó i l'exposició prenatal de PCB (dades hipotètiques)]]

En la gràfica adjunta es presenten unes dades hipotètiques de la relació entre exposició prenatal a bifenils policlorats (PCB) i el pes del nadó <ref>Dades inspirades en Kezios 2012 [http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3411417/?tool=pubmed Prenatal polychlorinated biphenyl exposure is associated with decreased gestational length but not birth weight: archived samples from the Child Health and Development Studies pregnancy cohort]. Tanmateix, en aquest estudi, a diferència de l'exemple presentat aquí, la relació no era estadísticament significativa.</ref>. En aquest cas la relació és negativa (al augmentar el grau d'exposició als PCB, disminueix el pes del nado: la pendent de la recta és menor de zero). La pregunta que es planteja és si és plausible que, en la població d'on s'ha tret aquesta mostre de nadons, la variació entre el pes del nadó en funció del grau d'exposició prenatal al PCB és horitzontal (pendent de zero: hipòtesi nul·la). Si la resposta és negativa, es conclourà que existeix una relació significativa entre el pes del nadó i el grau d'exposició prenatal als PCB. El mètode de regressió lineal simple permet fer l'avaluació d'aquesta hipòtesi.

L'anàlisi de regressió és pot plantejar en els següents punts:

* És un problema de regressió i no de correlació? És a dir, existeix una variable dependent (que se sol denotar amb la lletra ''Y'') i una independent (que se sol denotar amb la lletra ''X'')?
* Les dades són independents? És una de les condicions d'aplicació de la regressió lineal simple.
* La variable dependent té una distribució normal en la població d'on s'ha tret la mostra?
* Representar les dades: diagrama dispersió: la distribució del núvol de punts segueix una línia recta (relació entre la variable dependent i la independent és lineal) o presenta altres patrons?
* La variància de la variable dependent segons els diferents valors de la variable independent és la mateixa o varia?
* Ajustar una recta a les dades per estimar
** Intercepte.
** Pendent de la recta.
* Avaluar l'ajust: La recta s'ajusta a les dades?: R<sup>2</sup>, anàlisi dels residus.
* La variable independent prediu el valor de la variable dependent? És a dir, la pendent de la recta és significativament diferent de zero?.

A continuació es desenvoluparà cada un d'aquestes etapes. Per il·lustrar-ho, es seguirà l'exemple de la relació (hipotètica) entre l'edat i la pressió arterial diastòlica presentat anteriorment (veure figura). 

Clarament és un problema de regressió, ja que existeix una variable dependent (la pressió arterial diastòlica) i una independent (l'edat). La hipòtesi és que al variar l'edat, varia la mitjana de la pressió. No te sentit el plantejament de si al variar la pressió arterial, varia l'edat.

Una de les condicions d'aplicació de la regressió lineal és que la variable dependent és distribueixi en la població segons una llei normal. Com s'ha comentat en l'aparat de l'[[#Estadística descriptiva | estadística descriptiva]] de la [[#Relació entre una variable qualitativa i una quantitativa | Relació entre una variable qualitativa i una quantitativa]], per avaluar gràficament això, es pot fer amb un diagrama de comparació quartil-quartil normal (''qqplot''):

[[File:Reg ls qq edat pad.svg|thumb|Diagrama de comparació quartil-quartil normal de la pressió arterial diastòlica (dades hipotètiques)]]

<source lang="rsplus">
# Generació de les dades hipotètiques de 100 persones
N         = 100
U         = rnorm(N)
EDAT      = 50 + rnorm(N)
PAD       = 20 + EDAT + U                     # Pressió arterial diastòlica
DADES_RLS = data.frame(PAD, EDAT)             # Es crea una taula de dades
save(DADES_RLS, file="dades/dades_rls.RData") # Se desa la taula de dades anàlisi posteriors
qqnorm(DADES_RLS$PAD)                         # Diagrama QQ normal
qqline(DADES_RLS$PAD)                         # Línia de 45º
</source>

En els valors extrems, les dades es s'allunyen de la normalitat. Es pot fer una prova de de normalitat:

<source lang="rsplus">
> shapiro.test(DADES_RLS$PAD)

	Shapiro-Wilk normality test

data:  DADES_RLS$PAD 
W = 0.9886, p-value = 0.5572
</source>

El valor de ''p'' major de 0,05 (0,557). Per tant no es pot rebutjar la hipòtesis de que en la població la variable dependent és distribueix segons una llei normal.

Tanmateix, si la distribució de la variable dependent no és normal, però la mida de la mostra és gran (p. ex., per sobre de 60), és pot utilitzar la regressió sense cap problema. <ref>Lumley T, Diehr P, Emerson S, Chen L. The importance of the normality assumption in large public health data sets. Annu. Rev Public Health 2002; 23:151–69</ref>

És molt important visualitzar el núvol de punts dels valors de la pressió arterial segons l'edat de cada persona. El diagrama de punts ja s'ha presentat anteriorment. La relació entre les dues variables és lineal i no mostre cap patró estrany.

{| class="wikitable" style="float:right;"  
|-
! Relació sinusoïdal !! Relació parabòlica
|-
| [[File:Reg ls sinusoidal.svg|thumb|Recta de regressió amb dades amb una relació entre la Y i la X sinusoïdal]] || [[File:Reg ls parabolica.svg|thumb|Recta de regressió amb dades amb una relació entre la Y i la X parabòlica]]
|}

Es pot donar el cas que la línia de regressió de ''Y'' sobre ''X'' sigui pràcticament horitzontal, suggerint una relació molt feble entre les dues variables, i ser degut al fet que la relació no sigui lineal. P, ex. en la figura, es mostra una relació sinusoïdal i un altre parabòlica. Si s'ajusta una recta de regressió a les dues, es pot arribar a concloure que les dues variables no estan relacionades, quan en realitat, en els dos casos, el valor de la ''Y'' depèn clarament del valor de la ''X''.

[[File:Reg ls curvilínia.svg|thumb|Recta de regressió amb dades amb una relació entre la Y i la X curvilínia. La línia vermella és la recta de regressió de ''Y'' sobre ''X''. La verda mostra la relació curvilínia]]

En altres ocasions, la pendent pot ser important (veure figura de dades curvilínies), però la recta de regressió no és una bona descripció de la relació entre les dues variables.

Una altra de les condicions d'aplicació de la regressió lineal, és que per cada un dels valors de la variable independent (la ''X''), la distribució de la dependent (la ''Y'') en a població és normal i amb la mateixa variància. És a dir, per les dades hipotètiques anteriors de la relació entre l'edat i la pressió arterial diastòlica, per cada edat, la distribució de la pressió en la població ha de ser normal i ha de tenir la mateixa variància sigui quina sigui l'edat.

Habitualment, es sol avaluar si es compleix aquesta assumpció després d'estimar el model de regressió. Però es pot ver una inspecció visual de si la variació del valor de la variable dependent (''Y'') és similar per qualsevol valor de la variable independent (''X''). Per exemple, en el diagrama de punts dels valors de la pressió arterial segons l'edat de cada persona presentat anteriorment, es pot observar que la dispersió dels punts de la pressió arterial per diferents valors de l'edat és similar.

[[File:Reg ls heteroscedasticitat.svg|thumb|Recta de regressió amb dades amb una relació entre la ''Y'' i la ''X'' amb variància no constant de la ''Y'' per diferents valors de la ''X'' (heteroscedasticitat)]]

En ocasions, la variació de la variable dependent és diferent segons els valors de la variable independent. P. ex., en la figura del costat ("Recta de regressió amb dades amb una relació entre la ''Y'' i la ''X'' amb variància no constant"), la dispersió dels valors de la ''Y'' sobre eixos paral·lels l'eix de les ordenades (el vertical) és més gran per valors elevats de la ''X'' que pels valors més petits de la ''X''.

El paquet '''car''' (Companion to Applied Regression) <ref>John Fox and Sanford Weisberg (2011). An {R} Companion to Applied Regression, Second Edition. Thousand Oaks CA: Sage. URL: http://socserv.socsci.mcmaster.ca/jfox/Books/Companion</ref> es troba la prova <code>ncvTest</code> que permet avaluar la hipòtesi de que les variàncies són constants. Per les dades de la relació entre edat i pressió arterial, la prova és:

<source lang="rsplus">
> library(car)
> load("dades/dades_rls.RData")    # Es carreguen les dades desades anteriorment
> ncvTest(lm(PAD ~ EDAT, data = DADES_RLS))
Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 0.03459208    Df = 1     p = 0.852453 
</source>

El valor de ''p'' és bastant més gran de 0,05 (''p'' = 0,85). Per tant, no és pot rebutjar la hipòtesi nul·la de que la variància de la pressió arterial és igual segons els valors de l'edat. No es viola la assumpció de homoscedasticitat.

=== Estimació de la recta de regressió ===

Si existeix una relació lineal entre la variable dependent i la independent (com, p. ex., en el cas la relació entre l'edat i la pressió arterial diastòlica), aquesta relació es pot descriure amb la recta:

:<math>Y = a + b \times X</math>

Aquesta és la línia de regressió, on ''a'' és l'intercepte (el valor de la ''Y'' quan la ''X'' val zero) i ''b'' la pendent de la línia.

Es pot realitzar amb l'ordre <code>lm</code>. Pes les dades anteriors de l'edat i la pressió arterial diastòlica, la recta de regressió i la representació de les dades es pot realitzar amb les ordres:

<source lang="rsplus">
> load("dades/dades_rls.RData")                        # Es carregen les dades desades anteriorment
> lm(PAD ~ EDAT, data = DADES_RLS)                     # Estimació de la regressió

Call:
lm(formula = PAD ~ EDAT, data = DADES_RLS)

Coefficients:
(Intercept)         EDAT  
    24.3954       0.9136  

> # Diagrama de dispersió
> svg("gra/reg_ls_diagr_disp_edat_pad.svg")            # Desar el diagrama en un fitxer tipus svg
> plot(PAD ~ EDAT, data = DADES_RLS, xlab= "Edat", ylab = "Pressió arterial diastòlica")
> abline(lm(PAD ~ EDAT, data = DADES_RLS),             # Afegir al diagrama la recta de regressió
+           col='red',                                 # Color vermell
+           lwd=2)                                     # Més ample del normal
> dev.off()                                            # Tancar i desar el diagrama
X11cairo 
       2 
</source>

[[File:Regressió lineal simple relació positiva.svg|thumb|Relació entre l'edat i la pressió arterial diastòlica (dades simulades de 100 homes adults). La línia vermella és la recta de regressió.]]

La sortida de l'ordre <code>lm</code> només informa del valor de l'intercepte (24.4 mmHg: el valor de la pressió arterial diastòlica quan l'edat és de zero any: no te cap sentit) i de la pendent (0,91 mmHg/any: l'increment de la pressió al augmentar l'edat en un any) de la recta de regressió.

La funció <code>coef()</code> mostra aquests coeficients. Primer s'emmagatzema el resultat de la regressió en una variable i per visualitzar els coeficients es pot utilitzar la funció <code>coef()</code>:

<source lang="rsplus">

> MODEL_RL = lm(PAD ~ EDAT, data = DADES_RLS)   # Estimació de la regressió
> coef(MODEL_RL)                                # Mostrar els coeficients de la recte de regressió
(Intercept)        EDAT 
 24.3954451   0.9135855 

</source>

Els resultats s'han emmagatzemat en la variable <code>MODEL_RL</code> que conté la informació del model de regressió.

Sovint per donar sentit a l'intercepte es realitzant un canvi d'escala de la variable independent consistent en restar, per cada un dels seus valors, la mitjana. Així els individus que tinguin com valor de ''X'' la mitjana, el valor de la nova ''X'' serà zero. I l'intercepte serà el valor de la ''Y'' quan la ''X'' pren el valor mitjà. Per l'exemple de l'edat i la pressió arterial, a l'edat de cada individuo se li restaria el valor de la mitjana de l'edat. Aquest pas previ no cal realitzar-lo ja que es pot fer directament al ajustar la regressió amb la funció <code>I()</code> que força que la resta es realitzi en el model lineal. La sintaxis és <code>I(''X'' - mean(''X''))</code> i la funció <code>I()</code> simplement vol dir que s'ha de realitzar l'operació <code>''X'' - mean(''X'')</code> al estimar la recta de regressió:

<source lang="rsplus">

> mean(DADES_RLS$EDAT)
[1] 50.03813
> lm(PAD ~ I(EDAT-mean(EDAT)), data = DADES_RLS)   # Regressió amb EDAT centrada a la mitjana

Call:
lm(formula = PAD ~ I(EDAT - mean(EDAT)), data = DADES_RLS)

Coefficients:
         (Intercept)  I(EDAT - mean(EDAT))  
             70.1095                0.9136  

</source>

La mitjana de l'edat és 50 anys. Per tant, ara l'intercepte (70,1 mmHg) correspon a la pressió arterial diastòlica d'un individu de 50 anys. Notis que la pendent (0,92 mmHg/any) no ha canviat respecte al model de regressió ajustat anteriorment amb la variable edat sense centrar. Al restar un valor a la variable independent, només varia el valor del intercepte. La pendent no canvia.

La pendent de 0,91 mmHg/any vol dir que per cada increment en un any en l'edat, la pressió arterial augmenta en 0,91 mmHg. Aquest valor és molt petit. Si es vol saber quan s'incrementa la pressió per cada augment de 10 anys en l'edat, es pot fer de dues maneres.

La primera és aprofitar la propietat de que si es multiplica la ''X'' per un factor ''A'', la pendent es divideix pel mateix factor ''A''. I viceversa: si es divideix la ''X'' per un factor ''A'', la pendent es multiplica pel mateix factor:

<source lang="rsplus">

> lm(PAD ~ I(EDAT/10), data = DADES_RLS)

Call:
lm(formula = PAD ~ I(EDAT/10), data = DADES_RLS)

Coefficients:
(Intercept)   I(EDAT/10)  
     24.395        9.136

</source>
 
Ara la pendent s'ha multiplicat per 10 i és de 9,1 mmHg/10 anys. És a dir, per cada increment en l'edat de 10 anys, la pressió arterial augmenta en 10mmHg.

La segona és simplement multiplicant el valor de la pendent pel factor que es desitgi. Si es multiplica per 5, la pendent 0,9136*5 = 4,57 correspondria al increment de la pressió arterial per un augment en l'edat de 5 anys. 

Per obtenir els intervals de confiança dels coeficients del model de regressió, es pot utilitzar l'ordre <code>confint</code>:

<source lang="rsplus">
> confint(lm(PAD ~ EDAT, data = DADES_RLS))      # Estimació IC del 95%
                 2.5 %    97.5 %
(Intercept) 14.6638432 34.127047
EDAT         0.7191395  1.108031
</source>

Per obtenir l'anàlisi de la variància del model, es pot utilitzar l'ordre <code>anova</code>:

<source lang="rsplus">

> MODEL_RL=lm(PAD ~ EDAT, data = DADES_RLS)
> anova(MODEL_RL)
Analysis of Variance Table

Response: PAD
          Df Sum Sq Mean Sq F value    Pr(>F)    
EDAT       1 81.317  81.317  86.934 3.553e-15 ***
Residuals 98 91.669   0.935                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

</source>

Per no tenir que repetir tota l'ordre <code>lm</code>, el resultat de la regressió s'emmagatzema a "<code>MODEL_RL</code>". Llavores a l'ordre <code>anova</code> només cal indicar "<code>MODEL_RL</code>". Això mateix s'hagués pogut fer en l'exemple anterior dels intervals de confiança.

=== Diagnòstics ===

Si els supòsits de regressió es compleixen, llavors els residus:

* No han de mostrar cap patró.
* S'han de distribuir normalment.
* Haurien de tenir una mitjana de zero amb variància constant <math>\sigma ^ {2}</math>.
* Han de ser independents.

Els tres darrers apartats s'ha inspirat en <ref>[http://ipsur.org/IPSUR.html#sec-Residual-Analysis-SLR 11.4 Residual Analysis en: Introduction to Probability and Statistics Using R, 11 Simple Linear Regression] (ara aquesta pàgina està desada a [http://archive.is/FUw0V#selection-1205.1-1205.25 archive.today]  (consultat 24-06-2014)). Una versió en pdf de la primera edició del llibre complet es troba [http://cran.r-project.org/web/packages/IPSUR/vignettes/IPSUR.pdf aquí] (consultat 24-06-2014).</ref> i <ref name="Quick-R: Regression Diagnostics">[http://www.statmethods.net/stats/rdiagnostics.html Quick-R: Regression Diagnostics] (consultat 24-06-2014).</ref>

==== L'assumpció de linealitat ====

La relació entre la variable dependent i la independents ha de ser lineal. El caràcter no lineal de la relació es pot posar en evidència amb un diagrama de dispersió dels valors observats contra els predits pel model de regressió, o dels residuals contra els predits o els valors observats contra el residuals. En l'exemple següent, es generen unes dades en què la relació entre les Y i la X no és lineal <ref name="Caffo-Coursera-R-Models">Exemple tret del curs Coursera [https://www.coursera.org/course/regmods ''Regression Models''] de Brian Caffo et al  consultat 29-06-2014</ref>


<source lang="rsplus">

> set.seed(123)                          # El fixa la llavor per generar nombres aleatoris                
> X = runif(100, -3, 3)                  # Valor de la variable independent
> Y = X + sin(X) + rnorm(100, sd = 0.2)  # Valor de la variable depenent
> svg("gra/reg_ls_sin2.svg")             # Desar el diagrama en un fitxer tipus svg
> plot(Y ~ X)                            # Diagrama de dispersió 
> abline(lm(Y ~ X))                      # Recta de regressió
> dev.off()                              # Tancar i desar el diagrama
X11cairo 
       2 

</source>

{| class="wikitable" style="float:right;"  
|-
! Relació no lineal (sinusoïdal) !! Residuals contra observats
|-
| [[File:Reg ls sin2.svg|thumb|Regressió lineal simple amb relació no lineal]]
 || [[File:Reg ls dgr resid x.svg|thumb|Diagrama de dispersió dels residuals contra els valors observats de la variable independent]]
|}


En el diagrama de dispersió del costat es pot veure que la relació entre la ''X'' i la ''Y'' no és lineal. El punts mostren una sinusoide al voltant de la recta de regressió.

Si la relació fos lineal, al fer un diagrama de dispersió dels valors del residuals pels diferents valors de la variable independent, els punts no haurien de mostrar cap patró i s'haurien de distribuir homogèniament per sobre i per sota del valor zero de l'eix vertical (la mitjana dels residuals és zero). En cas contrari, els punts es distribuiran seguint un patró clar.

Per posar en evidència la relació no lineal d'aquest exemple, es realitza el digrama de dispersió dels valors observats de la variable independent (les ''X'') contra els residuals:

<source lang="rsplus">

> svg("gra/reg_ls_dgr_resid_x.svg")      # Desar el diagrama en un fitxer tipus svg
> plot(resid(lm(Y ~ X)) ~ X)             # Diagrama de dispersió dels residuals sobre les X
> abline(h = 0)                          # S'afegeix una línia horitzontal a Y = 0  
> dev.off()                              # Tancar i desar el diagrama
X11cairo 
       2 

</source> 

Si la relació fos lineal, els punts es distribuirien simètricament al voltant de la línia horitzontal. Aquí es pot observar que els punts tenen un patró clarament sinusoïdal.

==== L'assumpció de normalitat dels residus ====

Es pot avaluar la normalitat dels residus amb mètodes gràfics i amb proves d'hipòtesis.

Per comprovar gràficament si els residus es distribueixen normalment se sol utilitzar un diagrama Q-Q.

Els diagrames de probabilitat normal són una eina gràfica per a la comparació d'un conjunt de dades amb la distribució normal. Es poden utilitzar amb els residus tipificats (el residu dividit per la seva desviació típica) del model de regressió lineal. <ref>[http://www.r-tutor.com/elementary-statistics/simple-linear-regression/normal-probability-plot-residuals Normal Probability Plot of Residuals]</ref>. Primer s'estima el model de regressió lineal amb l'ordre <code>lm</code> i posteriorment es calculen els residus tipificats amb l'ordre <code>rstandard</code>:

<source lang="rsplus">

# Es crea la base de dades amb nombres aleatoris:
> set.seed(1)
> N         = 100
> U         = rnorm(N)
> EDAT      = 50 + rnorm(N)
> PAD       = 20 + EDAT + U                      # Pressió arterial diastòlica
> DADES_RLS = data.frame(PAD, EDAT)              # Es crea la taula de dades

> # Regressío linal simple:
> PAD_LM = lm(PAD ~ EDAT, data = DADES_RLS)      # S'estima el model de regressió lineal
> PAD_LM_RESID_E = rstandard(PAD_LM)             # Es calculen els residu tipificats

</source>

Ara ja es pot crear el diagrama q-q amb l'ordre <code>qqnorm</code>, i s'afegeix una línia de referència amb l'ordre <code>qqline</code> per a comparació:

<source lang="rsplus">

> svg("gra/reg_ls_dgr_qq_resid_t_edat_pad.svg")  # Desar el diagrama en un fitxer tipus svg
> qqnorm(PAD_LM_RESID_E,                         # Crear el diagrama q-q
+        ylab="Residus tipificats",              # Etiqueta de l'eix de les y
+        xlab="Quantils de la llei normal",      # Etiqueta de l'eix de les x
+        main="Diagrama q-q de normalitat")      # Títol del gràfic
> qqline(PAD_LM_RESID_E)                         # S'afegeix una línia de referència
> dev.off()                                      # Tancar i desar el diagrama
X11cairo 
       2 

</source>

El gràfic s'ha desat em un fitxer tipus svg.
[[File:Reg ls dgr qq resid t edat pad.svg|thumb|Diagrama q-q dels residuals tipificats del model de regressió lineal de la pressió arterial sobre l'edat]]

Si el supòsit de normalitat és cert, llavors s'esperaria que els punts del diagrama es distribueixin a l'atzar a la proximitat de la línia recta de referència de la figura "Diagrama q-q dels residuals tipificat".

També es pot utilitzar una prova estadística per determinar si les desviacions del punts respecta a la línia de referència del diagrama q-q són estadísticament significativa, o si simplement són deguts a les fuctuacions del mostreig. Existeixen diverses proves estadístiques de normalitat. S'utilitzarà la prova de Shapiro-Wilk. Aquesta prova es realitzarà amb l'ordre <code>shapiro.test</code>. Les hipòtesi son: 
* <math>H_{0}:\mbox{ Els residus es distribueixen normalment }</math>
versus
* <math>H_{1}:\mbox{ Els residus no es distribueixen normalment }</math>

Els resultats de '''R''' són:

<source lang="rsplus">

> # Prova de Shapiro-Wilk)
> shapiro.test(residuals(PAD_LM))      # Amb els residus

	Shapiro-Wilk normality test

data:  residuals(PAD_LM) 
W = 0.9956, p-value = 0.9878

> shapiro.test(PAD_LM_RESID_E)         # Amb els residus tipificats

	Shapiro-Wilk normality test

data:  PAD_LM_RESID_E 
W = 0.9956, p-value = 0.9882

</source>

El valor de ''p'' és més gran que 0,05. Per tant no es rebutja pas la hipòtesi nul·la de que els residus es distribueixen normalment. Els residus compleixen el requeriment de normalitat.

==== L'assumpció de variància constant (homocedasticitat) ====

Es tornaran a analitzar els residuals tipificats per determinar si la variància es la mateixa per qualsevol valor de la variable independent. Es representen els residus tipificats enfront dels valors ajustats:

<source lang="rsplus">

> PAD_LM = lm(PAD ~ EDAT, data = DADES_RLS)      # S'estima el model de regressió lineal
> svg("gra/reg_ls_dgr_resid_t_ajust.svg")        # Desar el diagrama en un fitxer tipus svg
> plot(PAD_LM, which = 3)                        # Diagrama de dispersió dels residus tipificats
> dev.off()                                      # Tancar i desar el diagrama
X11cairo 
       2 
</source>

[[File:Reg ls dgr resid t ajust.svg|thumb|Residus tipificats enfront dels valors ajustats de la recta de regressió lineal de la pressió arterial sobre l'edat]]
El gràfic resultat es mostra aquí al costat ("Residus tipificats enfront dels valors ajustats de la recta de regressió lineal de la pressió arterial sobre l'edat").

Si l'assumpció de variància es compleix, en el diagrama de dispersió dels residuals contra els valors ajustats, s'esperaria que pels diferents valors de l'eix de les abscisses (l'horitzontal) la distància dels punts a la línia fos sempre similar. Aquesta distància no ha d'augmentar o disminuir al incrementar-se el valor de l'eix d'abscisses.

Es pot utilitzar la prova de Breusch-Paguen (paquet '''lmtest''') per decidir si la variància dels residus és o no constant. La hipòtesi nul·la és que la variància és la mateixa per a totes les observacions, i la hipòtesi alternativa és que la variància no és la mateixa per a totes les observacions. El resultat per la regressió de la pressió arterial sobre l'edat és:

<source lang="rsplus">

> library(lmtest)
> PAD_LM = lm(PAD ~ EDAT, data = DADES_RLS)      # S'estima el model de regressió lineal
> bptest(PAD_LM)                                 # Prova de Breusch-Pagan

	studentized Breusch-Pagan test

data:  PAD_LM 
BP = 1.982, df = 1, p-value = 0.1592

</source>

El valor de ''p'' és superior a 0,05 i per tant, no es pot rebutjar la hipòtesi nul·la de que la variància és la mateixa per totes les observacions.

Si el valor de la variància de la variable dependent és diferent per diferents valors de la variable independent, és viola una de les assumpcions per poder utilitzar un model de regressió. La violació de l'assumpció es posa en evidència al analitzar els residuals. Es veurà en un exemple. <ref name="Caffo-Coursera-R-Models" /> Es generen unes dades en que la variància de la ''Y'' depèn del valor de la ''X'' (<code>sd = 0.0001 * X)</code>):

<source lang="rsplus">
 
> set.seed(123)                                 # El fixa la llavor per generar nombres aleatoris
> X = runif(100, 0, 6)                          # Valors variable independent
> Y = X + rnorm(100, mean =0, sd = 0.0001 * X)  # Valors variable dependent
> svg("gra/reg_ls_homoc.svg")                   # Desar el diagrama en un fitxer tipus svg
> plot(Y ~ X)                                   # Diagrama de dispersió
> abline(lm(Y ~ X))                             # Línia de la recta de regressió
> dev.off()                                     # Tancar i desar el diagrama
X11cairo 
       2 

</source>

[[File:Reg ls homoc.svg|thumb|Regressió lineal simple amb violació de la homocedasticitat]]
A la figura del costat es mostra el diagrama de dispersió de les dades generades. Res fa sospitar que presenten una violació de l'assumpció de la homocedasticitat. Per evidenciar-ho es representen els residus tipificats contra els valor ajustats i els residus contre els 

<source lang="rsplus">

# Diagrama dels residus tipificats
svg("gra/reg_ls_dgr_resid_t_ajust2.svg")
plot(lm(Y ~ X), which = 3)
dev.off()

> # Diagrama dels residus sobre les X
> svg("gra/reg_ls_dgr_resid_x2.svg")
> plot(resid(lm(Y ~ X)) ~ X)
> abline(h = 0)  
> dev.off()    
X11cairo 
       2 

</source>

{| class="wikitable" style="float:right;"  
|-
! Residus tipificats !! Residus contra observats
|-
| [[File:Reg ls dgr resid t ajust2.svg|thumb|Residus tipificats enfront dels valors ajustats d'una recta de regressió lineal amb dades que violen l'assumpció d'homocedasticitat]]
 || [[File:Reg ls dgr resid x2.svg|thumb|Digrama de dispersió dels residuals contra els valors observats de la variable independent. Dades que violen l'homocedasticitat]]
|}

Als diagrames del costat és pot veure clarament que la distància dels punts a les respectives línies de referència van augmentant a mesura que augmenta el valor de les ''X''.

Al igual que amb les dades de la pressió arterial, també es pot utilitzar la prova de Breusch-Paguen (ordre <code>bptest</code>, paquet '''lmtest''') per realitzar una prova d'hipòtesi per avaluar si la variància dels residus és o no constant. Les ordres i el resultat amb '''R''' són:

<source lang="rsplus">

> # Es carrega la biblioteca
> library(lmtest)

> # Prova de Breusch-Pagan
> bptest(lm(Y ~ X))         

	studentized Breusch-Pagan test

data:  lm(Y ~ X) 
BP = 17.82, df = 1, p-value = 2.428e-05

</source>

El valor de la ''p'' és molt menor de 0,05. Per tant es rebutja la hipòtesi de que en la població la variància de la variable dependent és la mateixa per diferents valors de la variable independent.

==== L'assumpció d'independència dels residus ====

Un dels més forts dels supòsits de regressió és la relativa a la independència. La violació de la hipòtesi de la independència sovint es manifesta per la presència correlació en els residus (o autocorrelació). Pot existir correlació positiva o negativa.

Si existeix correlació positiva, els residus positius són seguits de residus positius i residus negatius seguits de residus negatius. Observant els residus d'esquerra a dreta, això és manifesta per un patró cíclic dels gràfics de residus, amb llargues seqüències de residus positius seguits per llargues seqüències del negatius.

D'altra banda, correlació negativa es posa en evidència amb residus positius seguits per residus negatius, que després són de nou seguits per residus positius, i així successivament. Per tant, els residus correlacionats negativament sovint s'associen amb un patró alternatiu positiu/negatiu en els gràfics de residus.

Les ordres per representar el diagrama del dispersió dels residus sobre els valors ajustat són:

<source lang="rsplus">

> # Regressío linal simple:
> PAD_LM = lm(PAD ~ EDAT, data = DADES_RLS)      # S'estima el model de regressió lineal
> svg("gra/reg_ls_dgr_resid_ajust.svg")          # Desar el diagrama en un fitxer tipus svg
> #  Diagrama de dispersió dels valors ajustats sobre els residus
> plot(PAD_LM, which = 1,
+      main="Diagrama de dispersió dels residus sobre els valors ajustats")
> dev.off()                                      # Tancar i desar el diagrama
X11cairo 
       2

</source>

[[File:Reg ls dgr resid ajust.svg|thumb|Diagrama de dispersió dels residus sobre els valors ajustats de la recta de regressió lineal de la pressió arterial sobre l'edat]]

El gràfic resultant es mostra a la figura del costat ("Diagrama de dispersió dels residus sobre els valors ajustats de la recta de regressió lineal de la pressió arterial sobre l'edat")

Es pot utilitzar la prova Durbin-Watson per contrastar la hipòtesi de l'existència d'autocorrelació en els residus. Es realitza amb la funció <code>dwtest()</code> des del paquet '''lmtest'''. Es durà a terme una prova de dues cues de que la correlació no és zero, que no és pas la prova predeterminat (el valor per omissió és provar que l'autocorrelació és positiva). Els resultats d'aquesta prova són:

<source lang="rsplus">

> PAD_LM = lm(PAD ~ EDAT, data = DADES_RLS)      # S'estima el model de regressió lineal
> dwtest(PAD_LM, alternative = "two.sided")      # Prova de Durbin-Watso

	Durbin-Watson test

data:  PAD_LM 
DW = 1.9962, p-value = 0.992
alternative hypothesis: true autocorrelation is not 0 

</source>

En aquest cas no es rebutja la hipòtesi nul·la; hi ha molt poca evidència d'autocorrelació diferent de zero en els residus.

==== Validació global ====

La funció <code>gvlma()</code> del paquet '''gvlma'''<ref>[http://cran.r-project.org/web/packages/gvlma/index.html gvlma: Global Validation of Linear Models Assumptions] (consultat 24-06-2014).</ref>, realitza una validació global de les hipòtesis del model lineals.<ref name="Quick-R: Regression Diagnostics" /> Els resultats de la validació per la regressió de la pressió  arterial sobre l'edat és:
 

<source lang="rsplus">

> ## Validació global
> # Es crea la base de dades amb nombres aleatoris:
> set.seed(1)
> N         = 100
> U         = rnorm(N)
> EDAT      = 50 + rnorm(N)
> PAD       = 20 + EDAT + U                      # Pressió arterial diastòlica
> DADES_RLS = data.frame(PAD, EDAT)              # Es crea la taula de dades

> # Regressío linal simple:
> PAD_LM = lm(PAD ~ EDAT, data = DADES_RLS)      # S'estima el model de regressió lineal
> library(gvlma)                                 # Es carrega el paquet gvlma
> summary(gvlma(PAD_LM))                         # Mostra els resultats de la validació

Call:
lm(formula = PAD ~ EDAT, data = DADES_RLS)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.32416 -0.60361  0.00536  0.58305  2.29316 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 20.15547    4.73336   4.258 4.73e-05 ***
EDAT         0.99907    0.09472  10.547  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.9028 on 98 degrees of freedom
Multiple R-squared: 0.5317,	Adjusted R-squared: 0.5269 
F-statistic: 111.2 on 1 and 98 DF,  p-value: < 2.2e-16 


ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance =  0.05 

Call:
 gvlma(x = PAD_LM) 

                       Value p-value                Decision
Global Stat        0.1740911  0.9964 Assumptions acceptable.
Skewness           0.0883964  0.7662 Assumptions acceptable.
Kurtosis           0.0002619  0.9871 Assumptions acceptable.
Link Function      0.0007187  0.9786 Assumptions acceptable.
Heteroscedasticity 0.0847141  0.7710 Assumptions acceptable.

</source>


==== Mesures correctives ====

Sovint es troben problemes amb el model que suggereixen que almenys es viola una de les assumpcions. Què fer llavors? S'han proposat moltes mesures, algunes de les quals s'esmenten a continuació.

* Violació de l'assumpció de normalitat dels residus. Les mateixes transformacions que es proposen a continuació per a l'estabilització de la variància poder servir per solucionar aquest problema. Sovint es maten dos ocells d'un sol tret.

* Variància no constant. De vegades, una transformació de la variable dependent es farà càrrec del problema. Hi ha una gran classe dels anomenats transformacions de Box-Cox. Prenen la forma:

:: <math>Y^{\ast}=Y^{\lambda}</math>

:on <math>Y^{\ast}</math> és la nova variable dependent calculada a partir de la variable dependent original i λ és una constant. El mètode proposat per Box i Cox determinar un valor adequat de  λ automàticament per màxima versemblança. Algunes de les transformacions del la ''Y'' mes comunes són:

<center>
{| class="wikitable"
|-
! Valor λ !! Transformació
|-
| λ = 2 || <math>Y^{\ast}=Y^{2}</math>
|-
| λ = 0,5 || <math>Y^{\ast}=\sqrt{Y}</math>
|-
| λ = 0 || <math>Y^{\ast}=\ln Y</math>
|-
| λ = -1|| <math>Y^{\ast}=1/Y</math>
|}
</center>

:Com a alternativa, es pot utilitzar el mètode dels mínims quadrats ponderats.

* No independència dels residus. Es poden utilitzar models autoregressius.

* Relació no lineal. Podem modificar directament el model per a una millor aproximació a la resposta. En particular, potser una funció de regressió polinòmica de la forma:

:: <math>Y = \beta_{0} + \beta_{1}x_{1} + \beta_{2}x_{1}^{2}</math>

=== Predicció ===

La funció <code>predict()</code> permet fer prediccions del valor de la variable dependent ''Y'' per determinats valors de la variable independent ''X''. Per exemple, per estimar els valors de la pressió arterial per les persones amb edats de 49, 50 i 51 anys, és pot fer:

<source lang="rsplus">

> DADES_EDAT = data.frame(EDAT=c(49,50,51))       # Taula de dades amb els valors de l'edat
> MODEL_RL = lm(PAD ~ EDAT, data = DADES_RLS)     # Estimació de la regressió
> predict(MODEL_RL, DADES_EDAT)                   # Predicció dels valors de la pressió arterial
       1        2        3 
69.16113 70.07472 70.98830 

</source>

La taula de dades ''DADES_EDAT'' te només una variable (EDAT) amb tres valors: 49, 50 i 51 anys. La pressió arterial predita pel model per les persones de 49, 50 i 51 anys és 69 70 i 71 mmHg respectivament.

També és podrien obtenir les prediccions utilitzant les fórmula del model de regressió:

:<math>Pressi\acute{o}  \, arterial = 24.4 + 0.914 \times Edat</math>

L'intercepte es pot obtenir amb l'ordre <code>coef(MODEL_RL)[1]</code> i la pendent amb l'ordre <code>coef(MODEL_RL)[2]</code>. Per tant, amb aquestes ordres es poden calcular els valors de la pressió arterial per les edats 49, 50 i 51 anys:

<source lang="rsplus">

> coef(MODEL_RL)[1]                               # Torna l'intercepte del model de regressió
(Intercept) 
   24.39545 
> coef(MODEL_RL)[2]                               # Torna la pendent del model de regressió 
     EDAT 
0.9135855 
> coef(MODEL_RL)[1] + coef(MODEL_RL)[2] * c(49,50,51)
[1] 69.16113 70.07472 70.98830

</source>

Les prediccions són les mateixes que les obtingudes anteriorment.

=== Resum ===

{| class="wikitable"
|-
! Resum de les ordres utilitzades en l'anàlisi de regressió lineal simple
|-
| <source lang="rsplus">
#---------------------------------------------------------------------------
# Regressió de la pressió arterial diastòlica sobre l'edat
#---------------------------------------------------------------------------
# Generació de les dades hipotètiques de 100 persones
set.seed(1)
N         = 100
U         = rnorm(N)
EDAT      = 50 + rnorm(N)
PAD       = 20 + EDAT + U                      # Pressió arterial diastòlica
DADES_RLS = data.frame(PAD, EDAT)              # Es crea una taula de dades

## Visualització dades: diagrama de dispersió de les dades
svg("gra/reg_ls_diagr_disp_edat_pad.svg")
plot(PAD ~ EDAT, data = DADES_RLS, xlab= "Edat", ylab = "Pressió arterial diastòlica")
dev.off()

## Recta de regressió de la relació Edat i pAd
lm(PAD ~ EDAT, data = DADES_RLS)               # Estimació de la regressió
# Diagrama de dispersió i recta de regressió
svg("gra/reg_ls_diagr_disp_edat_pad.svg")      # Desar el diagrama en un fitxer tipus svg
plot(PAD ~ EDAT, data = DADES_RLS, xlab= "Edat", ylab = "Pressió arterial diastòlica")
abline(lm(PAD ~ EDAT, data = DADES_RLS),       # Afegir al diagrama la recta de regressió
          col='red',                           # Color vermell
          lwd=2)                               # Més ample del normal
dev.off()                                      # Tancar i desar el diagrama

## Mostrar els coeficients, resum dels resultats i anova
MODEL_RL = lm(PAD ~ EDAT, data = DADES_RLS)    # Estimació de la regressió
coef(MODEL_RL)                                 # Mostrar els coeficients de la recte de regressió
summary(MODEL_RL)                              # Mostrar Resum dels resultats
anova(MODEL_RL)                                # Mostrar l'ANOVA

## I de C del 95%
confint(lm(PAD ~ EDAT, data = DADES_RLS))      # Estimació IC del 95%

## Diagnòstic
# Normalitat dels residus
# Diagrama q-q
MODEL_RL_RESID_E = rstandard(MODEL_RL)         # Es calculen els residu tipificats
svg("gra/reg_ls_dgr_qq_resid_t_edat_pad.svg")  # Desar el diagrama en un fitxer tipus svg
qqnorm(MODEL_RL_RESID_E,                       # Crear el diagrama q-q
       ylab="Residus tipificats",              # Etiqueta de l'eix de les y
       xlab="Quantils de la llei normal",      # Etiqueta de l'eix de les x
       main="Diagrama q-q de normalitat")      # Títol del gràfic
qqline(MODEL_RL_RESID_E)                       # S'afegeix una línia de referència
dev.off()                                      # Tancar i desar el diagrama
# Prova de Shapiro-Wilk)
shapiro.test(residuals(MODEL_RL))              # Amb els residus
shapiro.test(MODEL_RL_RESID_E)                 # Amb els residus tipificats

# Variància constant
# Diagrama dels residus tipificats
svg("gra/reg_ls_dgr_resid_t_ajust.svg")        # Desar el diagrama en un fitxer tipus svg
plot(MODEL_RL, which = 3)                      # Diagrama de dispersió dels residus tipificats
dev.off()     
# Prova d'hipòtesi
library(lmtest)
MODEL_RL = lm(PAD ~ EDAT, data = DADES_RLS)    # S'estima el model de regressió lineal
bptest(MODEL_RL)                               # Prova de Breusch-Pagan

## Independència
#  Diagrama de dispersió dels valors ajustats sobre els residus
svg("gra/reg_ls_dgr_resid_ajust.svg")          # Desar el diagrama en un fitxer tipus svg
plot(MODEL_RL, which = 1,
     main="Diagrama de dispersió dels residus sobre els valors ajustats")
dev.off()                                      # Tancar i desar el diagrama
# Prova d'hipòtesi
dwtest(MODEL_RL, alternative = "two.sided")    # Prova de Durbin-Watson

## Validació global
library(gvlma)                                 # Es carrega el paquet gvlma
summary(gvlma(MODEL_RL))                       # Mostra els resultats de la validació

</source>
|}

== Correlació ==
=== Coeficient de Correlació de Pearson ===

XXX
http://ww2.coastal.edu/kingw/statistics/R-tutorials/simplelinear.html
There is also a formula interface for cor.test( ), but it's tricky. Both variables should be listed after the tilde...

> with(cats, cor.test(~ Bwt + Hwt))    # output not shown

Using the formula interface makes it easy to subset the data by rows of the data frame...

> with(cats, cor.test(~ Bwt + Hwt, subset=(Sex=="F")))

The "subset=" option is not available unless you use the formula interface.

cor(DADES[3:8]) # Columns 3 to 8 contain the 6 baseline measures

round(cor(DADES[3:8]), 2) # Round to 2 decimal places

round(cor(DADES$S1.pre, DADES$S2.pre), 2) # Round to 2 decimal places

cor.test(DADES$EDAT, PE$TA)
XXX
[[Image:Anscombe's quartet 3.svg|thumb|325px|right|Quatre conjunt de dades amb la mateixa correlació de 0,816.]]

La imatge de la dreta mostra quatre diagrames de dispersió de un conjunt de quatre diferents parells de variables creades per Francis Anscombe <ref>{{cite journal | last=Anscombe | first=Francis J. | year=1973 | title=Graphs in statistical analysis | journal=The American Statistician | volume=27 | pages=17–21 | jstor=2682899 | doi=10.2307/2682899}}</ref>. Les quatre variables Y tenen la mateixa mitjana (7,5), variància (4,12), correlació (0,816) i regressió la lineal (''y'' = 3 + 0,5 ''x''). No obstant això, com es pot veure en els diagrames de dispersió, la distribució de les variables és molt diferent. La primera (a dalt a l'esquerra) sembla que es distribueix normalment, i es correspon amb el que es podria esperar quan les dues variables considerades estan correlacionades i seguint el supòsit de normalitat. El segon (dalt a la dreta) no es distribueix normalment, mentre que existeix una relació òbvia entre les dues variables que no és lineal. En aquest cas, el coeficient de correlació de Pearson no indica que hi ha una relació funcional exacta: només la mesura que aquesta relació es pot aproximar per una relació lineal. En el tercer cas (part inferior esquerra), la relació lineal és perfecta, a excepció d'un valor atípic, que exerceix una influència suficient per reduir el coeficient de correlació 1-,816. Finalment, el quart exemple (baix a la dreta) mostra un altre exemple, en el que un valor atípic és suficient per produir un alt coeficient de correlació, tot i que la relació entre les dues variables no és lineal.<ref>Traduït de: [https://en.wikipedia.org/wiki/Correlation_and_dependence#Correlation_and_linearity Correlation and linearity]</ref>

=== Coeficient de correlació d’Spearman (no paramètric) ===

XXX
:cor.test(X1, X2, use = "complete.obs", method = "spearman")
XXX

= Referències =
<references/>

= Enllaços externs =
== Regressió lineal simple ==
* [http://www.uam.es/personal_pdi/ciencias/joser/paginaR/regresion.html Regresión lineal simple con R] {{es}}
* [http://zoonek2.free.fr/UNIX/48_R/09.html Regression] per Vincent Zoonekynd. {{en}}
* [http://zoonek2.free.fr/UNIX/48_R/10.html Other regressions] per Vincent Zoonekynd. {{en}}
* [http://zoonek2.free.fr/UNIX/48_R/11.html Regression Problems -- and their Solutions] per Vincent Zoonekynd. {{en}}
[[Categoria:Manual de R (Estadística)]]