library(mnormt) # generació de normals 'multivariants'
library(car)
library(rgl)
library(Rcmdr)
library(grDevices)# recursos gràfics, entre altres
library(xtable) # taules amb LaTeX
library(fBasics) # funcions matricials
set.seed(666)

1. Simulació de dades multivariants

1.a. Esfericitat. Genera una mostra Gaussiana multivariant data1 \(\sim \mathrm{MN}(\mu, \sum)\) de mida \(10^{3}\), seguint els passos següents:

  • Guarda el vector de mitjanes \(\mu_{1}\) i la matriu de covariàncies teòriques \(\sum_{1}\) com a objectes de R: \[\mu_{1} = \left( \begin{array}{c} -3 \\ 1 \\ 13 \end{array} \right) \qquad \sum\nolimits_{1} = \left( \begin{array}{ccc}4 & 0 & 0 \\ 0 & 4 & 0 \\ 0 & 0 & 4 \end{array} \right) \]
mu1 <- c(-3, 1, 13)
sigma1 <- matrix(data = c(4, 0, 0, 0, 4, 0, 0, 0, 4), nrow = 3, byrow = T)
nsample <- 10^3
  • Genera una mostra Gaussiana multivariant amb la funció rmnorm() del paquet mnormt i assigna-la a l’objecte de R que anomenaràs data1 (ha de tenir dimensions \(10^{3} \times 3\)). Comprova que data1 és una matriu i converteix-la a data frame (les funcions as.data.frame(), as.matrix(), etc. canvien el tipus de l’objecte quan és possible). Assigna els noms X1, X2 i X3 a les columnes de data1. Mostra els 5 primers casos del fitxer data1.
data1 <- rmnorm(n = nsample, mean = mu1, varcov = sigma1)
dim(data1) # Ha de ser 10^3 * 3
class(data1) # Comprovem que és una matriu
data1 <- as.data.frame(data1) # Covertim l'objecte data1 en data.frame
colnames(data1) <- c("X1", "X2", "X3") # Canviem el nom de les columnes
data1[1:5,] # Mostrem els 5 primers casos
## [1] 1000    3
## [1] "matrix"
  • Per què parlem d’esfericitat?

Perquè ens trobem en un cas de multimostres, ja que tenim 3 mostres diferents (cada columna del objecte dades1 fa referència a una mostra) i hi ha homogeneïtat en les matrius de covariàncies, és a dir, les covariàncies són 0 (incorrelació) i les variàncies són iguals.

1.b. El·lipticitat i angles amb els eixos. Fes el mateix per a una nova mostra Gaussiana data2, amb vector de mitjanes i matriu de covariàncies: \[\mu_{2} = \left( \begin{array}{c} -7 \\ 2 \\ 10 \end{array} \right) \qquad \sum\nolimits_{2} = \left( \begin{array}{c} 3 & 1.385 & -2.939 \\ 1.385 & 1 & -1.979 \\ -2.939 & -1.979 & 8 \end{array} \right) \]

mu2 <- c(-7, 2, -10)
sigma2 <- matrix(data = c(3, 1.385, -2.939, 1.385, 1, -1.979, -2.939, -1.979, 8), nrow = 3, byrow = T)
data2 <- rmnorm(n = nsample, mean = mu2, varcov = sigma2)

Converteix data2 a data.frame amb els mateixos noms de columna que abans, X1, X2, i X3. Mostra els 5 últims casos del fitxer data2.

data2 <- as.data.frame(data2) # Covertim l'objecte en data.frame
colnames(data2) <-  paste0("X", 1:3) # Canviem el nom de les columnes
data2[c(nsample-5):nsample,] # Mostrem els 5 ultims casos

1.c. Per què parlem d’el·lipticitat i angles amb els eixos? Com seria un cas amb forma el·líptica sense angles amb els eixos? Dona una matriu de covariàncies que satisfaci aquesta condició i un vector de mitjanes qualsevol (no cal que generis dades).

