Воркшоп по факторному анализу (EFA, CFA)
Баранов Дмитрий
Используемые пакеты:
library(lavaan)
#install.packages("lavaan")
library(formattable)
#install.packages("formattable")
library(dplyr)
#install.packages("dplyr")
library(psych)
#install.packages("psych")
library(polycor)
#install.packages("polycor")
library(nFactors)
#install.packages("nFactors")
library(ggpubr)
#install.packages("ggpubr")
library(semTools)
#install.packages("semTools")
library(lavaanPlot)
#install.packages("lavaanPlot")
library(EFA.dimensions)
#install.packages("EFA.dimensions")Перед началом:
Факторный анализ - достаточно сложная техника, по каждому шагу проведения которой ведется большое количество споров. Поэтому к заявлениям во время воркшопа нельзя относиться как к единственно верным. Чтобы больше погрузиться в тему, есть большое количество статей, блогов и хэндбуков, которые позволяют намного более полно понять эту сложную тему.
Вот некоторые книги для изучения:
Brown, T. A. (2015). Confirmatory factor analysis for applied research. Guilford publications.
Kline, R. B. (2023). Principles and practice of structural equation modeling. Guilford publications.
Классическая теория тестирования (Classical test theory, CTT)
Bandalos, D. L. (2018). Measurement theory and applications for the social sciences. Guilford Publications.
Тесты играют большую роль в нашей жизни… Для повышения надежности (reliability) оценки какого-либо показателя (например, эмоциональный интеллект, самооценка, депрессивность) задается несколько вопросов или дается несколько заданий похожего содержания, чтобы минимизировать потенциальную ошибку отдельного вопроса/задания (непонятно задан, ответ был “очень личным” и т.д.).
КТТ предполагает:
X = T + E
Где Х - полученное значение (observed score), T - реальное значение (true score), E - ошибка (error)
Например, Вася отвечает на вопрос “Насколько вы удовлетворены своей жизнью?” по шкале от 1 до 5 (где 1 - совсем не удовлетворен, 5 - полностью удовлетворен). Васин ответ: 4
Этот ответ складывается из того, насколько Вася действительно в целос доволен жизнью (T) и некоторых ошибок (E): как Вася понял конкретно этот вопрос, выспался ли он сегодня, какое у него настроение и т.д.
Если несколько похожих вопросов дают коррелирующие ответы, то их надёжность высокая, и из ответов на эти вопросы можно попробовать собрать фактор или метрику. Методы оценки надежности - отдельная тема, которую мы коснёмся частично и только в контексте факторного анализа (ФА)
Многие идеи факторного анализа строятся на идее КТТ, к чему мы вернемся чуть позже.
Факторный анализ
Факторный анализ - стат. метод, который позволяет из набора отдельных переменных создавать факторы - метрики, которые трудно измерить напрямую одним вопросом или тестовым заданием (в случае социальных наук, например, удовлетворенность жизнью, интеллект, самооценка)
Пример из реальной практики: общая самооценка (по Розенбергу) - это фактор, который “собирается” из ответов на следующие вопросы:
https://www.apa.org/obesity-guideline/rosenberg-self-esteem.pdf
Please record the appropriate answer for each item, depending on whether you Strongly agree, agree, disagree, or strongly disagree with it.
1 = Strongly agree, 2 = Agree, 3 = Disagree, 4 = Strongly disagree
On the whole, I am satisfied with myself.
At times I think I am no good at all.
I feel that I have a number of good qualities.
I am able to do things as well as most other people.
I feel I do not have much to be proud of.
I certainly feel useless at times.
I feel that I’m a person of worth.
I wish I could have more respect for myself.
All in all, I am inclined to think that I am a failure.
I take a positive attitude toward myself.
Таким образом, можно создавать тесты (опросы), которые с помощью некоторых переменных (items, это 10 вопросов от Розенберга) позволяют создавать факторы (или латентные переменные, это общая самооценка в нашем случае)
Существует два основных подхода к факторному анализу:
Эксплораторный ФА (EFA)
Не предполагает изначально заданной факторной структуры*
Конфирматорный ФА (CFA)
Проверяет насколько хорошо уже предполагаемая нами факторная структура ложится на данные
Основные шаги проведения факторного анализа:
Изучить теорию по теме, познакомиться с уже созданными исследователями шкалами, проведенными ФА по теме
Подготовить опрос/тест и собрать данные*
Выбрать метод анализа (EFA или CFA отдельно? оба на одних данных? и т.д.) (сейчас по этому поводу будет долгий разговор)
Провести анализ, доработать получившиеся факторы, оценить качество модели
Провести дальнейший анализ, ответить на поставленные исследовательские вопросы* (например, посмотреть, как получившиеся значения факторов распределены относительно возраста/пола)
Важно: по сути своей факторный анализ - это никакой не интеллект, который решит все ваши проблемы и всё за вас придумает, это набор чудес линейной алгебры, не более
В факторном анализе компьютер создаёт идеальную математическую модель и накладывает на эту идеальную модель наши данные
Построение теста/опроса
Если у вас нет данных для анализа и вы хотите их собрать:
- Поймите, что хотите изучить, поставьте исслед. вопросы, проработайте дизайн исследования
Допустим, вы хотите изучить, как общая самооценка взаимосвязана с коммуникабельностью у подростков. Составьте для каждого фактора пул вопросов (не меньше трёх на фактор, лучше четыре), вопросы можно найти в уже существующих и проверенных шкалах из других исследований (что лучше) или создать самостоятельно (велик риск, что получится не очень хорошо)
- Составьте опрос и проведите его (собрать лучше от 300 наблюдений*)
Возможно, опрос будет выглядеть так:
- Социально демографическая информация (пол, возраст и т.д.)
- 10 вопросов из Розенберга
- 4 придуманных вами вопроса на коммуникабельность (или найденных в других исследованиях, например, тут:
Можно: а) перевернуть некоторые вопросы, относящиеся к одному фактору (вместо “Как часто вы чувствуете себя счастливым?” спросить “Как часто вы чувствуете себя совершенно несчастным?”), б) определить порядок вопросов (перемешивать ли между собой вопросы, которые относятся к разным факторам?)
Важно: если у кого-то тот или иной инструментарий (опрос) показал хорошие результаты после факторного анализа, это не значит, что у вас получится так же хорошо, ввиду сложностей перевода, культурных особенностей, выборки или других проблем.
Данные для анализа
Количество наблюдений
Точного правила о количестве наблюдений для факторного анализа нет.
Лучше иметь более 300 полных наблюдений. Если какие-то наблюдения неполные, можно:
исключить неполные наблюдения (но лучше быть уверенным, что пропуски случайные, а не систематические)
импутировать пропущенные наблюдения (напр., есть пакеты mice и missForrest, но выбирать лучший метод импутации и оценивать её качество - сложная задача)
Негативно заданные вопросы
В ординальных данных будет удобно “перевернуть” негативно заданные вопросы (в нашем примере про счастье, если на вопрос “Как часто вы чувствуете себя совершенно несчастным человеком?” ответили “скорее согласен”, значение лучше перевернуть в “скорее не согласен”)
Проверьте распределение данных
Это поможет провалидировать данные, понять, что с ними уж точно можно дальше работать
Данные
Для анализа будем использовать сокращенный классический датасет Холцингера и Свайнфорда.
https://search.r-project.org/CRAN/refmans/psychTools/html/holzinger.swineford.html
Данные содержат информацию о школьниках 7-8 классов, которые проходили различные тесты. Каждая переменная - результат отдельного теста на те или иные интеллектуальные способности, также есть данные о поле, возрасте и учебном классе участников исследования.
x1 (тест 1) — Тест зрительного восприятия Spearman VPT (оценка анализа геометрических паттернов)
x2 (тест 2) — Упрощённый тест пространственных отношений Brigham (сборка кубов по образцу)
x3 (тест 4) — Тест «Лозенги» Thorndike (идентификация перевёрнутых фигур среди похожих)
x4 (тест 6) — Тест понимания текста (ответы на вопросы по прочитанному абзацу)
x5 (тест 7) — Тест завершения предложений (подбор слов по контексту)
x6 (тест 9) — Тест значения слов (поиск синонимов/антонимов)
x7 (тест 10) — Тест скоростного сложения (быстрые арифметические вычисления)
x8 (тест 12) — Тест счёта точек (подсчёт точек внутри фигур на скорость)
x9 (тест 13) — Тест дискриминации линий (разделение прямых и кривых линий за время)
library(lavaan)
data = HolzingerSwineford1939
library(formattable)
formattable(head(data, n = 10))| id | sex | ageyr | agemo | school | grade | x1 | x2 | x3 | x4 | x5 | x6 | x7 | x8 | x9 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 1 | 13 | 1 | Pasteur | 7 | 3.333333 | 7.75 | 0.375 | 2.333333 | 5.75 | 1.2857143 | 3.391304 | 5.75 | 6.361111 |
| 2 | 2 | 13 | 7 | Pasteur | 7 | 5.333333 | 5.25 | 2.125 | 1.666667 | 3.00 | 1.2857143 | 3.782609 | 6.25 | 7.916667 |
| 3 | 2 | 13 | 1 | Pasteur | 7 | 4.500000 | 5.25 | 1.875 | 1.000000 | 1.75 | 0.4285714 | 3.260870 | 3.90 | 4.416667 |
| 4 | 1 | 13 | 2 | Pasteur | 7 | 5.333333 | 7.75 | 3.000 | 2.666667 | 4.50 | 2.4285714 | 3.000000 | 5.30 | 4.861111 |
| 5 | 2 | 12 | 2 | Pasteur | 7 | 4.833333 | 4.75 | 0.875 | 2.666667 | 4.00 | 2.5714286 | 3.695652 | 6.30 | 5.916667 |
| 6 | 2 | 14 | 1 | Pasteur | 7 | 5.333333 | 5.00 | 2.250 | 1.000000 | 3.00 | 0.8571429 | 4.347826 | 6.65 | 7.500000 |
| 7 | 1 | 12 | 1 | Pasteur | 7 | 2.833333 | 6.00 | 1.000 | 3.333333 | 6.00 | 2.8571429 | 4.695652 | 6.20 | 4.861111 |
| 8 | 2 | 12 | 2 | Pasteur | 7 | 5.666667 | 6.25 | 1.875 | 3.666667 | 4.25 | 1.2857143 | 3.391304 | 5.15 | 3.666667 |
| 9 | 2 | 13 | 0 | Pasteur | 7 | 4.500000 | 5.75 | 1.500 | 2.666667 | 5.75 | 2.7142857 | 4.521739 | 4.65 | 7.361111 |
| 11 | 2 | 12 | 5 | Pasteur | 7 | 3.500000 | 5.25 | 0.750 | 2.666667 | 5.00 | 2.5714286 | 4.130435 | 4.55 | 4.361111 |
Посмотрим на структуру данных, чтобы убедиться, что типы переменных стоят правильно
str(data)## 'data.frame': 301 obs. of 15 variables:
## $ id : int 1 2 3 4 5 6 7 8 9 11 ...
## $ sex : int 1 2 2 1 2 2 1 2 2 2 ...
## $ ageyr : int 13 13 13 13 12 14 12 12 13 12 ...
## $ agemo : int 1 7 1 2 2 1 1 2 0 5 ...
## $ school: Factor w/ 2 levels "Grant-White",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ grade : int 7 7 7 7 7 7 7 7 7 7 ...
## $ x1 : num 3.33 5.33 4.5 5.33 4.83 ...
## $ x2 : num 7.75 5.25 5.25 7.75 4.75 5 6 6.25 5.75 5.25 ...
## $ x3 : num 0.375 2.125 1.875 3 0.875 ...
## $ x4 : num 2.33 1.67 1 2.67 2.67 ...
## $ x5 : num 5.75 3 1.75 4.5 4 3 6 4.25 5.75 5 ...
## $ x6 : num 1.286 1.286 0.429 2.429 2.571 ...
## $ x7 : num 3.39 3.78 3.26 3 3.7 ...
## $ x8 : num 5.75 6.25 3.9 5.3 6.3 6.65 6.2 5.15 4.65 4.55 ...
## $ x9 : num 6.36 7.92 4.42 4.86 5.92 ...Проверим, полные ли данные
data = data %>% select(x1:x9)
library(dplyr)
colSums(is.na(data))## x1 x2 x3 x4 x5 x6 x7 x8 x9
## 0 0 0 0 0 0 0 0 0Оценим распределение данных
library(psych)
formattable(describe(data))| vars | n | mean | sd | median | trimmed | mad | min | max | range | skew | kurtosis | se | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| x1 | 1 | 301 | 4.935770 | 1.167432 | 5.000000 | 4.964730 | 1.235500 | 0.6666667 | 8.500000 | 7.833333 | -0.2543455 | 0.30753382 | 0.06728967 |
| x2 | 2 | 301 | 6.088040 | 1.177451 | 6.000000 | 6.017635 | 1.111950 | 2.2500000 | 9.250000 | 7.000000 | 0.4700766 | 0.33239397 | 0.06786712 |
| x3 | 3 | 301 | 2.250415 | 1.130979 | 2.125000 | 2.199170 | 1.297275 | 0.2500000 | 4.500000 | 4.250000 | 0.3834294 | -0.90752645 | 0.06518857 |
| x4 | 4 | 301 | 3.060908 | 1.164116 | 3.000000 | 3.024896 | 0.988400 | 0.0000000 | 6.333333 | 6.333333 | 0.2674867 | 0.08012676 | 0.06709855 |
| x5 | 5 | 301 | 4.340532 | 1.290472 | 4.500000 | 4.395228 | 1.482600 | 1.0000000 | 7.000000 | 6.000000 | -0.3497961 | -0.55253689 | 0.07438158 |
| x6 | 6 | 301 | 2.185572 | 1.095603 | 2.000000 | 2.088322 | 1.059000 | 0.1428571 | 6.142857 | 6.000000 | 0.8579486 | 0.81655717 | 0.06314952 |
| x7 | 7 | 301 | 4.185902 | 1.089534 | 4.086957 | 4.163630 | 1.095835 | 1.3043478 | 7.434783 | 6.130435 | 0.2490881 | -0.30740386 | 0.06279967 |
| x8 | 8 | 301 | 5.527076 | 1.012615 | 5.500000 | 5.492946 | 0.963690 | 3.0500000 | 10.000000 | 6.950000 | 0.5252580 | 1.17155564 | 0.05836617 |
| x9 | 9 | 301 | 5.374123 | 1.009152 | 5.416667 | 5.366067 | 0.988400 | 2.7777778 | 9.250000 | 6.472222 | 0.2038709 | 0.28990791 | 0.05816654 |
Данные полные, в распределениях нет ничего, что вызывало бы подозрения
Эксплораторный факторный анализ
Используется, когда мы не определили факторную структуру, и пытаемся “нащупать” её в данных (у нас есть много вопросов и мы пытаемся найти среди них факторы). Иногда используется как пререквезит КФА (хотя некоторые исследователи считают, что для этих целей лучше подходит анализ главных компонент (PCA)).
Основные шаги:
Определить, что данные подходят для анализа
Определить количество факторов, которые есть в данных
Провести анализ
Оценить надежность получившихся факторов
Проверка того, насколько данные подходят для ЭФА
ЭФА опирается на вычисление собственных векторов матрицы:
X = ΛF+e
где:
X — матрица наблюдаемых данных,
Λ — матрица факторных нагрузок
F — матрица латентных факторов
e — матрица уникальных ошибок (шум)
Итоговая модель чаще всего основывается на методе максимального правдоподобия (maximum likelihood, ML)
Матрица корреляций
Важно: в зависимости от типа данных, могут быть необходимы различные корреляции. Для ординальных данных используется полихорическая корреляция. У нас все переменные числовые, поэтому будем пользоваться обычной корреляцией Пирсона
library(polycor)
cors = hetcor(data)
#?hetcor
cor_mat = cors$correlations
corPlot(cor_mat, type = "Pearson", numbers=T, upper=FALSE, main = "Матрица корреляций", show.legend = F)
Каждая переменная должна с кем-нибудь коррелировать
KMO
Этот индекс позволяет оценить насколько данные подходят для факторного анализа. Он проверяет, коррелируют ли между собой отдельные переменные (можно ли выделить факторы)
#?KMO
KMO(data)## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = data)
## Overall MSA = 0.75
## MSA for each item =
## x1 x2 x3 x4 x5 x6 x7 x8 x9
## 0.81 0.78 0.73 0.76 0.74 0.81 0.59 0.68 0.79Оценка MSA для каждой переменной (measure of sampling adequacy)
0.8 - очень здорово
- < 0.5 - неприемлемо
В случае, если MSA переменной низкий, можно исключить переменную её анализа. У нас результаты приемлемы (хотя х7 ведёт себя не очень хорошо)
Тест Бартлетта
Проверяет, отличается ли наша матрица корреляций от матрицы идентичности (в которой ни одна переменная не коррелирует ни с одной другой). Для ФА хорошо, когда он стат. значим
cortest.bartlett(data)## $chisq
## [1] 904.0971
##
## $p.value
## [1] 1.912079e-166
##
## $df
## [1] 36В целом, данные неплохо подходят для ЭФА
Определение количества факторов
Дальше, нам нужно понять, сколько факторов мы можем выделить из наших 9 переменных
Есть 3 основных метода:
Правило Кайзера
Основная идея: если собственное значение фактора меньше 1, то дисперсия, объясняемая фактором меньше, чем дисперсия, объясняемая одной переменной (Brown, 2015)
library(nFactors)
ev = eigen(cor(data))
ev$values## [1] 3.2163442 1.6387132 1.3651593 0.6989185 0.5843475 0.4996872 0.4731021
## [8] 0.2860024 0.23772573 фактора
Scree plot
nS <- nScree(x=ev$values, model = 'factors')
plotnScree(nS) 
3 фактора
Parallel analysis
#?fa.parallel
fa.parallel(data, fa = 'fa', fm = "ml")
## Parallel analysis suggests that the number of factors = 3 and the number of components = NA3 фактора
Можете заметить, что тут при анализе мы задали метод экстракции факторов: ml (максимальное правдоподобие)
Важно: бывает так, что факторов выделяется больше, чем мы можем ожидать. Иногда это обусловлено тем, что ошибки определенных факторов коррелируют (например, такое часто бывает с вопросами, которые заданы негативно, негативные вопросы выделяются в отдельный фактор, тогда как с точки зрения теории могли бы быть объединены с другим). Однако, если вы получили не ту факторную структуру, которую ожидали - это не всегда ошибка
Итог: Мы будем выделять 3 фактора
Анализ
Чтобы извлечь факторы из наших переменных, нужно ответить на 4 главных вопроса:
какой метод использовать? (чаще всего используют ml или pf (principal factor))
какой тип корреляции использовать для наших данных? (у нас все переменные непрерывные, поэтому используем Пирсона)
сколько факторов извлекать? (мы извлекаем 3)
какой тип ротации использовать? (если получилось более одного фактора)
Ротация используется, чтобы повысить объясняемость факторов и повысить нагрузки (корреляции переменных с факторами). Опять же идейно возвращаемся к собственным векторам :)
Ротация бывает трёх видов:
orthogonal (ортогональная) - если мы предполагаем, что полученные факторы НЕ БУДУТ коррелировать между собой (чаще всего используют метод varimax)
oblique (наклонная) - если мы предполагаем, что полученные факторы БУДУТ коррелировать между собой (чаще всего используют метод oblimin)
отсутствие ротации
Судя по переменным на которых мы работаем, можно ожидать, что полученные факторы будут коррелировать (они, вероятно, оценивают разные “измерения” интеллекта - пространственное мышление, вербальные навыки и арифметические способности), поэтому будем приоритетно пользоваться наклонной ротацией и методом oblimin
В идеальном мире можно попробовать сравнить как 3 разных метода работают на ваших данных (хотя Brown отмечает, что в социальных науках использование наклонных ротаций приоритетно)
#?fa
fa1 = fa(data, rotate = "oblimin", fm = "ml", nfactors = 3, cor = "cor")
fa1## Factor Analysis using method = ml
## Call: fa(r = data, nfactors = 3, rotate = "oblimin", fm = "ml", cor = "cor")
## Standardized loadings (pattern matrix) based upon correlation matrix
## ML1 ML3 ML2 h2 u2 com
## x1 0.19 0.60 0.03 0.49 0.51 1.2
## x2 0.04 0.51 -0.12 0.25 0.75 1.1
## x3 -0.07 0.69 0.02 0.46 0.54 1.0
## x4 0.84 0.02 0.01 0.72 0.28 1.0
## x5 0.89 -0.07 0.01 0.76 0.24 1.0
## x6 0.81 0.08 -0.01 0.69 0.31 1.0
## x7 0.04 -0.15 0.72 0.50 0.50 1.1
## x8 -0.03 0.10 0.70 0.53 0.47 1.0
## x9 0.03 0.37 0.46 0.46 0.54 1.9
##
## ML1 ML3 ML2
## SS loadings 2.24 1.34 1.28
## Proportion Var 0.25 0.15 0.14
## Cumulative Var 0.25 0.40 0.54
## Proportion Explained 0.46 0.28 0.26
## Cumulative Proportion 0.46 0.74 1.00
##
## With factor correlations of
## ML1 ML3 ML2
## ML1 1.00 0.33 0.22
## ML3 0.33 1.00 0.27
## ML2 0.22 0.27 1.00
##
## Mean item complexity = 1.2
## Test of the hypothesis that 3 factors are sufficient.
##
## df null model = 36 with the objective function = 3.05 with Chi Square = 904.1
## df of the model are 12 and the objective function was 0.08
##
## The root mean square of the residuals (RMSR) is 0.02
## The df corrected root mean square of the residuals is 0.03
##
## The harmonic n.obs is 301 with the empirical chi square 8.03 with prob < 0.78
## The total n.obs was 301 with Likelihood Chi Square = 22.38 with prob < 0.034
##
## Tucker Lewis Index of factoring reliability = 0.964
## RMSEA index = 0.053 and the 90 % confidence intervals are 0.015 0.088
## BIC = -46.11
## Fit based upon off diagonal values = 1
## Measures of factor score adequacy
## ML1 ML3 ML2
## Correlation of (regression) scores with factors 0.94 0.84 0.85
## Multiple R square of scores with factors 0.89 0.71 0.72
## Minimum correlation of possible factor scores 0.78 0.42 0.45Совместно факторы объясняют 54% дисперсии в данных
- Оценка модели:
Хи-квадрат модели должен быть незначимым (хотя на больших данных такое случается очень редко, если так случилось - скорее всего модель очень хорошая) - у нас значимый
RMSR должен быть < 0.08 - у нас 0.02
RMSEA должен быть < 0.08 (с доверительными интервалами) - у нас 0.088 по верхней границе интервала
Индекс Такера-Льюиса должен быть > 0.9 - у нас 0.96
Модель неплохая, но не идеальная
2.Значения ML1, ML2 и ML3 показывают степень корреляции между переменной и фактором (факторную нагрузку). В идеале, это корреляция должна быть высокой у переменной и одного из факторов (идеально >0.7, приемлемо > 0.4), при этом кросс-нагрузки должны быть низкими
В нашем случае хорошо ведут себя все переменные кроме х9, у которой низкая нагрузка с фактором ML2, дополняемая высокой кросс-нагрузкой на фактор ML3
h2 означает общность (communality) - это процент дисперсии переменной, объясненной всеми факторами, ну или что-то похожее на значение true score из КТТ :) (идеально, если это значение больше 0.5)
u2 означает уникальность (uniqueness), ну или что-то вроде значения ошибки из CTT :)
h2 + u2 = 1
График
Мы можем построить график факторных нагрузок для всех факторов
factor.plot(fa1)
Факторная диаграмма
График отображает корреляции (нагрузки) между факторами и переменными, а также межфакторную корреляцию
У нас получилась следующая факторная структура:
Фактор 1 (вербальные навыки):
x5 (тест 7) — Тест завершения предложений (подбор слов по контексту)
x4 (тест 6) — Тест понимания текста (ответы на вопросы по прочитанному абзацу)
x6 (тест 9) — Тест значения слов (поиск синонимов/антонимов)
Фактор 2 (пространственные навыки):
x3 (тест 4) — Тест «Лозенги» Thorndike (идентификация перевёрнутых фигур среди похожих)
x1 (тест 1) — Тест зрительного восприятия Spearman VPT (оценка анализа геометрических паттернов)
x2 (тест 2) — Упрощённый тест пространственных отношений Brigham (сборка кубов по образцу)
Фактор 3 (арифметические навыки??):
x7 (тест 10) — Тест скоростного сложения (быстрые арифметические вычисления)
x8 (тест 12) — Тест счёта точек (подсчёт точек внутри фигур на скорость)
x9 (тест 13) — Тест дискриминации линий (разделение прямых и кривых линий за время)
fa.diagram(fa1)
Анализ надежности (reliability tests)
Для оценки надёжности в ЭФА часто используют Альфу Кронбаха. Значения выше 0.7 считают приемлемыми, это будет значить, что фактор получился надежным
Для вербальных навыков:
verbal = data %>% select(x4, x5, x6)
psych::alpha(verbal)##
## Reliability analysis
## Call: psych::alpha(x = verbal)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.88 0.88 0.84 0.72 7.7 0.011 3.2 1.1 0.72
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.86 0.88 0.90
## Duhachek 0.86 0.88 0.91
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## x4 0.83 0.84 0.72 0.72 5.1 0.019 NA 0.72
## x5 0.83 0.83 0.70 0.70 4.8 0.020 NA 0.70
## x6 0.84 0.85 0.73 0.73 5.5 0.018 NA 0.73
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## x4 301 0.90 0.90 0.82 0.78 3.1 1.2
## x5 301 0.92 0.91 0.84 0.79 4.3 1.3
## x6 301 0.89 0.90 0.81 0.77 2.2 1.1Для пространственных навыков:
space = data %>% select(x1, x2, x3)
psych::alpha(space)##
## Reliability analysis
## Call: psych::alpha(x = space)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.63 0.63 0.54 0.36 1.7 0.037 4.4 0.88 0.34
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.55 0.63 0.69
## Duhachek 0.55 0.63 0.70
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## x1 0.51 0.51 0.34 0.34 1.03 0.057 NA 0.34
## x2 0.61 0.61 0.44 0.44 1.58 0.045 NA 0.44
## x3 0.46 0.46 0.30 0.30 0.85 0.062 NA 0.30
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## x1 301 0.77 0.77 0.58 0.45 4.9 1.2
## x2 301 0.73 0.72 0.47 0.37 6.1 1.2
## x3 301 0.78 0.78 0.62 0.48 2.3 1.1Для арифметических навыков:
arif = data %>% select(x7, x8, x9)
psych::alpha(arif)##
## Reliability analysis
## Call: psych::alpha(x = arif)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.69 0.69 0.6 0.43 2.2 0.031 5 0.81 0.45
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.62 0.69 0.74
## Duhachek 0.63 0.69 0.75
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## x7 0.62 0.62 0.45 0.45 1.6 0.044 NA 0.45
## x8 0.51 0.51 0.34 0.34 1.0 0.057 NA 0.34
## x9 0.65 0.65 0.49 0.49 1.9 0.040 NA 0.49
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## x7 301 0.79 0.78 0.59 0.49 4.2 1.1
## x8 301 0.82 0.82 0.69 0.57 5.5 1.0
## x9 301 0.75 0.76 0.55 0.46 5.4 1.0Однако, сам Кронбах отмечал, что использование альфы не лучшее решение:
“The numerous citations to my paper by no means indicate that the person who cited it had read it, and does not even demonstrate that he had looked at it” (Cronbach, L. J., & Shavelson, R. J. (2004). My current thoughts on coefficient alpha and successor procedures. Educational and Psychological Measurement, 64, 391–418. doi:10.1177/0013164404266386
Альтернативой альфе может быть омега:
Dunn, T. J., Baguley, T., & Brunsden, V. (2014). From alpha to omega: A practical solution to the pervasive problem of internal consistency estimation. British journal of psychology, 105(3), 399-412.
omega(data, rotate = "oblimin", fm = "ml", nfactors = 3, cor = "cor")
## Omega
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip,
## digits = digits, title = title, sl = sl, labels = labels,
## plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option,
## covar = covar, cor = "cor")
## Alpha: 0.76
## G.6: 0.81
## Omega Hierarchical: 0.45
## Omega H asymptotic: 0.53
## Omega Total 0.85
##
## Schmid Leiman Factor loadings greater than 0.2
## g F1* F2* F3* h2 u2 p2
## x1 0.49 0.46 0.49 0.51 0.50
## x2 0.30 0.39 0.25 0.75 0.35
## x3 0.41 0.53 0.46 0.54 0.38
## x4 0.45 0.72 0.72 0.28 0.28
## x5 0.41 0.76 0.76 0.24 0.23
## x6 0.46 0.69 0.69 0.31 0.30
## x7 0.23 0.65 0.50 0.50 0.11
## x8 0.35 0.64 0.53 0.47 0.23
## x9 0.45 0.28 0.42 0.46 0.54 0.44
##
## With Sums of squares of:
## g F1* F2* F3*
## 1.46 1.62 0.75 1.02
##
## general/max 0.9 max/min = 2.15
## mean percent general = 0.31 with sd = 0.12 and cv of 0.38
## Explained Common Variance of the general factor = 0.3
##
## The degrees of freedom are 12 and the fit is 0.08
## The number of observations was 301 with Chi Square = 22.38 with prob < 0.034
## The root mean square of the residuals is 0.02
## The df corrected root mean square of the residuals is 0.03
## RMSEA index = 0.053 and the 10 % confidence intervals are 0.015 0.088
## BIC = -46.11
##
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 27 and the fit is 1.75
## The number of observations was 301 with Chi Square = 517.19 with prob < 4.4e-92
## The root mean square of the residuals is 0.2
## The df corrected root mean square of the residuals is 0.23
##
## RMSEA index = 0.246 and the 10 % confidence intervals are 0.228 0.265
## BIC = 363.1
##
## Measures of factor score adequacy
## g F1* F2* F3*
## Correlation of scores with factors 0.68 0.84 0.67 0.78
## Multiple R square of scores with factors 0.46 0.71 0.44 0.61
## Minimum correlation of factor score estimates -0.07 0.43 -0.11 0.22
##
## Total, General and Subset omega for each subset
## g F1* F2* F3*
## Omega total for total scores and subscales 0.85 0.89 0.65 0.72
## Omega general for total scores and subscales 0.45 0.24 0.28 0.19
## Omega group for total scores and subscales 0.35 0.65 0.37 0.53Извлечение факторных значений
Основываясь на результатах анализа мы можем сказать, что он более-менее приемлем (хотя и далеко не идеален). Теперь мы можем попробовать извлечь факторные значения для дальнейшего анализа (т.е. попробовать оценить значение латентной переменной для каждого наблюдения)
Есть 2 основных способа (Brown, 2015):
грубые (coarse) оценки (суммирование или усреднение значений переменных, относящихся к фактору, часто получаются biased)
очищенные (refined) оценки (могут быть извлечены разными методами, менее biased)
Попробуем извлечь очищенные оценки для каждого из факторов (по умолчанию они вычисляются с помощью простой регрессии, но если углубиться в тему, то можно выбрать метод, который лучше подходит по вашему мнению)
data_with_scores = HolzingerSwineford1939
data_with_scores$verbal = as.data.frame(fa1$scores)$ML1
data_with_scores$space = as.data.frame(fa1$scores)$ML2
data_with_scores$arif = as.data.frame(fa1$scores)$ML3
#получаем примерно стандартное нормальное распределение
library(ggpubr)
hist(data_with_scores$verbal)
ggqqplot(data_with_scores$verbal)
Проведём статистический тест. Отличаются ли вербальные способности (этот фактор получился наиболее хорошим) у представителей разного пола?
С вашего позволения проведём тест без предварительной проверки на assumptions.
t.test(verbal ~ sex, data_with_scores)##
## Welch Two Sample t-test
##
## data: verbal by sex
## t = -1.1887, df = 298.81, p-value = 0.2355
## alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
## 95 percent confidence interval:
## -0.34245396 0.08453608
## sample estimates:
## mean in group 1 mean in group 2
## -0.06640743 0.06255151Вербальные способности у представителей разного пола на наших данных не отличаются
Конфирматорный факторный анализ
Используется, когда мы на основе собственных размышлений самостоятельно определили факторную структуру. Таким образом, проверяется насколько данные подтверждают наши предположения о существующей факторной структуре. Часто является первым шагом при моделировании структурными уравнениями (SEM)
При проверке каких-либо шкал ЭФА и КФА иногда используют совместно на одних и тех же данных, хотя и по отдельности эти методы очень полезны.
Важно: факторный анализ - непростая техника в ряде случаев. Порой результаты проводимых на одних и тех же данных ЭФА и КФА не совпадают, что может послужить очень хорошей пищей для размышлений для исследователя о существующей факторной структуре в данных
Основные шаги:
Подготовить данные и провести их разведку
Теоретически обосновать предполагаемую факторную структуру
Провести анализ (общая факторная структура и каждый фактор отдельно)
2.1 Доработать модель (исключить переменные или использовать рекомендации из индексов модицифации)
- Оценить надежность получившихся факторов
Теоретическое обоснование модели
Почитав литературу по теме и пользуясь логикой здравого смысла мы решили определить следующую факторную структуру:
Holzinger, K. J., & Swineford, F. (1939). A study in factor analysis: The stability of a bi-factor solution. Supplementary educational monographs.
Фактор 1 (визуальные навыки):
x1 (тест 1) — Тест зрительного восприятия Spearman VPT (оценка анализа геометрических паттернов)
x2 (тест 2) — Упрощённый тест пространственных отношений Brigham (сборка кубов по образцу)
x3 (тест 4) — Тест «Лозенги» Thorndike (идентификация перевёрнутых фигур среди похожих)
Фактор 2 (текстовые навыки):
x4 (тест 6) — Тест понимания текста (ответы на вопросы по прочитанному абзацу)
x5 (тест 7) — Тест завершения предложений (подбор слов по контексту)
x6 (тест 9) — Тест значения слов (поиск синонимов/антонимов)
Фактор 3 (навыки скорости):
x7 (тест 10) — Тест скоростного сложения (быстрые арифметические вычисления)
x8 (тест 12) — Тест счёта точек (подсчёт точек внутри фигур на скорость)
x9 (тест 13) — Тест дискриминации линий (разделение прямых и кривых линий за время)
Анализ
Факторную структуру мы определили сами, поэтому мы знаем сколько факторов хотим извлекать и какие переменные в них входят, поэтому можно приступать к анализу сразу.
Важно: если у вас ординальные данные, добавьте в функцию cfa аргумент ordered = True. Кроме того, в функции summary мы добавили аргумент standardized = TRUE, благодаря чему все значения будут стандартизированы (напр., все ковариации превратятся в корреляции)
Общая модель
library(lavaan)
library(semTools)
model = '
visual =~ x1 + x2 + x3
textual =~ x4 + x5 + x6
speed =~ x7 + x8 + x9
'
fit_first = cfa(model, data)
summary(fit_first, fit.measures = TRUE, standardized = TRUE, rsquare = T)## lavaan 0.6-19 ended normally after 35 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 21
##
## Number of observations 301
##
## Model Test User Model:
##
## Test statistic 85.306
## Degrees of freedom 24
## P-value (Chi-square) 0.000
##
## Model Test Baseline Model:
##
## Test statistic 918.852
## Degrees of freedom 36
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.931
## Tucker-Lewis Index (TLI) 0.896
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -3737.745
## Loglikelihood unrestricted model (H1) -3695.092
##
## Akaike (AIC) 7517.490
## Bayesian (BIC) 7595.339
## Sample-size adjusted Bayesian (SABIC) 7528.739
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.092
## 90 Percent confidence interval - lower 0.071
## 90 Percent confidence interval - upper 0.114
## P-value H_0: RMSEA <= 0.050 0.001
## P-value H_0: RMSEA >= 0.080 0.840
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.065
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## visual =~
## x1 1.000 0.900 0.772
## x2 0.554 0.100 5.554 0.000 0.498 0.424
## x3 0.729 0.109 6.685 0.000 0.656 0.581
## textual =~
## x4 1.000 0.990 0.852
## x5 1.113 0.065 17.014 0.000 1.102 0.855
## x6 0.926 0.055 16.703 0.000 0.917 0.838
## speed =~
## x7 1.000 0.619 0.570
## x8 1.180 0.165 7.152 0.000 0.731 0.723
## x9 1.082 0.151 7.155 0.000 0.670 0.665
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## visual ~~
## textual 0.408 0.074 5.552 0.000 0.459 0.459
## speed 0.262 0.056 4.660 0.000 0.471 0.471
## textual ~~
## speed 0.173 0.049 3.518 0.000 0.283 0.283
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .x1 0.549 0.114 4.833 0.000 0.549 0.404
## .x2 1.134 0.102 11.146 0.000 1.134 0.821
## .x3 0.844 0.091 9.317 0.000 0.844 0.662
## .x4 0.371 0.048 7.779 0.000 0.371 0.275
## .x5 0.446 0.058 7.642 0.000 0.446 0.269
## .x6 0.356 0.043 8.277 0.000 0.356 0.298
## .x7 0.799 0.081 9.823 0.000 0.799 0.676
## .x8 0.488 0.074 6.573 0.000 0.488 0.477
## .x9 0.566 0.071 8.003 0.000 0.566 0.558
## visual 0.809 0.145 5.564 0.000 1.000 1.000
## textual 0.979 0.112 8.737 0.000 1.000 1.000
## speed 0.384 0.086 4.451 0.000 1.000 1.000
##
## R-Square:
## Estimate
## x1 0.596
## x2 0.179
## x3 0.338
## x4 0.725
## x5 0.731
## x6 0.702
## x7 0.324
## x8 0.523
## x9 0.442- Оценка модели, global fit (почти такая же, как в ЭФА):
Хи-квадрат модели должен быть незначимым (хотя на больших данных такое случается очень редко, если так случилось - скорее всего модель очень хорошая) - у нас значимый
SRMR должен быть < 0.08 - у нас 0.065
RMSEA должен быть < 0.08 (с доверительными интервалами) - у нас 0.114 по верхней границе интервала
Индекс Такера-Льюиса (TLI) должен быть > 0.9 - у нас 0.896
Comparative fit index (CFI) должен быть > 0.9 - у нас 0.931
TLI и CFI очень похожи по своей сути, но часто TLI можно считать более приоритетным, так как при расчетах он также оценивает сложность модели и “штрафует” более сложные
Модель неплохая, но далеко не идеальная
2.Значения std.all в latent variables показывает факторную нагрузку (в идеале должна быть у каждой переменной больше 0,7)
R-квадрат показывает долю дисперсии, объясненную фактором внутри той или иной переменной (в идеале больше 0,5)
Факторная структура
library(lavaanPlot)
lavaanPlot(model = fit_first,
stand = TRUE, coefs = T)Далее, перейдём к оценке того, как мы можем улучшить нашу модель
Ковариации остатков
Ковариация остатков показывает степень взаимосвязанности остатков (ошибок) между собой. Считается, что ковариации по модулю превосходящие 0.1 должны вызывать опасения (часто это ещё может быть сопряжено с высоким RSMR).
В идеале все ковариации по модулю должны быть менее 0.1
resid_cov = lavResiduals(fit_first)$cov
(resid_cov > 0.1) | (resid_cov < -0.1)## x1 x2 x3 x4 x5 x6 x7 x8 x9
## x1 FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE TRUE
## x2 FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## x3 FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE TRUE
## x4 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## x5 FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE
## x6 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## x7 TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## x8 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## x9 TRUE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSEИндексы модификации
Индексы модицикации показывают на какое значение уменьшится хи-квадрат модели, если применить соответствующую модификацию. MI больше 10 считаются важными
В случае, если какая-то переменная встречается часто среди MI, можно попробовать её исключить из модели (мы так сделать не можем, у нас по 3 переменных на фактор)
Либо мы можем попробовать использовать рекомендацию ОДНОГО из индексов за один раз, после чего следует заново оценить качество модели
Важно: MI нельзя использовать “просто так”. Они должны быть подкреплены какой-либо теоретической предпосылкой. И если вы используете эту предпосылку один раз, то её нужно использовать везде (например, если у двух негативно заданных переменных коррелируют остатки, и вы решили учесть это в модели, то учесть это лучше не в одном случае, а во всех)
modindices(fit_first, sort. = T)## lhs op rhs mi epc sepc.lv sepc.all sepc.nox
## 30 visual =~ x9 36.411 0.577 0.519 0.515 0.515
## 76 x7 ~~ x8 34.145 0.536 0.536 0.859 0.859
## 28 visual =~ x7 18.631 -0.422 -0.380 -0.349 -0.349
## 78 x8 ~~ x9 14.946 -0.423 -0.423 -0.805 -0.805
## 33 textual =~ x3 9.151 -0.272 -0.269 -0.238 -0.238
## 55 x2 ~~ x7 8.918 -0.183 -0.183 -0.192 -0.192
## 31 textual =~ x1 8.903 0.350 0.347 0.297 0.297
## 51 x2 ~~ x3 8.532 0.218 0.218 0.223 0.223
## 59 x3 ~~ x5 7.858 -0.130 -0.130 -0.212 -0.212
## 26 visual =~ x5 7.441 -0.210 -0.189 -0.147 -0.147
## 50 x1 ~~ x9 7.335 0.138 0.138 0.247 0.247
## 65 x4 ~~ x6 6.220 -0.235 -0.235 -0.646 -0.646
## 66 x4 ~~ x7 5.920 0.098 0.098 0.180 0.180
## 48 x1 ~~ x7 5.420 -0.129 -0.129 -0.195 -0.195
## 77 x7 ~~ x9 5.183 -0.187 -0.187 -0.278 -0.278
## 36 textual =~ x9 4.796 0.138 0.137 0.136 0.136
## 29 visual =~ x8 4.295 -0.210 -0.189 -0.187 -0.187
## 63 x3 ~~ x9 4.126 0.102 0.102 0.147 0.147
## 67 x4 ~~ x8 3.805 -0.069 -0.069 -0.162 -0.162
## 43 x1 ~~ x2 3.606 -0.184 -0.184 -0.233 -0.233
## 45 x1 ~~ x4 3.554 0.078 0.078 0.173 0.173
## 35 textual =~ x8 3.359 -0.121 -0.120 -0.118 -0.118
## 27 visual =~ x6 2.843 0.111 0.100 0.092 0.092
## 64 x4 ~~ x5 2.534 0.186 0.186 0.457 0.457
## 57 x2 ~~ x9 1.895 0.075 0.075 0.094 0.094
## 60 x3 ~~ x6 1.855 0.055 0.055 0.100 0.100
## 38 speed =~ x2 1.580 -0.198 -0.123 -0.105 -0.105
## 70 x5 ~~ x7 1.233 -0.049 -0.049 -0.083 -0.083
## 25 visual =~ x4 1.211 0.077 0.069 0.059 0.059
## 72 x5 ~~ x9 0.999 0.040 0.040 0.079 0.079
## 44 x1 ~~ x3 0.935 -0.139 -0.139 -0.203 -0.203
## 69 x5 ~~ x6 0.916 0.101 0.101 0.253 0.253
## 54 x2 ~~ x6 0.785 0.039 0.039 0.062 0.062
## 39 speed =~ x3 0.716 0.136 0.084 0.075 0.075
## 61 x3 ~~ x7 0.638 -0.044 -0.044 -0.054 -0.054
## 49 x1 ~~ x8 0.634 -0.041 -0.041 -0.079 -0.079
## 52 x2 ~~ x4 0.534 -0.034 -0.034 -0.052 -0.052
## 46 x1 ~~ x5 0.522 -0.033 -0.033 -0.067 -0.067
## 71 x5 ~~ x8 0.347 0.023 0.023 0.049 0.049
## 74 x6 ~~ x8 0.275 0.018 0.018 0.043 0.043
## 42 speed =~ x6 0.273 0.044 0.027 0.025 0.025
## 73 x6 ~~ x7 0.259 -0.020 -0.020 -0.037 -0.037
## 41 speed =~ x5 0.201 -0.044 -0.027 -0.021 -0.021
## 68 x4 ~~ x9 0.196 -0.016 -0.016 -0.035 -0.035
## 58 x3 ~~ x4 0.142 -0.016 -0.016 -0.028 -0.028
## 34 textual =~ x7 0.098 -0.021 -0.021 -0.019 -0.019
## 75 x6 ~~ x9 0.097 -0.011 -0.011 -0.024 -0.024
## 62 x3 ~~ x8 0.059 -0.012 -0.012 -0.019 -0.019
## 56 x2 ~~ x8 0.054 -0.012 -0.012 -0.017 -0.017
## 47 x1 ~~ x6 0.048 0.009 0.009 0.020 0.020
## 53 x2 ~~ x5 0.023 -0.008 -0.008 -0.011 -0.011
## 32 textual =~ x2 0.017 -0.011 -0.011 -0.010 -0.010
## 37 speed =~ x1 0.014 0.024 0.015 0.013 0.013
## 40 speed =~ x4 0.003 -0.005 -0.003 -0.003 -0.003Далее нам нужно провести CFA не только для общей факторной структуры, но и для каждого фактора в отдельности (это может помочь улучшить общую модель, найти скрытые в отдельных факторах проблемы)
Обратите внимание: так как у нас только 3 переменных в каждом факторе, то все однофакторные модели будут just-identified и мы не сможем оценить local fit (SRMR и RMSEA равны 0, CFI и TLI равны 1)
Если бы в каком-то из факторов у нас было больше 4 переменных, мы могли бы потенциально использовать рекомендации MI для отдельного фактора
Визуальные навыки
library(lavaan)
library(semTools)
model = '
visual =~ x1 + x2 + x3
'
fit = cfa(model, data)
summary(fit, fit.measures = TRUE, standardized = TRUE, rsquare = T)## lavaan 0.6-19 ended normally after 23 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 6
##
## Number of observations 301
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 0
##
## Model Test Baseline Model:
##
## Test statistic 111.271
## Degrees of freedom 3
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 1.000
## Tucker-Lewis Index (TLI) 1.000
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1356.977
## Loglikelihood unrestricted model (H1) -1356.977
##
## Akaike (AIC) 2725.955
## Bayesian (BIC) 2748.197
## Sample-size adjusted Bayesian (SABIC) 2729.169
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.000
## P-value H_0: RMSEA <= 0.050 NA
## P-value H_0: RMSEA >= 0.080 NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.000
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## visual =~
## x1 1.000 0.724 0.621
## x2 0.778 0.141 5.532 0.000 0.563 0.479
## x3 1.107 0.214 5.173 0.000 0.801 0.710
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .x1 0.835 0.118 7.064 0.000 0.835 0.614
## .x2 1.065 0.105 10.177 0.000 1.065 0.771
## .x3 0.633 0.129 4.899 0.000 0.633 0.496
## visual 0.524 0.130 4.021 0.000 1.000 1.000
##
## R-Square:
## Estimate
## x1 0.386
## x2 0.229
## x3 0.504library(lavaanPlot)
lavaanPlot(model = fit,
stand = TRUE, coefs = T)Не имеет смысла, так как модель just-identified и SRMR = 0
resid_cov = lavResiduals(fit)$cov
(resid_cov > 0.1) | (resid_cov < -0.1)## x1 x2 x3
## x1 FALSE FALSE FALSE
## x2 FALSE FALSE FALSE
## x3 FALSE FALSE FALSEНе будет работать, так как модель just-identified
modindices(fit, sort. = T)## [1] lhs op rhs mi epc sepc.lv sepc.all sepc.nox
## <0 строк> (или 'row.names' нулевой длины)Текстовые навыки
library(lavaan)
library(semTools)
model = '
textual =~ x4 + x5 + x6
'
fit = cfa(model, data)
summary(fit, fit.measures = TRUE, standardized = TRUE, rsquare = T)## lavaan 0.6-19 ended normally after 15 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 6
##
## Number of observations 301
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 0
##
## Model Test Baseline Model:
##
## Test statistic 497.430
## Degrees of freedom 3
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 1.000
## Tucker-Lewis Index (TLI) 1.000
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1181.065
## Loglikelihood unrestricted model (H1) -1181.065
##
## Akaike (AIC) 2374.130
## Bayesian (BIC) 2396.372
## Sample-size adjusted Bayesian (SABIC) 2377.344
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.000
## P-value H_0: RMSEA <= 0.050 NA
## P-value H_0: RMSEA >= 0.080 NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.000
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## textual =~
## x4 1.000 0.984 0.847
## x5 1.133 0.067 16.906 0.000 1.115 0.866
## x6 0.924 0.056 16.391 0.000 0.910 0.832
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .x4 0.382 0.049 7.805 0.000 0.382 0.283
## .x5 0.416 0.059 7.038 0.000 0.416 0.251
## .x6 0.369 0.044 8.367 0.000 0.369 0.308
## textual 0.969 0.112 8.640 0.000 1.000 1.000
##
## R-Square:
## Estimate
## x4 0.717
## x5 0.749
## x6 0.692library(lavaanPlot)
lavaanPlot(model = fit,
stand = TRUE, coefs = T)Навыки скорости
library(lavaan)
library(semTools)
model = '
speed =~ x7 + x8 + x9
'
fit = cfa(model, data)
summary(fit, fit.measures = TRUE, standardized = TRUE, rsquare = T)## lavaan 0.6-19 ended normally after 21 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 6
##
## Number of observations 301
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 0
##
## Model Test Baseline Model:
##
## Test statistic 156.623
## Degrees of freedom 3
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 1.000
## Tucker-Lewis Index (TLI) 1.000
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1233.814
## Loglikelihood unrestricted model (H1) -1233.814
##
## Akaike (AIC) 2479.627
## Bayesian (BIC) 2501.870
## Sample-size adjusted Bayesian (SABIC) 2482.841
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.000
## P-value H_0: RMSEA <= 0.050 NA
## P-value H_0: RMSEA >= 0.080 NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.000
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## speed =~
## x7 1.000 0.661 0.608
## x8 1.225 0.190 6.460 0.000 0.810 0.801
## x9 0.854 0.121 7.046 0.000 0.565 0.561
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .x7 0.746 0.086 8.650 0.000 0.746 0.631
## .x8 0.366 0.097 3.794 0.000 0.366 0.358
## .x9 0.696 0.072 9.640 0.000 0.696 0.686
## speed 0.437 0.097 4.520 0.000 1.000 1.000
##
## R-Square:
## Estimate
## x7 0.369
## x8 0.642
## x9 0.314library(lavaanPlot)
lavaanPlot(model = fit,
stand = TRUE, coefs = T)Модификация модели
Оценив global и local fit модели, попробуем её улучшить
Индексы модификации нашей первой модели следующие:
modindices(fit_first, sort. = T)## lhs op rhs mi epc sepc.lv sepc.all sepc.nox
## 30 visual =~ x9 36.411 0.577 0.519 0.515 0.515
## 76 x7 ~~ x8 34.145 0.536 0.536 0.859 0.859
## 28 visual =~ x7 18.631 -0.422 -0.380 -0.349 -0.349
## 78 x8 ~~ x9 14.946 -0.423 -0.423 -0.805 -0.805
## 33 textual =~ x3 9.151 -0.272 -0.269 -0.238 -0.238
## 55 x2 ~~ x7 8.918 -0.183 -0.183 -0.192 -0.192
## 31 textual =~ x1 8.903 0.350 0.347 0.297 0.297
## 51 x2 ~~ x3 8.532 0.218 0.218 0.223 0.223
## 59 x3 ~~ x5 7.858 -0.130 -0.130 -0.212 -0.212
## 26 visual =~ x5 7.441 -0.210 -0.189 -0.147 -0.147
## 50 x1 ~~ x9 7.335 0.138 0.138 0.247 0.247
## 65 x4 ~~ x6 6.220 -0.235 -0.235 -0.646 -0.646
## 66 x4 ~~ x7 5.920 0.098 0.098 0.180 0.180
## 48 x1 ~~ x7 5.420 -0.129 -0.129 -0.195 -0.195
## 77 x7 ~~ x9 5.183 -0.187 -0.187 -0.278 -0.278
## 36 textual =~ x9 4.796 0.138 0.137 0.136 0.136
## 29 visual =~ x8 4.295 -0.210 -0.189 -0.187 -0.187
## 63 x3 ~~ x9 4.126 0.102 0.102 0.147 0.147
## 67 x4 ~~ x8 3.805 -0.069 -0.069 -0.162 -0.162
## 43 x1 ~~ x2 3.606 -0.184 -0.184 -0.233 -0.233
## 45 x1 ~~ x4 3.554 0.078 0.078 0.173 0.173
## 35 textual =~ x8 3.359 -0.121 -0.120 -0.118 -0.118
## 27 visual =~ x6 2.843 0.111 0.100 0.092 0.092
## 64 x4 ~~ x5 2.534 0.186 0.186 0.457 0.457
## 57 x2 ~~ x9 1.895 0.075 0.075 0.094 0.094
## 60 x3 ~~ x6 1.855 0.055 0.055 0.100 0.100
## 38 speed =~ x2 1.580 -0.198 -0.123 -0.105 -0.105
## 70 x5 ~~ x7 1.233 -0.049 -0.049 -0.083 -0.083
## 25 visual =~ x4 1.211 0.077 0.069 0.059 0.059
## 72 x5 ~~ x9 0.999 0.040 0.040 0.079 0.079
## 44 x1 ~~ x3 0.935 -0.139 -0.139 -0.203 -0.203
## 69 x5 ~~ x6 0.916 0.101 0.101 0.253 0.253
## 54 x2 ~~ x6 0.785 0.039 0.039 0.062 0.062
## 39 speed =~ x3 0.716 0.136 0.084 0.075 0.075
## 61 x3 ~~ x7 0.638 -0.044 -0.044 -0.054 -0.054
## 49 x1 ~~ x8 0.634 -0.041 -0.041 -0.079 -0.079
## 52 x2 ~~ x4 0.534 -0.034 -0.034 -0.052 -0.052
## 46 x1 ~~ x5 0.522 -0.033 -0.033 -0.067 -0.067
## 71 x5 ~~ x8 0.347 0.023 0.023 0.049 0.049
## 74 x6 ~~ x8 0.275 0.018 0.018 0.043 0.043
## 42 speed =~ x6 0.273 0.044 0.027 0.025 0.025
## 73 x6 ~~ x7 0.259 -0.020 -0.020 -0.037 -0.037
## 41 speed =~ x5 0.201 -0.044 -0.027 -0.021 -0.021
## 68 x4 ~~ x9 0.196 -0.016 -0.016 -0.035 -0.035
## 58 x3 ~~ x4 0.142 -0.016 -0.016 -0.028 -0.028
## 34 textual =~ x7 0.098 -0.021 -0.021 -0.019 -0.019
## 75 x6 ~~ x9 0.097 -0.011 -0.011 -0.024 -0.024
## 62 x3 ~~ x8 0.059 -0.012 -0.012 -0.019 -0.019
## 56 x2 ~~ x8 0.054 -0.012 -0.012 -0.017 -0.017
## 47 x1 ~~ x6 0.048 0.009 0.009 0.020 0.020
## 53 x2 ~~ x5 0.023 -0.008 -0.008 -0.011 -0.011
## 32 textual =~ x2 0.017 -0.011 -0.011 -0.010 -0.010
## 37 speed =~ x1 0.014 0.024 0.015 0.013 0.013
## 40 speed =~ x4 0.003 -0.005 -0.003 -0.003 -0.003Предлагаю скоррелировать тесты x7 и x8 (это сложно обосновать теоретически, так как не очень понятно в чём эти тесты заключаются, но в рамках примера давайте предположим, что у них очень-очень похожая процедура, что может привести к корреляции ошибок при проведении этих тестов)
Общая модель
Добавим взаимосвязь x7 и x8 в модель и оценим, что у нас получается в таком случае
library(lavaan)
library(semTools)
model = '
visual =~ x1 + x2 + x3
textual =~ x4 + x5 + x6
speed =~ x7 + x8 + x9
x7 ~~ x8
'
fit_final = cfa(model, data)
summary(fit_final, fit.measures = TRUE, standardized = TRUE, rsquare = T)## lavaan 0.6-19 ended normally after 43 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 22
##
## Number of observations 301
##
## Model Test User Model:
##
## Test statistic 53.272
## Degrees of freedom 23
## P-value (Chi-square) 0.000
##
## Model Test Baseline Model:
##
## Test statistic 918.852
## Degrees of freedom 36
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.966
## Tucker-Lewis Index (TLI) 0.946
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -3721.728
## Loglikelihood unrestricted model (H1) -3695.092
##
## Akaike (AIC) 7487.457
## Bayesian (BIC) 7569.013
## Sample-size adjusted Bayesian (SABIC) 7499.242
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.066
## 90 Percent confidence interval - lower 0.043
## 90 Percent confidence interval - upper 0.090
## P-value H_0: RMSEA <= 0.050 0.118
## P-value H_0: RMSEA >= 0.080 0.175
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.047
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## visual =~
## x1 1.000 0.885 0.759
## x2 0.576 0.098 5.898 0.000 0.509 0.433
## x3 0.752 0.103 7.289 0.000 0.665 0.589
## textual =~
## x4 1.000 0.989 0.851
## x5 1.115 0.066 17.015 0.000 1.103 0.856
## x6 0.926 0.056 16.682 0.000 0.916 0.838
## speed =~
## x7 1.000 0.383 0.352
## x8 1.244 0.194 6.414 0.000 0.477 0.471
## x9 2.515 0.641 3.924 0.000 0.963 0.956
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .x7 ~~
## .x8 0.353 0.067 5.239 0.000 0.353 0.389
## visual ~~
## textual 0.400 0.073 5.511 0.000 0.457 0.457
## speed 0.184 0.054 3.423 0.001 0.544 0.544
## textual ~~
## speed 0.102 0.036 2.854 0.004 0.270 0.270
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .x1 0.576 0.101 5.678 0.000 0.576 0.424
## .x2 1.122 0.100 11.171 0.000 1.122 0.812
## .x3 0.832 0.087 9.552 0.000 0.832 0.653
## .x4 0.372 0.048 7.791 0.000 0.372 0.276
## .x5 0.444 0.058 7.600 0.000 0.444 0.267
## .x6 0.357 0.043 8.287 0.000 0.357 0.298
## .x7 1.036 0.090 11.501 0.000 1.036 0.876
## .x8 0.795 0.080 9.988 0.000 0.795 0.778
## .x9 0.088 0.188 0.466 0.641 0.088 0.086
## visual 0.783 0.135 5.810 0.000 1.000 1.000
## textual 0.978 0.112 8.729 0.000 1.000 1.000
## speed 0.147 0.056 2.615 0.009 1.000 1.000
##
## R-Square:
## Estimate
## x1 0.576
## x2 0.188
## x3 0.347
## x4 0.724
## x5 0.733
## x6 0.702
## x7 0.124
## x8 0.222
## x9 0.914Модель стала лучше!!!
library(lavaanPlot)
lavaanPlot(model = fit_final,
stand = TRUE, coefs = T)resid_cov = lavResiduals(fit_final)$cov
(resid_cov > 0.1) | (resid_cov < -0.1)## x1 x2 x3 x4 x5 x6 x7 x8 x9
## x1 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## x2 FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## x3 FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE
## x4 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## x5 FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE
## x6 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## x7 FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## x8 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## x9 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSEmodindices(fit_final, sort. = T)## lhs op rhs mi epc sepc.lv sepc.all sepc.nox
## 34 textual =~ x3 9.057 -0.257 -0.254 -0.225 -0.225
## 32 textual =~ x1 8.734 0.313 0.309 0.265 0.265
## 60 x3 ~~ x5 8.712 -0.136 -0.136 -0.223 -0.223
## 67 x4 ~~ x7 8.117 0.112 0.112 0.181 0.181
## 56 x2 ~~ x7 7.266 -0.160 -0.160 -0.148 -0.148
## 52 x2 ~~ x3 6.536 0.181 0.181 0.187 0.187
## 27 visual =~ x5 6.076 -0.191 -0.169 -0.131 -0.131
## 29 visual =~ x7 5.438 -0.231 -0.204 -0.188 -0.188
## 66 x4 ~~ x6 5.367 -0.222 -0.222 -0.609 -0.609
## 30 visual =~ x8 5.049 0.276 0.244 0.242 0.242
## 78 x8 ~~ x9 4.806 -0.293 -0.293 -1.113 -1.113
## 77 x7 ~~ x9 4.806 0.236 0.236 0.783 0.783
## 46 x1 ~~ x4 4.093 0.083 0.083 0.179 0.179
## 68 x4 ~~ x8 3.615 -0.066 -0.066 -0.121 -0.121
## 44 x1 ~~ x2 2.967 -0.148 -0.148 -0.184 -0.184
## 28 visual =~ x6 2.501 0.105 0.093 0.085 0.085
## 65 x4 ~~ x5 2.276 0.180 0.180 0.443 0.443
## 73 x5 ~~ x9 2.025 0.055 0.055 0.281 0.281
## 61 x3 ~~ x6 1.994 0.056 0.056 0.103 0.103
## 49 x1 ~~ x7 1.981 -0.074 -0.074 -0.096 -0.096
## 71 x5 ~~ x7 1.374 -0.051 -0.051 -0.075 -0.075
## 50 x1 ~~ x8 1.258 0.056 0.056 0.083 0.083
## 75 x6 ~~ x8 0.916 0.032 0.032 0.060 0.060
## 37 textual =~ x9 0.914 -0.127 -0.126 -0.125 -0.125
## 31 visual =~ x9 0.914 2.819 2.494 2.475 2.475
## 26 visual =~ x4 0.892 0.066 0.058 0.050 0.050
## 76 x6 ~~ x9 0.873 -0.031 -0.031 -0.178 -0.178
## 47 x1 ~~ x5 0.859 -0.042 -0.042 -0.083 -0.083
## 55 x2 ~~ x6 0.773 0.039 0.039 0.061 0.061
## 35 textual =~ x7 0.753 0.054 0.053 0.049 0.049
## 69 x4 ~~ x9 0.751 -0.031 -0.031 -0.169 -0.169
## 70 x5 ~~ x6 0.749 0.093 0.093 0.234 0.234
## 40 speed =~ x3 0.668 0.210 0.080 0.071 0.071
## 53 x2 ~~ x4 0.521 -0.033 -0.033 -0.051 -0.051
## 72 x5 ~~ x8 0.505 0.027 0.027 0.046 0.046
## 64 x3 ~~ x9 0.475 0.039 0.039 0.145 0.145
## 57 x2 ~~ x8 0.436 0.035 0.035 0.037 0.037
## 45 x1 ~~ x3 0.434 -0.079 -0.079 -0.115 -0.115
## 39 speed =~ x2 0.373 -0.148 -0.057 -0.048 -0.048
## 63 x3 ~~ x8 0.218 0.023 0.023 0.028 0.028
## 62 x3 ~~ x7 0.195 -0.024 -0.024 -0.026 -0.026
## 48 x1 ~~ x6 0.165 0.016 0.016 0.035 0.035
## 74 x6 ~~ x7 0.137 -0.014 -0.014 -0.023 -0.023
## 51 x1 ~~ x9 0.128 -0.025 -0.025 -0.111 -0.111
## 41 speed =~ x4 0.125 -0.044 -0.017 -0.015 -0.015
## 38 speed =~ x1 0.097 -0.103 -0.039 -0.034 -0.034
## 59 x3 ~~ x4 0.089 -0.012 -0.012 -0.022 -0.022
## 54 x2 ~~ x5 0.066 -0.013 -0.013 -0.018 -0.018
## 42 speed =~ x5 0.057 0.033 0.013 0.010 0.010
## 33 textual =~ x2 0.041 -0.017 -0.017 -0.015 -0.015
## 36 textual =~ x8 0.036 0.011 0.011 0.011 0.011
## 58 x2 ~~ x9 0.015 0.007 0.007 0.022 0.022
## 43 speed =~ x6 0.013 0.014 0.005 0.005 0.005Нет ни одного MI больше 10 (после всего лишь одного маленького изменения)
Визуальные навыки
library(lavaan)
library(semTools)
model = '
visual =~ x1 + x2 + x3
'
fit = cfa(model, data)
summary(fit, fit.measures = TRUE, standardized = TRUE, rsquare = T)## lavaan 0.6-19 ended normally after 23 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 6
##
## Number of observations 301
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 0
##
## Model Test Baseline Model:
##
## Test statistic 111.271
## Degrees of freedom 3
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 1.000
## Tucker-Lewis Index (TLI) 1.000
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1356.977
## Loglikelihood unrestricted model (H1) -1356.977
##
## Akaike (AIC) 2725.955
## Bayesian (BIC) 2748.197
## Sample-size adjusted Bayesian (SABIC) 2729.169
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.000
## P-value H_0: RMSEA <= 0.050 NA
## P-value H_0: RMSEA >= 0.080 NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.000
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## visual =~
## x1 1.000 0.724 0.621
## x2 0.778 0.141 5.532 0.000 0.563 0.479
## x3 1.107 0.214 5.173 0.000 0.801 0.710
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .x1 0.835 0.118 7.064 0.000 0.835 0.614
## .x2 1.065 0.105 10.177 0.000 1.065 0.771
## .x3 0.633 0.129 4.899 0.000 0.633 0.496
## visual 0.524 0.130 4.021 0.000 1.000 1.000
##
## R-Square:
## Estimate
## x1 0.386
## x2 0.229
## x3 0.504library(lavaanPlot)
lavaanPlot(model = fit,
stand = TRUE, coefs = T)resid_cov = lavResiduals(fit)$cov
(resid_cov > 0.1) | (resid_cov < -0.1)## x1 x2 x3
## x1 FALSE FALSE FALSE
## x2 FALSE FALSE FALSE
## x3 FALSE FALSE FALSEТекстовые навыки
library(lavaan)
library(semTools)
model = '
textual =~ x4 + x5 + x6
'
fit = cfa(model, data)
summary(fit, fit.measures = TRUE, standardized = TRUE, rsquare = T)## lavaan 0.6-19 ended normally after 15 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 6
##
## Number of observations 301
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 0
##
## Model Test Baseline Model:
##
## Test statistic 497.430
## Degrees of freedom 3
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 1.000
## Tucker-Lewis Index (TLI) 1.000
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1181.065
## Loglikelihood unrestricted model (H1) -1181.065
##
## Akaike (AIC) 2374.130
## Bayesian (BIC) 2396.372
## Sample-size adjusted Bayesian (SABIC) 2377.344
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.000
## P-value H_0: RMSEA <= 0.050 NA
## P-value H_0: RMSEA >= 0.080 NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.000
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## textual =~
## x4 1.000 0.984 0.847
## x5 1.133 0.067 16.906 0.000 1.115 0.866
## x6 0.924 0.056 16.391 0.000 0.910 0.832
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .x4 0.382 0.049 7.805 0.000 0.382 0.283
## .x5 0.416 0.059 7.038 0.000 0.416 0.251
## .x6 0.369 0.044 8.367 0.000 0.369 0.308
## textual 0.969 0.112 8.640 0.000 1.000 1.000
##
## R-Square:
## Estimate
## x4 0.717
## x5 0.749
## x6 0.692library(lavaanPlot)
lavaanPlot(model = fit,
stand = TRUE, coefs = T)resid_cov = lavResiduals(fit)$cov
(resid_cov > 0.1) | (resid_cov < -0.1)## x4 x5 x6
## x4 FALSE FALSE FALSE
## x5 FALSE FALSE FALSE
## x6 FALSE FALSE FALSEНавыки скорости
Обратите внимание, тут тоже добавляем взаимосвязь в модель!(теперь она, судя по всему, over-identified)
library(lavaan)
library(semTools)
model = '
speed =~ x7 + x8 + x9
x7 ~~ x8
'
fit = cfa(model, data)
summary(fit, fit.measures = TRUE, standardized = TRUE, rsquare = T)## lavaan 0.6-19 ended normally after 20 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 7
##
## Number of observations 301
##
## Model Test User Model:
##
## Test statistic NA
## Degrees of freedom -1
## P-value (Unknown) NA
##
## Model Test Baseline Model:
##
## Test statistic NA
## Degrees of freedom NA
## P-value NA
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) NA
## Tucker-Lewis Index (TLI) NA
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1233.814
## Loglikelihood unrestricted model (H1) -1233.814
##
## Akaike (AIC) 2481.627
## Bayesian (BIC) 2507.577
## Sample-size adjusted Bayesian (SABIC) 2485.377
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000
## 90 Percent confidence interval - lower NA
## 90 Percent confidence interval - upper NA
## P-value H_0: RMSEA <= 0.050 NA
## P-value H_0: RMSEA >= 0.080 NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.000
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## speed =~
## x7 1.000 0.805 0.740
## x8 1.225 NA 0.986 0.975
## x9 0.577 NA 0.464 0.461
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .x7 ~~
## .x8 -0.258 NA -0.258 -1.566
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .x7 0.536 NA 0.536 0.453
## .x8 0.051 NA 0.051 0.049
## .x9 0.800 NA 0.800 0.788
## speed 0.647 NA 1.000 1.000
##
## R-Square:
## Estimate
## x7 0.547
## x8 0.951
## x9 0.212library(lavaanPlot)
lavaanPlot(model = fit,
stand = TRUE, coefs = T)resid_cov = lavResiduals(fit)$cov
(resid_cov > 0.1) | (resid_cov < -0.1)## x7 x8 x9
## x7 FALSE FALSE FALSE
## x8 FALSE FALSE FALSE
## x9 FALSE FALSE FALSEАнализ надежности
Возвращаемся к нашей большой финальной модели с тремя факторами. Оценим её надежность (omega должна быть больше 0.7 для каждого фактора)
library(semTools)
reliability(fit_final)## visual textual speed
## alpha 0.6261171 0.8827069 0.6884550
## omega 0.6262675 0.8852443 0.5586804
## omega2 0.6262675 0.8852443 0.5586804
## omega3 0.6152789 0.8851594 0.5581609
## avevar 0.3697524 0.7211745 0.4041205Извлечение факторных значений
Теперь попробуем извлечь факторные значения
Обратите внимание, что вновь можно выбрать интересующиё вас метод извлечения. Попробуем вновь воспользоваться регрессией.
data_cfa_scores = HolzingerSwineford1939
data_cfa_scores$visual = as.data.frame(lavPredict(fit_final, type = "lv", method = "regression"))$visual
data_cfa_scores$textual = as.data.frame(lavPredict(fit_final, type = "lv", method = "regression"))$textual
data_cfa_scores$speed = as.data.frame(lavPredict(fit_final, type = "lv", method = "regression"))$speedПовторим наш предыдущий анализ
hist(data_cfa_scores$textual)
ggqqplot(data_cfa_scores$textual)
Отличаются ли “текстовые” способности у представителей разного пола?
t.test(textual ~ sex, data_cfa_scores)##
## Welch Two Sample t-test
##
## data: textual by sex
## t = -1.1099, df = 298.85, p-value = 0.268
## alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
## 95 percent confidence interval:
## -0.33040797 0.09211722
## sample estimates:
## mean in group 1 mean in group 2
## -0.06135393 0.05779144Вербальные способности у представителей разного пола на наших данных не отличаются
Вывод
СFA и EFA показали схожие результаты. Из данных можно выделить 3 фактора:
Фактор 1 (визуальные навыки):
x1 (тест 1) — Тест зрительного восприятия Spearman VPT (оценка анализа геометрических паттернов)
x2 (тест 2) — Упрощённый тест пространственных отношений Brigham (сборка кубов по образцу)
x3 (тест 4) — Тест «Лозенги» Thorndike (идентификация перевёрнутых фигур среди похожих)
Фактор 2 (текстовые навыки):
x4 (тест 6) — Тест понимания текста (ответы на вопросы по прочитанному абзацу)
x5 (тест 7) — Тест завершения предложений (подбор слов по контексту)
x6 (тест 9) — Тест значения слов (поиск синонимов/антонимов)
Фактор 3 (навыки скорости):
x7 (тест 10) — Тест скоростного сложения (быстрые арифметические вычисления)
x8 (тест 12) — Тест счёта точек (подсчёт точек внутри фигур на скорость)
x9 (тест 13) — Тест дискриминации линий (разделение прямых и кривых линий за время)
При этом в обоих случаях модели получились не самыми лучшими. Скорее такие факторы в исследованиях лучше не использовать (хотя может и можно, учитывая возможные ошибки при оценке значений факторов)
Задание
Проведите эксплораторный и конфирматорный факторный анализы на данных шкалы Розенберга
Данные представлены ниже:
library(EFA.dimensions)
rosenberg = data_RSE
str(rosenberg)## 'data.frame': 300 obs. of 10 variables:
## $ Q1 : int 4 4 4 4 2 2 1 4 1 3 ...
## $ Q2 : int 4 3 4 4 2 3 3 4 2 4 ...
## $ Q3 : int 2 3 2 3 3 3 4 3 4 2 ...
## $ Q4 : int 3 3 3 4 2 2 2 3 2 3 ...
## $ Q5 : int 1 2 1 3 4 3 3 2 4 2 ...
## $ Q6 : int 3 2 3 3 2 2 2 2 2 3 ...
## $ Q7 : int 2 2 3 3 2 2 2 1 1 3 ...
## $ Q8 : int 2 3 2 2 3 3 4 3 3 4 ...
## $ Q9 : int 3 4 2 3 3 3 4 4 4 4 ...
## $ Q10: int 2 3 2 2 2 3 4 4 4 4 ...Узнать больше о шкале самооценки Розенберга можно по ссылке:
https://www.apa.org/obesity-guideline/rosenberg-self-esteem.pdf
Обратите внимание:
Хотя данные и закодированы как числа, они ординальные (вспомните об этом при выборе метода построения корреляций и в ряде других случаев)
Вам не известно, какая переменная в датасете соответствует какому вопросу в оригинальном опросе, который можно посмотреть по ссылке выше
Среди переменных есть вопросы, ответы на которые можно инвертировать (например, I feel I do not have much to be proud of)
Кроме того, такие вопросы (заданые в негативном ключе) могут иметь коррелирующие ошибки, что может отразиться на выдаваемом количестве факторов в EFA (подсказка: фактор всё-таки должен быть один). Также обратите внимание на такие вопросы при анализе индексов модицикации в CFA
Узнать чуть больше о данных можно попробовать в документации пакета EFA.dimensions:
?EFA.dimensions:: data_RSE