Почему существуют различия в psych::principal между Varimax и Varimax?

В связанном вопросе я спросил, почему существуют различия между stats::varimax и GPArotation::Varimax, оба из которых вызываются psych::principal, в зависимости от параметра, установленного для rotate =.

Различия между этими двумя (см. Другой вопрос) объясняют некоторые, но не все отличия от psych::principal. Похоже, что эти различия каким-то образом усугубляются psych::principal. (У меня есть простая теория, почему, и я хотел бы получить ее подтверждение).

library(GPArotation)
library(psych)
data("Thurstone")

principal.unrotated <- principal(r = Thurstone, nfactors = 4, rotate = "none")  # find unrotated PCs first
loa <- unclass(principal.unrotated$loadings)

# just to compare that the rot.mat is correct
varimax.stats <- stats::varimax(x = loa, normalize = TRUE)
varimax.GPA <- GPArotation::Varimax(L = loa, normalize = TRUE)
# notice we're here NOT interested in the difference between stats and GPA, that's the other question

diff.from.rot.meth <- unclass(varimax.stats$loadings - varimax.GPA$loadings)  # very small differences, see this question: https://stackoverflow.com/questions/32350891/why-are-there-differences-between-gparotationvarimax-and-statsvarimax
mean(abs(diff.from.rot.meth))
#> [1] 8.036863e-05

principal.varimax.stats <- principal(r = Thurstone, nfactors = 4, rotate = "varimax")
principal.Varimax.GPA <- principal(r = Thurstone, nfactors = 4, rotate = "Varimax")

diff.from.princ <- principal.Varimax.GPA$rot.mat - principal.varimax.stats$rot.mat  # quite a substantial change, because Theta is NOT a rotmat, that makes sense
mean(abs(diff.from.princ))
#> [1] 0.021233

mean(abs(diff.from.rot.meth)) - mean(abs(diff.from.princ))  # principal has MUCH bigger differences
#> [1] -0.02115263

Это кажется слишком большим для артефакта с плавающей запятой или чего-то в этом роде.