Per una banda es parla d’el·lipcticitat perquè les variàncies de les tres mostres (X1, X2 i X3) són diferents, per tant, en comptes de projectar una esfera generen una forma el·liptica. Per altra banda, parlem d’angles amb els eixos perquè les covariàncies entre aquestes variables no és zero, sinó que tenen correlació.

Una forma el·líptica sense angles amb els eixos seria una multinormal amb variàncies diferents i covariàncies igual a 0, és a dir, a la diagonal de la matiru de covariàncies hi hauria tres valors diferents i a la resta zeros. A continuació es mostra un exemple:

\[ \sum\nolimits_{2} = \left( \begin{array}{c} 3 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 8 \end{array} \right) \]

1.d. Adjunta els dos fitxers en un nou data.frame, data de dimensió \(2 \cdot 10^{3} \times 3\) (una opció és utilitzar la funció rbind()). Afegeix una columna G (grup) al fitxer data, que prengui el valor 1 si la fila pertany al fitxer data1 i el valor 2 si és de data2. Els noms de les columnes del fitxer data han de ser X1, X2, X3 i G. Mostra els 5 primers i els 5 últims casos del fitxer data.

data <- rbind(data1, data2)
dim(data)
data$G <- c( rep(1, nsample), rep(2, nsample) )
## [1] 2000    3
data[1:5,]
data[c(nrow(data)-5):nrow(data), ]

2. Estimacions. Centrament i tipificació. Considerem les 3 variables numèriques del fitxer data2:

2.a. Calcula el vector amb les mitjanes mostrals de les variables amb la funció apply() i guarda’l en un objecte \(\mathrm{m}_{2}\). Calcula la matriu de covariàncies amb cov() i guarda-la en l’objecte \(\mathrm{S}_{2}\). Compara \(\mathrm{m}_{2}\) amb \(\mu_{2}\) i \(\mathrm{S}_{2}\) amb \(\sum_{2}\): comenta-ho.

m2 <- apply(X = data2, MARGIN = 2, FUN = mean)
S2 <- cov(data2)
all.equal(mu2, m2, check.attributes = F)
all.equal(S2, sigma2, check.attributes = F)
## [1] "Mean relative difference: 0.01054343"
## [1] "Mean relative difference: 0.1090626"

Observem que la diferència es gairebé zero, això es degut a que la mostra es genera aleatòriament i no encaixa a la perfecció amb la teòrica, però veiem que la diferència és mínima. A més, la diferència entre el vector de mitjanes teòriques i les mostrals la diferència és menor que entre la matriu de covariàncies teòriques i la mostral, això es deu a que la variància és una mesura de dispersó i valors molt atípics poden influenciar en el càlcul.

2.b. Calcula la matriu de correlacions mostrals de les variables de dues maneres. Comprova que et donen igual.

  • Amb la funció cor()
cor(data2)
##            X1         X2         X3
## X1  1.0000000  0.8087147 -0.6425905
## X2  0.8087147  1.0000000 -0.7219981
## X3 -0.6425905 -0.7219981  1.0000000
  • Usant \(R = D^{-1}SD^{-1}\) on \(D\) és la matiu diagonal amb les desviacions típiques mostrals de les variables.
d <- sqrt(diag(S2))
D <- diag(d, length(d), length(d))
R <- solve(D) %*% S2 %*% solve(D)
R
##            [,1]       [,2]       [,3]
## [1,]  1.0000000  0.8087147 -0.6425905
## [2,]  0.8087147  1.0000000 -0.7219981
## [3,] -0.6425905 -0.7219981  1.0000000
all.equal(cor(data2), R, check.attributes = F)
## [1] TRUE

2.c. Calcula la dimensió (nombre de files i nombre de columnes) i el rang de la matriu data2 i el mateix per a la matriu \(\mathrm{S}_{2}\).

dim(data2)
rk(data2)
dim(S2)
rk(S2)
## [1] 1000    3
## [1] 3
## [1] 3 3
## [1] 3

