---
title: "T4. Taller encuestas. Análisis factorial AF^[*doc:`r knitr::current_input()`*]"
author: Álvaro Hernández Vicente, Elvira Ferre Jaén, Antonio José Perán Orcajada,
  Ana Belén Marín Valverde, Antonio Maurandi López^[Servicio de Apoyo Estadístico;
  <alvarohv@um.es>, <elvira@um.es>, <antoniojose.peran@um.es>, <anabelen.marin4@um.es>,
  <amaurandi@um.es>]
date: "17 de noviembre de 2016"
output:
  pdf_document:
    fig_caption: yes
    includes:
      in_header: ./plantillas/header-sae1.tex
    number_sections: yes
    toc: yes
    toc_depth: 3
  html_document:
    toc: yes
    toc_depth: '3'
header-includes: \usepackage{booktabs}
subtitle: VIII Jornadas R Albacete 2016^[<http://r-es.org/8jornadasR/>]
bibliography: ./referencias/ref.bib
---

```{r setup, include=FALSE}
library( ggplot2 )
library( knitr )

# install.packages("corpcor")
# install.packages("GPArotation")
# install.packages("psych")

library( corpcor )
library( GPArotation )
library( psych )

knitr::opts_chunk$set( echo = TRUE, warning = FALSE )



```


# Lectura de datos

Para llevar a cabo todo el procedimiendo que se realiza en este documento seguimos el capítulo 17: Exploratory factor analysis de *Discovering Statistics Using R* [@field12].

```{r datos}
load( "saeraq.RData" )
```

# Comprobaciones previas

Es aconsejable realizar las siguientes comprobaciones antes de proceder con el análisis factorial.

* Observar la matriz de correlaciones: en busca de variables poco o demasiado correlacionadas.
* Test de Bartlett
* Coeficiente Kaiser-Meyer-Olkin (KMO)
* Determinante de la matriz de correlaciones.

## Matriz de correlaciones

Podemos trabajar indistintamente con el conjunto de datos o con la matriz de correlaciones, pero en cualquier caso, siempre con variables numéricas, no con factores. Se aconseja trabajar con la matriz de correlaciones, ya que así evitamos recalcular dicha matriz cada vez que llamemos a alguna función del análisis.

```{r preparardatos, echo=FALSE}

df <- as.data.frame(lapply(df[, 7:ncol(df)], as.numeric))
names(df) <- dicc$item[7:nrow(dicc)]
```

En nuestro caso, calculamos la matriz de correlaciones:

\footnotesize

```{r corrmatrix}
corrMatrix <- cor(df)
round(corrMatrix, 3)
```

\normalsize

Solo con esta matriz ya podemos extraer información importante de los datos que nos servirá para realizar el análisis correctamente. El propósito principal del `Análisis Factorial` y del `Análisis de Componentes Principales` es el de encontrar constructos formados por variables que correlacionan bien entre ellas, pero no perfectamente. 

Por tanto, antes de llevar a cabo el análisis, se recomienta examinar la matriz de correlaciones en busca de variables que no correlacionen bien con ninguna otra, es decir, con coeficientes de correlación todos menor que 0.3, y variables que correlacionen demasiado bien con otras, es decir, variables que tengan algún coeficiente de correlación mayor que 0.9.
Las primeras deberíamos eliminarlas del análisis, y las segundas podremos mantenerlas, pero teniendo en cuentra que quizá causen problemas de multicolinealidad.

```{r}
cat("Nº de variables que no correlacionan bien: ",
            sum( sapply(as.data.frame(corrMatrix), function(x) all(x < 0.3)) ) )

cat("Nº de variables que podrían causar multicolinealidad: ",
            sum( sapply(as.data.frame(corrMatrix), function(x) any(x >= 0.9 & x != 1)) ) )

```

En nuestro caso, no tenemos problemas de este tipo.



## Test de Bartlett

Buscamos que el test de Bartlett salga significativo, es decir, que nuestra matriz no sea similar a una matriz identidad.

```{r bartlett}
library(psych)
cortest.bartlett(corrMatrix, nrow(df))
```



## Coeficiente Kaiser-Meyer-Olkin (KMO)

**Ver:** http://tolstoy.newcastle.edu.au/R/e2/help/07/08/22816.html

```{r kmo, eval=TRUE, echo=FALSE}
# KMO Kaiser-Meyer-Olkin Measure of Sampling Adequacy 
kmo <- function( data ){

  library(MASS)
  X <- cor(as.matrix(data))
  iX <- ginv(X)
  S2 <- diag(diag((iX^-1)))

  AIS <- S2%*%iX%*%S2                      # anti-image covariance matrix
  IS <- X+AIS-2*S2                         # image covariance matrix
  Dai <- sqrt(diag(diag(AIS)))
  IR <- ginv(Dai)%*%IS%*%ginv(Dai)         # image correlation matrix
  AIR <- ginv(Dai)%*%AIS%*%ginv(Dai)       # anti-image correlation matrix

  a <- apply((AIR - diag(diag(AIR)))^2, 2, sum)
  AA <- sum(a)
  b <- apply((X - diag(nrow(X)))^2, 2, sum)
  BB <- sum(b)

  MSA <- b/(b+a)                        # indiv. measures of sampling adequacy

  AIR <- AIR-diag(nrow(AIR))+diag(MSA)  # Examine the anti-image of the
                                        # correlation matrix. That is the
                                        # negative of the partial correlations,
                                        # partialling out all other variables.

  kmo <- BB/(AA+BB)                     # overall KMO statistic

  # Reporting the conclusion
    if (kmo >= 0.00 && kmo < 0.50){
      test <- 'The KMO test yields a degree of common variance unacceptable for FA.'

} else if (kmo >= 0.50 && kmo < 0.60){
      test <- 'The KMO test yields a degree of common variance miserable.'
} else if (kmo >= 0.60 && kmo < 0.70){
      test <- 'The KMO test yields a degree of common variance mediocre.'
} else if (kmo >= 0.70 && kmo < 0.80){
      test <- 'The KMO test yields a degree of common variance middling.'
} else if (kmo >= 0.80 && kmo < 0.90){
      test <- 'The KMO test yields a degree of common variance meritorious.'
} else {

      test <- 'The KMO test yields a degree of common variance marvelous.'     }

    ans <- list( overall = kmo,

                  report = test,
                  individual = MSA,
                  AIS = AIS,
                  AIR = AIR )


    return(ans)

} # end of kmo() 

```

* Es una medida de la comparación de los coeficientes de correlación observados con los coeficientes de correlación parcial que asume valores entre 0 y 1.

* No deben aceptarse valores menores a 0.5 según Kaiser(1974); si el valor está entre 0.5 y 0.7 se considera mediocre, si está entre 0.7 y 0.9 se considera bueno, y si es mayor que 0.9 se considera muy bueno, según Hutcheson & Sofroniou (1999).

* El KMO puede calcularse de forma global o para cada una de las variables. Deberemos considerar eliminar las variables con KMO individual inferior a 0.5, y una vez sustraídas, recalcularlo.

```{r, kmocalculo}

kmo(df)[1:3]

```

En nuestro caso, el KMO es muy bueno, y todas las variables tienen un KMO individual considerablemente alto.


## Determinante de la matriz de correlaciones

Para poder realizar el análisis factorial, la matriz de correlaciones ha de ser definida positiva, es decir, al menos su determinante debe ser mayor que 0. Se recomienta que sea mayor que 0.00001 (Andy Field?)

```{r determinante}
cat("Valor del Determinante: ", det(corrMatrix))
```

# Análisis de Componentes Principales

## Determinación del número de componentes principales

Para determinar el número de componentes principales a mantener, realizamos un análisis *a priori* en el que el número de componentes sea igual al número de variables, y elegimos el número de componentes a extraer de acuerdo a alguno de los siguientes criterios.

* Criterio de Kaiser: en azul, escogeríamos tantas componentes como valores propios superen el valor 1.
* Criterio de Jolliffe: en rojo,  las que superen el 0.7.

```{r screeplot, echo=FALSE}
pc1 <- principal(corrMatrix, nfactors = ncol(df), rotate = "none")

ggplot(mapping = aes(x = 1:23, y=pc1$values) ) +
  geom_point() +
  geom_hline( aes( yintercept = 1, colour = "Kaiser" ), show.legend = TRUE ) +
  geom_hline( aes( yintercept = 0.7, colour = "Jolliffe" ), show.legend = TRUE ) +
  labs(title = "Scree Plot", x = "Componente Principal",
       y = "Valores Propios", colour = "Criterio" )
```

En nuestro caso, retendremos 4 componentes principales de acuerdo al Criterio de Kaiser.

\newpage

## Análisis

```{r analisis}
pc1 <- principal(corrMatrix, nfactors = 4, rotate = "none")
pc1
```

En la primera matriz que muestra el resultado anterior tenemos, en las distintas columnas:

* Las cargas de cada una de las componentes en cada una de las variables (PC1, PC2, PC3, PC4).
* Las comunalidades (h2), que se interpreta como el porcentaje de varianza de cada variable que explican las componentes seleccionadas.
* La varianza única o uniquenesses (u2), que se interpreta como lo que queda sin explicar.
* Los índices de complejidad de Hoffman (com).

Hay una serie de pautas que pueden ayudarnos a la hora de comprobar que nuestra elección del *Criterio de Kaiser* para la extracción de factores ha sido buena. Por ejemplo, si el tamaño muestral es es menor de 30 y las comunalidades son todas mayores a 0.7; o cuando el tamaño muestral es superior a 250 y la media de las comunalidades es superior a 0.6.

En nuestro caso, dicha media es `r mean(pc1$communality)`, próxima a 0.6, luego podemos considaderarla buena debido al gran tamaño de nuestra muestra: `r nrow(df)` individuos. 

## Medidas de ajuste del modelo

Las medidas de ajuste del modelo se basan en estimar de distintas maneras el tamaño de los residuos que resultan al comparar la matriz de correlaciones original con la matriz de correlaciones reproducida por el modelo.

Lo ideal sería que la matriz de correlaciones reproducida por el modelo fuese lo más parecida posible a la calculada con los datos.

Una buena medida de ajuste viene dada por la función `principal`,

```{r, echo=FALSE}
cat(" Fit based upon off diagonal values = ", pc1$fit.off)
```

y resulta de dividir la suma de los cuadrados de los residuos por la suma de los cuadrados de los coeficientes de la matriz de correlación original.

Otras formas de medir el ajuste interesantes podemos realizarlas con la siguiente función provista en *Discovering Statistics Using R*, Andy Field (2012).

```{r medidas-de-ajuste, eval=TRUE}
# library(psych)

residual.stats <- function(corrMatrix, model){
  matrix <- factor.residuals(corrMatrix, model$loadings)
  residuals <- as.matrix(matrix[upper.tri(matrix)])
  large.resid <- abs(residuals) > 0.05
  numberLargeResids <- sum(large.resid)
  propLargerResid <- numberLargeResids/nrow(residuals)
  rmsr <- sqrt(mean(residuals^2))
  
  cat("Root means squared residual = ", rmsr, "\n")
  cat("Number of absolute residuals > 0.05 = ", numberLargeResids, "\n")
  cat("Proportion of absolute residuals > 0.05 = ", propLargerResid, "\n")
  hist(residuals)
}

```

```{r calculomedidasajuste}
library(psych)
residual.stats(corrMatrix, pc1)

```

## Rotación

El procedimiento de rotar los factores permite una mejor interpretación de la estructura subyacente en los datos. El método consiste en aumentar las cargas de las variables sobre alguna componente y minimizarlas sobre otros sin alterar las comunalidades, aunque sí la varianza explicada por cada componente.

```{r rotacion}
pc1 <- principal(corrMatrix, nfactors = 4, rotate = "varimax")
pc1
```

Haciendo uso de la función `print.psych` del paquete `psych` podemos mostrarlos de forma ordenada y seleccionando únicamente las cargas que nos interesen, por ejemplo, las mayores que 0.3.

```{r rotacionprint}
library(psych)
print.psych(pc1, cut = 0.3, sort = TRUE)

```

## Interpretación

Basándonos en la salida anterior, podemos ver que los ítems 6, 18, 13, 7, 14, 10, y 15, que se corresponden con los enunciados

```{r constructo1, echo=FALSE, eval=TRUE}
dicc$spanish[c(6, 18, 13, 7, 14, 10, 15)+6]
```

los cuales tienen en común que todos representan alguna forma de expresar cierto miedo hacia los computadores, por ello, podemos pensar que esta componente representa exactamente eso, el `Miedo a los computadores`.

La siguiente componente carga fuertemente sobre los ítems 20, 21, 3, 12, 4, 16, 1 y 5, que se corresponden con los enunciados

```{r constructo2, echo=FALSE, eval=TRUE}
dicc$spanish[c(20, 21, 3, 12, 4, 16, 1, 5)+6]
```

los cuales representan el `Miedo a la estadística`.

La tercera componente carga sobre los ítems 8, 17 y 11, que se corresponden con los enunciados

```{r constructo3, echo=FALSE, eval=TRUE}
dicc$spanish[c(8, 17, 11)+6]
```

los cuales representan el `Miedo a las matemáticas`.

Y por último, la cuarta componente está formada por los ítems 9, 22, 23, 2, y 19, que se corresponden con los enunciados

```{r constructo4, echo=FALSE, eval=TRUE}
dicc$spanish[c(9, 22, 23, 2, 19)+6]
```

los cuales representan la `Evaluación por pares`, pues todos contienen cierta componente de opinión que los demás puedan tener sobre uno mismo .

## Cargas de individuos

Para obtener las cargas de cada individuo sobre cada componente, es decir, las "coordenadas" de cada individuo, solo es necesario ejecutar

```{r puntuacionesindividuos}
pc1 <- principal(df, nfactors = 4, rotate = "varimax", scores = TRUE)
pc1$scores[1:10, ]
```

## Resumen

En la siguiente tabla se muestra un resumen de la interpretación que hemos llevado a cabo.

```{r resumen, echo=FALSE, eval=TRUE}
componentes <- c("Miedo a los computadores", "Miedo a la Estadística", "Evaluación por pares", "Miedo a las Matemáticas")
items <- dicc$english[7:nrow(dicc)]
ordenitems <- c(6,18,13,7,14,10,15,20,21,3,12,4,16,1,5,8,17,11,9,22,23,2,19)

kk <- as.data.frame(round(pc1$loadings[1:23, 1:4], 3))

colnames(kk) <- componentes
rownames(kk) <- items

kk <- kk[ordenitems, ]

kable( kk, caption = "Tabla Resumen de las Componentes Principales" )
```


# Análisis de Fiabilidad

El coeficiente de fiabilidad más ampliamente utilizado es el $\alpha$ de Cronbach. Suele hablarse, en ocasiones y de forma errónea, sobre un $\alpha$ de Cronbach general para una encuesta. si tenemos claro que en cierta encuesta subyacen distintos constructos, como es nuestro caso, debe calcularse un $\alpha$ para cada uno de ellos, teniendo en cuenta, además, que todos lo ítems apunten hacia el mismo sentido.

Los constructos de nuestra encuesta son

```{r constructos}
miedoComputadores <- c(6, 7, 10, 13, 14, 15, 18) 
miedoEstadistica  <- c(1, 3, 4, 5, 12, 16, 20, 21)
miedoMatematicas  <- c(8, 11, 17)
evaluacionPares   <- c(2, 9, 19, 22, 23)
```

## $\alpha$ de Cronbach para el `Miedo a los computadores`

```{r alpha1}
library( psych )
psych::alpha( df[ , miedoComputadores] )
```

La salida anterior nos ofrece mucha información sobre la fiabilidad de nuestro constructo midiendo lo que pretende. Por un lado, `raw_alpha` y `std.alpha` son los $\alpha$ de Cronbach basado en la covarianza y la correlación respectivamente. El coeficiente `G6(smc)` es el Guttman's Lambda 6 reliability?? y `average_r` es la correlación media inter-ítem.

La tabla `Reliability if an item is dropped` nos indica cómo cambiarían cada uno de los coeficientes anteriores al eliminar cierto ítem. Si alguno de estos coeficientes mejorase al eliminar algún ítem, deberíamos plantearnos excluírlo de este constructo, pues no lo estaría reflejándolo de forma precisa.

En la tabla `Item Statistics` tenemos, entre otras cosas, la media y desviación típica de las puntuaciones de cada ítem en las columnas `mean` y `sd` respectivamente. Por otro lado, `raw.r` y `std.r` nos dan la correlación de cada ítem con todos los demás ítems de la escala (incluyendo el ítem en cuestión), cuando los ítems están o no estandarizados respectivamente.  

La última tabla es simplemente una tabla de frecuencias, que indica el porcentaje de individuos que han respondido cada respuesta en cada uno de los ítems que conforman el constructo.

## Alfa de Cronbach para el `Miedo a la estadística`

Observar que aquí cambiamos la dirección de ítem 3, cuyo enunciado es

```{r item3}
dicc$spanish[9]
```


```{r alpha2}
# library( psych )
psych::alpha( df[ , miedoEstadistica], keys = c(1, -1, 1, 1, 1, 1, 1, 1) )
```

## $\alpha$ de Cronbach para el `Miedo a las matemáticas`

```{r alpha3}
# library( psych )
psych::alpha( df[ , miedoMatematicas] )
```

## $\alpha$ de Cronbach para el `Evaluación por pares`

```{r alpha4}
# library( psych )
psych::alpha( df[ , evaluacionPares] )
```

\newpage


# Referencias y bibliografía