Моя гипотеза заключается в том, что (дополнительная) разница связана с тем фактом, что GPArotation::Varimax по умолчанию имеет значение (Kaiser) normalize == FALSE, тогда как **stats::varimax по умолчанию имеет значение (Kaiser) normalize == TRUE, что не может быть установлен иначе, чем внутри `principal::psych``.

stats::varimax инструкция:

## varimax with normalize = TRUE is the default

GPArotation::Varimax / GPArotation::GPForth руководство:

Аргумент normalize дает указание на то, должна ли и как должна быть выполнена какая-либо нормализация перед поворотом, а затем отменена после поворота. Если normalize имеет значение FALSE (по умолчанию), нормализация не выполняется. Если normalize имеет значение TRUE, выполняется нормализация по Кайзеру. (Поэтому квадраты строк нормализованной суммы A составляют 1,0. Иногда это называют нормализацией Хорста.)

Кроме того, они psych::Kaiser руководство предупреждает:

Пакет GPArotation не нормализует (по умолчанию) и не работает fa. Затем, чтобы еще больше запутать, варимакс в статистике работает, а варимакс в GPArotation — нет.

Может ли кто-нибудь подтвердить, что разница на самом деле объясняется параметрами нормализации?


person maxheld    schedule 02.09.2015    source источник


Ответы (2)


Проблема разных нагрузок возникает из-за разной точности процедур (разумеется, поскольку в psych::principal параметр normalize не оценивается, все остальные процедуры должны использоваться с этим параметром, установленным в TRUE) . В то время как точность в stats::varimax и GPArotation::Varimax может быть настроена (параметр eps), она игнорируется в psych::principal и кажется неявно равной stats::varimax с eps=1e-5 .

Если мы увеличим точность в stats::varimax и GPArotations::Varimax до eps=1e-15, то получим тот же результат (по крайней мере, до 8 цифр), что и в параллельной реализации в моей программе MatMate, которую я проверил, чтобы быть очень точно равным к вычислениям SPSS тоже.

Отсутствующая обработка явной опции eps в psych::principal кажется ошибкой, и ее плохое неявное значение по умолчанию, безусловно, неудовлетворительно.

Интересно, что GPArotation::Varimax требует очень много вращений с eps=1e-15, см. вывод в конце; так что либо реализована другая внутренняя процедура, либо eps-параметр оценивается по-разному в решении, когда итерации можно остановить. Один пример, см. в конце этого ответа, может предложить такой эффект.

Ниже протокол сравнения; отображаются только первые две строки загрузок.

The first three computations with low accuracy (eps=1e-5) 
all procedures give equal or comparable results, 
except GPArotation, which is already near the exact answer
--------------------------------
>  principal(r = Thurstone, nfactors = 4, rotate = "varimax")$loadings    *1000000
>  stats::varimax(x = loa, normalize = TRUE,eps=1e-5)$loadings            *1000000
>  GPArotation::Varimax(L = loa, normalize = TRUE,eps=1e-5)$loadings      *1000000

The second three computations with (attempted) high accuracy (eps=1e-15) 
all procedures except psych::principal give equal results, and they
agree also with external crosscheck using MatMate
--------------------------------
>  principal(r = Thurstone, nfactors = 4, rotate = "varimax",eps=1e-15)$loadings*1000000
>  stats::varimax(x = loa, normalize = TRUE,eps=1e-15)$loadings*1000000
>  GPArotation::Varimax(L = loa, normalize = TRUE,eps=1e-10)$loadings*1000000

Попытки получить маленькие/по умолчанию точные результаты с большими eps или без них

# ===== Numerical documentation (only first two rows are displayed)================
> principal(r = Thurstone, nfactors = 4, rotate = "varimax")$loadings*1000000
                RC1       RC2       RC3       RC4      
Sentences       871655.72 216638.46 198427.07 175202.57
Vocabulary      855609.28 294166.99 153181.45 180525.99


> stats::varimax(x = loa, normalize = TRUE,eps=1e-5)$loadings         *1000000
                PC1       PC2       PC3       PC4      
Sentences       871655.72 216638.46 198427.07 175202.57
Vocabulary      855609.28 294166.99 153181.45 180525.99

> GPArotation::Varimax(L = loa, normalize = TRUE,eps=1e-5)$loadings   *1000000
                     PC1      PC2      PC3       PC4
Sentences       871717.3 216618.7 198176.3 175204.47
Vocabulary      855663.1 294146.3 152930.7 180517.21
# =============================================================================

Теперь попытки получить более точные результаты с меньшими eps

> principal(r = Thurstone, nfactors = 4, rotate = "varimax",eps=1e-15)$loadings  *1000000
                RC1       RC2       RC3       RC4      
Sentences       871655.72 216638.46 198427.07 175202.57
Vocabulary      855609.28 294166.99 153181.45 180525.99

> stats::varimax(x = loa, normalize = TRUE,eps=1e-15)$loadings                   *1000000
                PC1       PC2       PC3       PC4      
Sentences       871716.83 216619.69 198174.31 175207.86
Vocabulary      855662.58 294147.47 152928.77 180519.37

> GPArotation::Varimax(L = loa, normalize = TRUE,eps=1e-10)$loadings             *1000000
                     PC1      PC2       PC3       PC4
Sentences       871716.8 216619.7 198174.31 175207.86
Vocabulary      855662.6 294147.5 152928.77 180519.37

Warnmeldung:
In GPForth(L, Tmat = Tmat, method = "varimax", normalize = normalize,  :
  convergence not obtained in GPForth. 1000 iterations used.

# Result by MatMate: --------------------------------------------------------
 lad = cholesky(Thurstone) 
 pc = rot(lad,"pca")
 pc4 = pc[*,1..4]                           // arrive at the first four pc's
     t = gettrans( normzl(pc4),"varimax")   // get rotation-matrix for row-normalized pc's
 vmx = pc4 * t                              // rotate pc4 by rotation-matrix 
 display = vmx     * 1000000
                     PC1      PC2       PC3       PC4
Sentences       871716.83    216619.68   198174.31   175207.87
Vocabulary      855662.58    294147.46   152928.77   180519.37


# ===============================================================================> 

Гораздо лучшего совпадения результатов по stats::varimax и GPArotation::Varimax можно было бы добиться, установив eps на 1e-12 в stat:: и на 1e-6 в GPArotation, где одно значение равно квадрату другого. Мы получаем тогда

> GPArotation::Varimax(L = loa, normalize = TRUE,eps=1e-6)$loadings*1000000
                     PC1      PC2       PC3       PC4
Sentences       871716.8 216619.8 198174.49 175207.63
Vocabulary      855662.5 294147.6 152928.94 180519.21

> stats::varimax(x = loa, normalize = TRUE,eps=1e-12)$loadings*1000000
                PC1       PC2       PC3       PC4      
Sentences       871716.80 216619.74 198174.40 175207.85
Vocabulary      855662.55 294147.52 152928.86 180519.36
person Gottfried Helms    schedule 19.05.2016

Казалось бы, это подтверждает, что применение psych::kaiser (которое, я думаю, было создано для этой цели) сводит различия обратно к исходным различиям между stats::varimax и GPArotation::Varimax:

principal.Varimax.GPA.kaiser <- kaiser(f = principal.unrotated, rotate = "Varimax")
diff.statsvari.gpavar.bothkaiser <- unclass(principal.Varimax.GPA.kaiser$loadings - principal.varimax.stats$loadings)
mean(abs(diff.statsvari.gpavar.bothkaiser))
#> [1] 8.036863e-05

Это почти тот же результат, поэтому я думаю, что гипотеза подтверждена.

Большая разница, связанная с psych::principal, связана с разными значениями по умолчанию для normalize.


Обновить

Различия также (опять же) намного меньше для соответствующих матриц вращения (или что-то еще Th):

principal.Varimax.GPA.kaiser$Th - principal.varimax.stats$rot.mat  # those differences are very small now, too
#>               [,1]         [,2]          [,3]          [,4]
#> [1,]  1.380279e-04 1.380042e-05 -2.214319e-04 -2.279170e-06
#> [2,]  9.631517e-05 2.391296e-05  1.531723e-04 -3.371868e-05
#> [3,]  1.758299e-04 7.917460e-05  6.788867e-05  1.099072e-04
#> [4,] -9.548010e-05 6.500162e-05 -1.679753e-05 -5.213475e-05
person maxheld    schedule 02.09.2015