2.d. Dades centrades del fitxer data2. Digues X les dades centrades. S’obtenen restant a cada columna de data2 el valor mitjà d’aquesta columna. Nota: En aquest apartat, per operar amb data2, forceu as.matrix().

  • Què estem fent, geomètricament, quan centrem?

Es calcula la mitjana de cada variable de la matriu de les dades i es resta a cada punt de la columna. Per tant, les dades es desplacen i es distribueixen al voltant del punt d’origen de coordenades (el punt (0,0)).

  • Calculeu la matriu X amb la funció sweep() aplicada al fitxer de dades i al vector de mitjanes. Mostra la capçalera i la cua, calcula el rang de X.
X <- sweep(data.matrix(data2), 2, m2)
head(X)
tail(X)
rk(X)
##               X1           X2        X3
## [1,]  0.98771307  0.862465456 -3.933122
## [2,]  0.38062443 -0.007627684 -3.425084
## [3,] -0.02585107  0.347380760 -1.283411
## [4,]  1.42055887  1.138985860 -2.022635
## [5,]  3.31695823  1.643115908 -1.412015
## [6,] -0.16347744 -0.083140328  5.650118
##                  X1         X2        X3
##  [995,]  2.76541884  1.9198287 -1.112149
##  [996,]  2.32806520  1.5250842 -1.248813
##  [997,] -0.57933845 -0.4768471  1.903244
##  [998,] -1.56358699 -1.7085683  3.731422
##  [999,] -0.64361600 -0.4173417 -2.404586
## [1000,]  0.08092461 -0.7939769  5.247593
## [1] 3
  • Comprova que X es pot obtenir sense aplicar les funcions anteriors fent: \[ \mathrm{X} = \mathrm{data2} - 1_{n}m_{2}^t = \mathrm{data2} - \frac{1}{n}1_{n}1_{n}^{t} \mathrm{data2} \] on \(1_{n}\) és el vector columna de n uns. Nota: la solució d’aquest apartat és més general, atès que només utilitza funcions bàsiques, i es pot aplicar en qulsevol llenguatge de programació
n_1 <- as.matrix( rep(1, nsample) )
tm2 <- t(as.matrix(m2))

X2 <- as.matrix(data2) - n_1 %*% tm2 
all.equal(X, X2)
## [1] TRUE

Observem que el resultat és el mateix tant si es calcula amb la funció sweep() com si es fa manualment.

2.e. Comprova que la matriu de covariància es pot calcular amb les dades centrades: \(\frac{1}{n-1} \dot X^{t}X = S_{2}\).

S2_manual <- (1 / (nsample - 1)) * ( t(X) %*% X )
all.equal(S2, S2_manual)
## [1] TRUE

2.f. Com calcularies les dades tipificades, partint de les dades centrades X amb la fució sweep()? Fes-ho, pel fitxer data2. Digues Z la matriu de dades tipificades.

A partir de les dades centrades (guardades a la matriu X), dividint cada valor per la desviació típica. Anteriorment, he guardat a l’objecte d la desviació típica del fitxer dades.

Z <- sweep( data.matrix(X), 2, d, "/")

2.g. Comprova que la funció scale() dona tant les dades centrades com les dades tipificades, segons com fixis els arguments de la funció.

Z1 <- scale(data2, center = TRUE, scale = TRUE)

Per defecte a la funció scale() centra i tipifica les dades. En el cas que només vulguem centrar-les hem de posar dins de la funció center = T i scale = F. En el cas que només vulguem tipificar-les hem de posar dins de la funció center = F i scale = T. Si les volem centrar i tipificar, podem no indicar res ja que ho fa per defecte, però si ens volem assegurar de que ho faci hem de posar dins de la funció center = T i scale = T. En tots els casos anteriors s’ha d’entrar un objecte matricial on farà el càlcul especificat.

3. Visiualització

