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 ([1] pàg. 271):
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) informa sobre el grau d'associació.
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. [2] 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 [3]. Les dades són:
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 χ2 es calcula amb el sumatori:
on:
= és l'estadístic que s'apropa asimptòticament a una distribució χ2;
= la freqüència observada;
= la freqüència esperada;
= 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 [4] ha modificat una funció prèvia que funciona de forma molt similar als procediments del SAS i del SPSS [5]. S'utilitzarà per analitzar les dades anteriors. La descripció de les dades es pot fer:
## 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'analitzaranlibrary(MASS)## Es visualitza la taulaCrossTabs(Aids2$sex,Aids2$status)CellContents|-----------------||Count||-----------------|TotalObservationsinTable:2843|Aids2$statusAids2$sex|A|D|RowTotal|-------------|-----------|-----------|-----------|F|36|53|89|-------------|-----------|-----------|-----------|M|1046|1708|2754|-------------|-----------|-----------|-----------|ColumnTotal|1082|1761|2843|-------------|-----------|-----------|-----------|
Per realitzar la prova de la khi quadrat es pot executar el codi:
CrossTabs(Aids2$sex,Aids2$status,digits=2,chisq=T)CellContents|-----------------||Count||-----------------|TotalObservationsinTable:2843|Aids2$statusAids2$sex|A|D|RowTotal|-------------|-----------|-----------|-----------|F|36|53|89|-------------|-----------|-----------|-----------|M|1046|1708|2754|-------------|-----------|-----------|-----------|ColumnTotal|1082|1761|2843|-------------|-----------|-----------|-----------|============================================================StatisticsforAllTableFactorsPearson's Chi-squared test ------------------------------------------------------------Chi^2 = 0.2228163 d.f. = 1 p = 0.6369023 Pearson'sChi-squaredtestwithYates' continuity correction ------------------------------------------------------------Chi^2=0.1304118d.f.=1p=0.7180055Minimumexpectedfrequency:33.87197
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:
## Mostrar els percentatges de fila:CrossTabs(Aids2$sex,Aids2$status,digits=2,chisq=T,row=T)CellContents|-----------------||Count||RowPercent||-----------------|TotalObservationsinTable:2843|Aids2$statusAids2$sex|A|D|RowTotal|-------------|-----------|-----------|-----------|F|36|53|89||40.45% | 59.55%|3.13% | -------------|-----------|-----------|-----------| M | 1046 | 1708 | 2754 | | 37.98%|62.02% | 96.87%|-------------|-----------|-----------|-----------|ColumnTotal|1082|1761|2843|-------------|-----------|-----------|-----------|============================================================StatisticsforAllTableFactorsPearson's Chi-squared test ------------------------------------------------------------Chi^2 = 0.2228163 d.f. = 1 p = 0.6369023 Pearson'sChi-squaredtestwithYates' continuity correction ------------------------------------------------------------Chi^2=0.1304118d.f.=1p=0.7180055Minimumexpectedfrequency:33.87197
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:
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). ([6])
En el paquet Epicalc[7], 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:
>data(BP,package="epicalc")# Carregar les dades en memòria>head(BP)# Llistar el primer registresidsexsbpdbpsaltaddbirthdate11male11080yes1943-01-0622female8555no1969-01-0333male167112yes1933-06-1044female145110yes1946-11-2355female180120no1941-01-0366male11278no1942-04-16
La pregunta és si les persones que solen afegir sal al que mengen (saltadd = yes) tenen la pressió arterial més alta.
El primer pas seria visualitzar la pressió arterial segons el nivell de consum de sal. Es pot fer un Diagrama de caixa:
boxplot(dbp~saltadd,data=BP)
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ó tapply 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: tapply(variable quantitativa, variable qualitativa, funció resum):
>tapply(BP$dbp,BP$saltadd,length)# Nombre de personesnoyes3743>tapply(BP$dbp,BP$saltadd,mean)# Mitjananoyes88.91892104.86047>tapply(BP$dbp,BP$saltadd,median)# Mediananoyes86105>tapply(BP$dbp,BP$saltadd,var)# Variàncianoyes269.8544638.1705>tapply(BP$dbp,BP$saltadd,summary)$noMin.1stQu.MedianMean3rdQu.Max.55,0080,0086,0088,9297,00124,00$yesMin.1stQu.MedianMean3rdQu.Max.62,081,0105,0104,9121,0158,0
També amb la funció tapply es poden visualitzar els estadístics amb una sola ordre ([8] pàg. 81, [9]):
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ó: [10]
# 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'histogramaxfit=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'histogramaxfit=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
El codi és complicat a fi d'evitar que la corba de normal quedi tallada. L'explicació es pot trobar a llibre de Dalgaard ([11] pàg. 42-43). Si no és volgués la línia normal, seria suficient les ordres:
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 hist (help(hist)), 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:
# QQPlotqqnorm(BP$dbp[BP$saltadd=="no"],main="Grup no sal")# Grup no salqqline(BP$dbp[BP$saltadd=="no"])# Línia a 45 grausqqnorm(BP$dbp[BP$saltadd=="yes"],main="Grup no sal")# Grup salqqline(BP$dbp[BP$saltadd=="yes"])# Línia a 45 graus
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 shapiro.test:
># Proves de normalitat>shapiro.test(BP$dbp[BP$saltadd=="no"])# Grup no salShapiro-Wilknormalitytestdata:BP$dbp[BP$saltadd=="no"]W=0.9621,p-value=0.2349>shapiro.test(BP$dbp[BP$saltadd=="yes"])# Grup salShapiro-Wilknormalitytestdata:BP$dbp[BP$saltadd=="yes"]W=0.962,p-value=0.1639
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.
La comparació de dues mitjanes m1 i m2 observades en dues mostres independents de mida n1 i n2 es basa en el valor de:
on s2 és l'estimació de la variància comuna i que s'estima:
on s12 i s22 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 t.test:
>t.test(<Variablequantitativa>~<Variablequalitativa>,data=<DADES>,# Nom de la taula de dadespaired=FALSE,# Si FALSE: dades independentsvar.equal=TRUE,# Si TRUE assumeix igualtat variànciessubset=# Si es vol analitzar només una part de les dades )
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 shapiro.test).
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 tapply i la comparació de les dues es realitza amb l'ordre var.test:
>data(BP,package="epicalc")# Carregen les dades en memòria>tapply(BP$dbp,BP$saltadd,var)# Variància de cada grupnoyes269.8544638.1705# Variàncies de cada grup>var.test(dbp~saltadd,data=BP)# Prova comparació variànciesFtesttocomparetwovariancesdata:dbpbysaltaddF=0.4229,numdf=36,denomdf=42,p-value=0.009659alternativehypothesis:trueratioofvariancesisnotequalto195percentconfidenceinterval:0.22497930.8074667sampleestimates:ratioofvariances0.4228562
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ó var.equal = FALSE:
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.
Si es volgués només avaluar els homes, es pot fer amb l'opció subset=(sex=="male"):
># Seleccionar només les home>library(epicalc)# Es carrega la biblioteca>tab1(BP$sex)# Nombre d'homes i donesBP$sex:FrequencyPercentCum.percentmale454545female5555100Total100100100>t.test(dbp~saltadd,+data=BP,+var.equal=FALSE,+subset=(sex=="male")+)WelchTwoSamplet-testdata:dbpbysaltaddt=-1.8191,df=19.484,p-value=0.0843alternativehypothesis:truedifferenceinmeansisnotequalto095percentconfidenceinterval:-25.9272741.793941sampleestimates:meaningroupnomeaningroupyes90.33333102.40000
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' indexació:
># Selecció del homes amb indexació>t.test(BP$dbp[BP$sex=="male"]~BP$saltadd[BP$sex=="male"],+var.equal=FALSE+)WelchTwoSamplet-testdata:BP$dbp[BP$sex=="male"]byBP$saltadd[BP$sex=="male"]t=-1.8191,df=19.484,p-value=0.0843alternativehypothesis:truedifferenceinmeansisnotequalto095percentconfidenceinterval:-25.9272741.793941sampleestimates:meaningroupnomeaningroupyes90.33333102.40000
El resultat de l'ordre t.test s'exporta sense cap problema major a un fitxer HTML amb el R2HTML. El codi seria:
library(R2HTML)library(epicalc)library(Epi)HTMLStart(outdir=getwd(),# Deixa el fitxer al directori de treballfile="t-test",# Desa els resultats en el fitxer "t-test"extension="html",echo=FALSE,# No salvis les ordres al fitxer HTMLHTMLframe=F)# No crear un fitxer HTML amb marcsHTML.title("Relació entre consum de sal i hipertensió",HR=1)# Títol del informeHTML.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 horitzontalHTML.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
Es crea el fitxer t-test.html en el directori de treball que es pot compartir. Amb l'ordre stat.table 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 [12].
É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 ([1] pàg. 257).
Per realitzar la prova U de Mann-Whitney s'utilitza l'ordre wilcox.test:
>wilcox.test(<Variablequantitativa>~<Variablequalitativa>,data=<DADES>,# Nom de la taula de dadespaired=FALSE,# Si FALSE: dades independentsexact=FALSE,# Si s'ha de calcular el valor exacte de ''p''subset=# Si es vol analitzar només una part de les dades )
Els valors mancants els tracta segon com estigui definit a "na.action". Per determinar-ho es pot executar l'ordre getOption("na.action"). Per les dades anteriors de la relació entre consum de sal i hipertensió, l'ordre seria:
>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+)Wilcoxonranksumtestwithcontinuitycorrectiondata:dbpbysaltaddW=523,p-value=0.008631alternativehypothesis:truelocationshiftisnotequalto0Warningmessage:Inwilcox.test.default(x=c(55L,120L,78L,70L,72L,90L,118L,:cannotcomputeexactp-valuewithties
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 "na.action" té el valor "na.omit".
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ó "exact = F ". Encara que no fa falta, ja que és el valor per omissió, també es pot informar l'opció paired = F" a fi de que quedi clar que les dades són independents:
Si es volgués només avaluar les dones, es pot fer amb l'opció subset=(sex=="female"):
># Seleccionar només les dones>wilcox.test(dbp~saltadd,+data=BP,+paired=F,+exact=F,+subset=(sex=="female")+)Wilcoxonranksumtestwithcontinuitycorrectiondata:dbpbysaltaddW=104.5,p-value=0.03081alternativehypothesis:truelocationshiftisnotequalto0
La prova és estadísticament significativa (p = 0,031).
El resultat de l'ordre wilcox.test s'exporta sense cap problema major a un fitxer HTML amb el R2HTML. El codi seria:
library(R2HTML)HTMLStart(outdir=getwd(),# Deixa el fitxer al directori de treballfile="u-mann-whitney",# Desa els resultats en el fitxer "u-mann-whitney"extension="html",echo=FALSE,# No salvis les ordres al fitxer HTMLHTMLframe=F)# No crear un fitxer HTML amb marcsHTML.title("Relació entre consum de sal i hipertensió",HR=1)# Títol del informewilcox.test(dbp~saltadd,data=BP,paired=F,exact=F)HTMLStop()# Desa el fitxer HTML
Comparació de diverses mitjanes: anàlisi de la variància
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.
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 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).
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.
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ó [13]. 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?: R2, 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 de la 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):
Diagrama de comparació quartil-quartil normal de la pressió arterial diastòlica (dades hipotètiques)
# Generació de les dades hipotètiques de 100 personesN=100U=rnorm(N)EDAT=50+rnorm(N)PAD=20+EDAT+U# Pressió arterial diastòlicaDADES_RLS=data.frame(PAD,EDAT)# Es crea una taula de dadessave(DADES_RLS,file="dades/dades_rls.RData")# Se desa la taula de dades anàlisi posteriorsqqnorm(DADES_RLS$PAD)# Diagrama QQ normalqqline(DADES_RLS$PAD)# Línia de 45º
En els valors extrems, les dades es s'allunyen de la normalitat. Es pot fer una prova de de normalitat:
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. [14]
É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.
Relació sinusoïdal
Relació parabòlica
Recta de regressió amb dades amb una relació entre la Y i la X sinusoïdal
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.
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.
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) [15] es troba la prova ncvTest 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:
>library(car)>load("dades/dades_rls.RData")# Es carreguen les dades desades anteriorment>ncvTest(lm(PAD~EDAT,data=DADES_RLS))Non-constantVarianceScoreTestVarianceformula:~fitted.valuesChisquare=0.03459208Df=1p=0.852453
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.
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:
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 lm. 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:
>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)EDAT24.39540.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 diagramaX11cairo2
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 lm 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ó coef() mostra aquests coeficients. Primer s'emmagatzema el resultat de la regressió en una variable i per visualitzar els coeficients es pot utilitzar la funció coef():
>MODEL_RL=lm(PAD~EDAT,data=DADES_RLS)# Estimació de la regressió>coef(MODEL_RL)# Mostrar els coeficients de la recte de regressió(Intercept)EDAT24.39544510.9135855
Els resultats s'han emmagatzemat en la variable MODEL_RL 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ó I() que força que la resta es realitzi en el model lineal. La sintaxis és I(X - mean(X)) i la funció I() simplement vol dir que s'ha de realitzar l'operació X - mean(X) al estimar la recta de regressió:
>mean(DADES_RLS$EDAT)[1]50.03813>lm(PAD~I(EDAT-mean(EDAT)),data=DADES_RLS)# Regressió amb EDAT centrada a la mitjanaCall:lm(formula=PAD~I(EDAT-mean(EDAT)),data=DADES_RLS)Coefficients:(Intercept)I(EDAT-mean(EDAT))70.10950.9136
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:
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 confint:
>confint(lm(PAD~EDAT,data=DADES_RLS))# Estimació IC del 95%2.5% 97.5 %(Intercept)14.663843234.127047EDAT0.71913951.108031
Per obtenir l'anàlisi de la variància del model, es pot utilitzar l'ordre anova:
Per no tenir que repetir tota l'ordre lm, el resultat de la regressió s'emmagatzema a "MODEL_RL". Llavores a l'ordre anova només cal indicar "MODEL_RL". Això mateix s'hagués pogut fer en l'exemple anterior dels intervals de confiança.
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 [18]
>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 diagramaX11cairo2
Relació no lineal (sinusoïdal)
Residuals contra observats
Regressió lineal simple amb relació no lineal
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:
>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 diagramaX11cairo2
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.
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. [19]. Primer s'estima el model de regressió lineal amb l'ordre lm i posteriorment es calculen els residus tipificats amb l'ordre rstandard:
# 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
Ara ja es pot crear el diagrama q-q amb l'ordre qqnorm, i s'afegeix una línia de referència amb l'ordre qqline per a comparació:
>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 diagramaX11cairo2
El gràfic s'ha desat em un fitxer tipus svg.
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 shapiro.test. Les hipòtesi son:
versus
Els resultats de R són:
># Prova de Shapiro-Wilk)>shapiro.test(residuals(PAD_LM))# Amb els residusShapiro-Wilknormalitytestdata:residuals(PAD_LM)W=0.9956,p-value=0.9878>shapiro.test(PAD_LM_RESID_E)# Amb els residus tipificatsShapiro-Wilknormalitytestdata:PAD_LM_RESID_EW=0.9956,p-value=0.9882
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:
>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 diagramaX11cairo2
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:
>library(lmtest)>PAD_LM=lm(PAD~EDAT,data=DADES_RLS)# S'estima el model de regressió lineal>bptest(PAD_LM)# Prova de Breusch-PaganstudentizedBreusch-Pagantestdata:PAD_LMBP=1.982,df=1,p-value=0.1592
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. [18] Es generen unes dades en que la variància de la Y depèn del valor de la X (sd = 0.0001 * X)):
>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 diagramaX11cairo2
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
# Diagrama dels residus tipificatssvg("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()X11cairo2
Residus tipificats
Residus contra observats
Residus tipificats enfront dels valors ajustats d'una recta de regressió lineal amb dades que violen l'assumpció d'homocedasticitat
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 bptest, 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:
># Es carrega la biblioteca>library(lmtest)># Prova de Breusch-Pagan>bptest(lm(Y~X))studentizedBreusch-Pagantestdata:lm(Y~X)BP=17.82,df=1,p-value=2.428e-05
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.
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:
># 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 diagramaX11cairo2
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ó dwtest() 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:
>PAD_LM=lm(PAD~EDAT,data=DADES_RLS)# S'estima el model de regressió lineal>dwtest(PAD_LM,alternative="two.sided")# Prova de Durbin-WatsoDurbin-Watsontestdata:PAD_LMDW=1.9962,p-value=0.992alternativehypothesis:trueautocorrelationisnot0
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.
La funció gvlma() del paquet gvlma[20], realitza una validació global de les hipòtesis del model lineals.[17] Els resultats de la validació per la regressió de la pressió arterial sobre l'edat és:
>## 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:Min1QMedian3QMax-2.32416-0.603610.005360.583052.29316Coefficients:EstimateStd.ErrortvaluePr(>|t|)(Intercept)20.155474.733364.2584.73e-05***EDAT0.999070.0947210.547<2e-16***---Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1Residualstandarderror:0.9028on98degreesoffreedomMultipleR-squared:0.5317,AdjustedR-squared:0.5269F-statistic:111.2on1and98DF,p-value:<2.2e-16ASSESSMENTOFTHELINEARMODELASSUMPTIONSUSINGTHEGLOBALTESTON4DEGREES-OF-FREEDOM:LevelofSignificance=0.05Call:gvlma(x=PAD_LM)Valuep-valueDecisionGlobalStat0.17409110.9964Assumptionsacceptable.Skewness0.08839640.7662Assumptionsacceptable.Kurtosis0.00026190.9871Assumptionsacceptable.LinkFunction0.00071870.9786Assumptionsacceptable.Heteroscedasticity0.08471410.7710Assumptionsacceptable.
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:
on é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:
Valor λ
Transformació
λ = 2
λ = 0,5
λ = 0
λ = -1
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:
La funció predict() 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:
>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ó arterial12369.1611370.0747270.98830
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ó:
L'intercepte es pot obtenir amb l'ordre coef(MODEL_RL)[1] i la pendent amb l'ordre coef(MODEL_RL)[2]. Per tant, amb aquestes ordres es poden calcular els valors de la pressió arterial per les edats 49, 50 i 51 anys:
>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ó EDAT0.9135855>coef(MODEL_RL)[1]+coef(MODEL_RL)[2]*c(49,50,51)[1]69.1611370.0747270.98830
Les prediccions són les mateixes que les obtingudes anteriorment.
Resum de les ordres utilitzades en l'anàlisi de regressió lineal simple
#---------------------------------------------------------------------------# Regressió de la pressió arterial diastòlica sobre l'edat#---------------------------------------------------------------------------# Generació de les dades hipotètiques de 100 personesset.seed(1)N=100U=rnorm(N)EDAT=50+rnorm(N)PAD=20+EDAT+U# Pressió arterial diastòlicaDADES_RLS=data.frame(PAD,EDAT)# Es crea una taula de dades## Visualització dades: diagrama de dispersió de les dadessvg("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 pAdlm(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 svgplot(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 vermelllwd=2)# Més ample del normaldev.off()# Tancar i desar el diagrama## Mostrar els coeficients, resum dels resultats i anovaMODEL_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 resultatsanova(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-qMODEL_RL_RESID_E=rstandard(MODEL_RL)# Es calculen els residu tipificatssvg("gra/reg_ls_dgr_qq_resid_t_edat_pad.svg")# Desar el diagrama en un fitxer tipus svgqqnorm(MODEL_RL_RESID_E,# Crear el diagrama q-qylab="Residus tipificats",# Etiqueta de l'eix de les yxlab="Quantils de la llei normal",# Etiqueta de l'eix de les xmain="Diagrama q-q de normalitat")# Títol del gràficqqline(MODEL_RL_RESID_E)# S'afegeix una línia de referènciadev.off()# Tancar i desar el diagrama# Prova de Shapiro-Wilk)shapiro.test(residuals(MODEL_RL))# Amb els residusshapiro.test(MODEL_RL_RESID_E)# Amb els residus tipificats# Variància constant# Diagrama dels residus tipificatssvg("gra/reg_ls_dgr_resid_t_ajust.svg")# Desar el diagrama en un fitxer tipus svgplot(MODEL_RL,which=3)# Diagrama de dispersió dels residus tipificatsdev.off()# Prova d'hipòtesilibrary(lmtest)MODEL_RL=lm(PAD~EDAT,data=DADES_RLS)# S'estima el model de regressió linealbptest(MODEL_RL)# Prova de Breusch-Pagan## Independència# Diagrama de dispersió dels valors ajustats sobre els residussvg("gra/reg_ls_dgr_resid_ajust.svg")# Desar el diagrama en un fitxer tipus svgplot(MODEL_RL,which=1,main="Diagrama de dispersió dels residus sobre els valors ajustats")dev.off()# Tancar i desar el diagrama# Prova d'hipòtesidwtest(MODEL_RL,alternative="two.sided")# Prova de Durbin-Watson## Validació globallibrary(gvlma)# Es carrega el paquet gvlmasummary(gvlma(MODEL_RL))# Mostra els resultats de la validació
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
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 [21]. 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.[22]
Coeficient de correlació d’Spearman (no paramètric)
↑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