3.a. Pel fitxer final data, dibuixa un núvol de punts de les variables X2 (explicativa), X3 (resposta), amb colors segons el grup G: usar la funció scatterplot() de la llibreria car, posa noms als eixos i títols a la gràfica. Identifiques algun patró? Interpreta les línies que apareixen a la gràfica. És coherent la gràfica amb els paràmetres de les simulacions que has fet? Justifica les respostes.

scatterplot( X3 ~ X2 | as.factor(G), data = data, ellipse = T,
             xlab = "X2 (v. normal)",
             ylab = "X3 (v. normal)",
             main = "Núvol de punts",
             smooth = F,
             pch = c(2,2),
             legend = list(coords = "bottomleft", title= "Grup"),
             col=c("red", "blue") )


Les dues línies que travessen el núvol de punts representa el model lineal que li hem indicat a la funció scatterplot(), on X2 és l’explicativa i X3 la respsota. Mostra dos models lineals perquè li hem demanat que ens dibuixi els núvols de punts segons el grup en què pertanyen les dades, que pot ser del grup 1 o del 2.

Les dades de data1 no mostren una tendència perquè la covariància entre les dues variables és 0 (es recull en la matriu de covariàncies). A més, com que les variàncies són iguals s’hauria de veure una esfera, però a causa de l’escala del gràfic no es dibuixa bé. Els punts es concentren al voltant del punt de mitjanes $ (_{2},_3) = (1,3)$.

Les dades de data2 mostren una tendència decreixent, a causa que la covariància entre les dues variables (que es mostra en la matriu de covariàncies) és negativa. En aquest cas es dibuixa una el·lipse perquè les variàncies de les variables són diferents. Els punts es concentren al voltant del punt de mitjanes $ (_{2},_3) = (2,-10)$.

3.b. Utilitza la funció coplot() per obtenir els núvols de punts dels dos grups per separat i compara amb l’apartat anterior.

?coplot
coplot( X3 ~ X2 | as.factor(G), data = data)

Observem que ens dona els mateixos resultats que al gràfic anterior, però ho tenim en dues gràfiques per separat. En aquest cas les dades del grup 1 semblen seguir una estructura el·líptica en comptes de l’estructura esfèrica que s’espera, això és deu, de nou, a l’escala de la gràfica.

3.c. Matrius de gràfiques: dibuixa les parelles de gràfiques de les variables del fitxer data, fent que es vegin diferents colors dels punts segons el grup. Utilitza pairs() i scatterplotMatrix() de la llibreria car, etc. Indicacions:

  • Dins de pairs(), afegeix els arguments pch=21, `bg=rainbow(2)[unclass(data$G)]. Què fa aquesta opció?

Aquesta opció crea un vector de 2 colors que els adjudica a cada grup diferent que hi ha a la variable data$G.

*De vegades, convé especificar el paquet davant la funció `car::scatterplotMatrix().

pairs(data[1:3], pch = 21, bg = rainbow(2)[data$G], main = "Núvol de punts segons variables" )

En aquest gràfic s’observa més detalladament la dispersió de les dades segons les variables. En vermell es dibuixen les dades pertanyents al grup 1 i en blau les del grup 2.

En el gràfic del mig a baix de tot és el mateix que s’ha dibuixat anteriorment, on X2 és l’explicativa i X3 és la resposta. A l’esquerra d’aquest es dibuixa el núvol de punts quan X1 és l’explicativa i X3 la resposta. Veiem que és molt similar al de la seva dreta, però hi ha un canvi en la tendència dels punts blaus i la posició en el gràfic, això es deu a la correlació entre les variables i a la mitjana de X1, respectivament.

Els núvols de punts vermells hauríem d’observar una forma esfèrica, però per culpa de l’escala del gràfic veiem més aviat una el·lipse que no una esfera (recordem que les variàncies de les variables són iguals). En el cas on es veu més el dibuix de l’esfera és al gràfic on X1 és l’explicativa i X2 la resposta.

scatterplotMatrix( ~ X1 + X2 + X3 | as.factor(G), data = data, by.group = T, col = rainbow(2))

Aquest gràfic és molt similar al executat a l’apartat anterior perquè ens mostra el núvol de punts de les dades, però també mostra el dibuix del model linial entre les dues variables i la densitat de les dades a la diagonal.

Al primer gràfic de la diagonal s’observa com les dades referents al grup 2 en la variable X1 la mitjana de les dades es troba aproximadament al voltant de -3, valor de la mitjana teòrica. En general, les densitats dibuixades segueixen l’estructura de la campana de Gauss, però a causa dels diferents valors dels paràmetres les corbes varien.

En les variables X1 i X2 les campanes es creuen, ja que els paràmetres tenen valors més similars. En canvi, a la variable X3 les campanes no es creuen en cap moment, ja que els valors dels paràmetres són força diferents. Si dibuixessim els valors de les dades centrades sí que se superposarien.

3.d. Al fitxer data1, es vol estimar la densitat bivariant de X2 i X3. Usa la funció bivden() del fitxer bivden.R. Indicacions per a aquest fixer i el de data2 i data:

  • Carrega la funció o funcions que hi hagi dins de l’escript .R a l’espai de treball fent source(nomdelfitxer.extensio).
source("bivden.R")
  • Seguidament, aplica la funció bivden() a les dades d’interès i guarda l’output en un objecte que anomenaràs dens.
dens1 <- bivden(x = data1$X2, y = data1$X3, ngridx = 10, ngridy = 10)
  • Representa l’output dens$seqx, dens$seqy i dens$den d’aquesta densitat amb filled.contour() de la llibreria grDevices.
filled.contour(x = dens1$seqx, y = dens1$seqy, z = dens1$den, color.palette = terrain.colors,
               main = "Densitat data1")

El gràfic mostra una escala percentual de la concentració dels punts. El color més clar és on hi ha una major concetració dels punts de les dades (al voltant dels valors de les mitjanes teòriques), i a mesura que ens allunyem d’aquest punt la concentració disiminueix. La densitat dibuixada és la de dues variables amb la mateixa variància, per això veiem que l’estructura que segueix és circular.

dens2 <- bivden(x = data2$X2, y = data2$X3, ngridx = 10, ngridy = 10)
filled.contour(x = dens2$seqx, y = dens2$seqy, z = dens2$den, color.palette = terrain.colors,
               main = "Densitat data2")

La major densitat es troba al voltant del punt de les mitjanes teòriques ( \((\mu_2, \mu_3) = (2, -10)\)). En aquest cas la densitat segueix una estructura no circular, degut a que les variàncies tenen valors diferents. A més hi ha inclinació perquè les covariàncies són diferents de zero.

dens3 <- bivden(x = data$X2, y = data$X3, ngridx = 10, ngridy = 10)
filled.contour(x = dens3$seqx, y = dens3$seqy, z = dens3$den, color.palette = terrain.colors,
               main = "Densitat data")

En aquest cas, la gràfica de la densitat veiem que hi ha dos punts de concentració de dades, un respon a les dades del grup1 i les altres a les del grup2. En el cas que centréssim les dades, tindriem una concentració de punts al punt d’origen de coordenades.

  • Explora també la funció persp(), amb els seus paràmetres.
persp(x = dens1$seqx, y = dens1$seqy, z = dens1$den,
      theta = 120,
      col = "aquamarine3",
      xlab = "X", ylab = "Y", zlab="dens",
      main = "Densitat 3D de d'una multinormal de dos dimensions (data1)")

En el gràfic es pot observar com la curtosi que dibuixa la densitat de les variables és bastant simètrica i apuntalada, per tant es podria mirar si alguna de les dues variables pateix leptocurtosis.

persp(x = dens2$seqx, y = dens2$seqy, z = dens2$den,
      theta = 120,
      col = "aquamarine3",
      xlab = "X", ylab = "Y", zlab="dens",
      main = "Densitat 3D de d'una multinormal de dos dimensions (data2)")

En aquest cas la densitat en 3D de la multinomial dibuixa una curtosis no tan simètrica com la que hem vist anterioment, això es deu a que les variables no tenen la mateixa variància.

persp(x = dens3$seqx, y = dens3$seqy, z = dens3$den,
      theta = 120,
      col = "aquamarine3",
      xlab = "X", ylab = "Y", zlab="dens",
      main = "Densitat 3D de d'una multinormal de dos dimensions (data2)")

En aquest cas es dibuixen dues curtosis diferents que fa referència als dos grups que formen les dades.

3.e. Gràfiques interactives. Per a una representació en 3D, de les variables numèriques del fitxer data1, utilitza les funcions scatter3d() de la llibreria car i la funció plot3D() de la llibreria rgl, amb diveres opcions. Fes que X3 sigui la variable resposta, X1 i X2 les explicatives. Idem per al fitxer data2. Idem per al fitxer data. Interpreta’ls.

car::scatter3d(X3 ~ X1 + X2, data = data1, sphere.size=.4)
plot3d(data1$X1, data1$X2, data1$X3, box=F) 
m1 <- apply(data1, 2, mean)
S1 <- cov(data1)
plot3d( ellipse3d( S1, centre = m1), col="green", alpha = 0.2, add = TRUE)
rglwidget()

Les dades tenen variàncies iguals i covariància zero, per tant, el dibuix de les dades (que es pot veure en el gràfic) segueix la forma esfèrica. Podem dir que els valors en diferents moments no estan correlacionats.

car::scatter3d(X3 ~ X1 + X2, data = data2, sphere.size=.4)
plot3d(data2$X1, data2$X2, data2$X3, box=F) 
m2 <- apply(data2, 2, mean)
S2 <- cov(data2)
plot3d( ellipse3d( S2, centre = m2), col="green", alpha = 0.2, add = TRUE)
rglwidget()

En aquest cas veiem que la forma que dibuixen les dades és el·líptica i té tendència perquè les variàncies són diferents i les covariàncies no són iguals a zero.

car::scatter3d(X3 ~ X1 + X2, data = data, sphere.size=.4)
plot3d(data$X1, data$X2, data$X3, box=F) 
md <- apply(data, 2, mean)
Sd <- cov(data)
plot3d( ellipse3d( Sd, centre = md), col="green", alpha = 0.2, add = TRUE)
rglwidget()

En aquest gràfic podem observar com la concentració de punts tenen dos tendències diferents, una fa referència a les dades del grup1 i l’altre a les del grup2. En aquest cas la forma dibuixada és una el·lipse degut a la variància diferent i convariàncies diferents de zero de les variables de data2.

options(rgl.useNULL = TRUE) 
car::scatter3d(X3 ~ X1 + X2, data = data1, sphere.size=.4) 
rglwidget()

Observem que el pla que es dibuixa no té tendència o és gairebé nul·la degut a la coviarància igual a zero entre les variables.

options(rgl.useNULL = TRUE) 
car::scatter3d(X3 ~ X1 + X2, data = data2, sphere.size=.4) 
rglwidget()

En aquest cas el pla te tendència perquè les covariàncies són diferents de zero.

options(rgl.useNULL = TRUE) 
car::scatter3d(X3 ~ X1 + X2, data = data, sphere.size=.4) 
rglwidget()

De nou, s’observa la tendència del pla que es dibuixa degut a les característiques de data2. A més, com que es dibuixen els dos grups al mateix gràfic de nou hi ha dues concentracions de punts diferents en les dades.

options(rgl.useNULL = TRUE) 
car::scatter3d(X3 ~ X1 + X2 | as.factor(G), data = data, sphere.size=.4) 
rglwidget()

En cas de dibuixar el gràfic segons el grup de dades que pertanyen veiem que els dos plans es troben separats, això es deu als valors de les mitjanes i variàncies de les dades que fan que no es trobin els punts superposats.