Разница между proc princomp в SAS и командой princomp в R?

Question

В настоящее время я пытаюсь получить эквивалентные результаты с помощью команды proc princomp в SAS и команды princomp() в R (в пакете stats). Результаты, которые я получаю, очень похожи, что наводит меня на мысль, что это не проблема с различными настройками параметров в двух командах. Однако выход также достаточно различен, так что оценки компонентов для каждой строки данных заметно отличаются. Они также поменялись знаком, но это, конечно, не имеет значения.

Конечной целью этого анализа является получение набора коэффициентов из PCA для оценки данных вне процедуры PCA (то есть формула, которая может быть применена к новым наборам данных для легкого получения оцененных данных).

Не публикуя все мои данные, я надеюсь, что кто-то может предоставить некоторую информацию о том, как эти две команды могут отличаться в своих вычислениях. Я не знаю достаточно о математике PCA, чтобы определить, является ли это концептуальной разницей в процессах или просто чем-то вроде разницы во внутреннем округлении. Для простоты я опубликую собственные векторы только для ПК1 и ПК2.

В САС:

proc princomp data=climate out=pc_out outstat=pc_outstat; 
var MAT MWMT MCMT logMAP logMSP CMI cmiJJA DD_5 NFFD; 
run;

возвращается

Eigenvectors
       Prin1  Prin2  Prin3  Prin4  Prin5  Prin6  Prin7  Prin8  Prin9 
MAT    0.372  0.257  -.035  -.033  -.106  0.270  -.036  0.216  -.811 
MWMT   0.381  0.077  0.160  -.261  0.627  0.137  -.054  0.497  0.302 
MCMT   0.341  0.324  -.229  0.046  -.544  0.421  0.045  0.059  0.493 
logMAP -.184  0.609  -.311  -.357  -.041  -.548  0.183  0.183  0.000 
logMSP -.205  0.506  0.747  -.137  -.040  0.159  -.156  -.266  0.033 
CMI    -.336  0.287  -.451  0.096  0.486  0.499  0.050  -.318  -.031 
cmiJJA -.365  0.179  0.112  0.688  -.019  0.012  0.015  0.588  0.018 
DD_5   0.379  0.142  0.173  0.368  0.183  -.173  0.725  -.282  0.007 
NFFD   0.363  0.242  -.136  0.402  0.158  -.351  -.637  -.264  0.052 

В R:

PCA.model <- princomp(climate[,c("MAT","MWMT","MCMT","logMAP","logMSP","CMI","cmiJJA","DD.5","NFFD")], scores=T, cor=T)
PCA.model$loadings

возвращается

Eigenvectors
       Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8 Comp.9
MAT    -0.372 -0.269         0.126        -0.250         0.270  0.789
MWMT   -0.387        -0.171         0.675                0.494 -0.325
MCMT   -0.339 -0.332  0.250  0.164 -0.500 -0.414               -0.510
logMAP  0.174 -0.604  0.309  0.252         0.619 -0.213  0.125       
logMSP  0.202 -0.501 -0.727  0.223        -0.162  0.175 -0.268       
CMI     0.334 -0.293  0.459 -0.222  0.471 -0.495        -0.271       
cmiJJA  0.365 -0.199 -0.174 -0.612 -0.247                0.590       
DD.5   -0.382 -0.143 -0.186 -0.421               -0.695 -0.360       
NFFD   -0.368 -0.227        -0.487         0.309  0.655 -0.205  

Как видите, значения похожи (знак обратный), но не идентичны. Различия имеют значение в оцененных данных, первая строка которых выглядит следующим образом:

     Prin1  Prin2  Prin3  Prin4  Prin5  Prin6  Prin7  Prin8  Prin9 
SAS  -1.95   1.68  -0.54   0.72  -1.07   0.10  -0.66  -0.02   0.05
R     1.61  -1.99   0.52  -0.42  -1.13  -0.16   0.79   0.12  -0.09

Если я использую GLM (в SAS) или lm() (в R) для вычисления коэффициентов по оцененным данным, я получаю очень похожие числа (обратный знак), за исключением перехвата. Вот так:

в САС:

proc glm order=data data=pc_out;
model Prin1 = MAT MWMT MCMT logMAP logMSP CMI cmiJJA DD_5 NFFD;
run;

в R:

scored <- cbind(PCA.model$scores, climate)
pca.lm <- lm(Comp.1~MAT+MWMT+MCMT+logMAP+logMSP+CMI+cmiJJA+DD.5+NFFD, data=scored)

возвращается

    Coefficients:
    (Int)  MAT    MWMT   MCMT   logMAP  logMSP  CMI     cmiJJA  DD.5     NFFD 
SAS  0.42   0.04   0.06   0.03  -0.65   -0.69   -0.003  -0.01    0.0002   0.004
R   -0.59  -0.04  -0.06  -0.03   0.62    0.68    0.004   0.02   -0.0002  -0.004

Поэтому может показаться, что пересечение модели меняет значение в оцененных данных. Любые мысли о том, почему это происходит (почему перехват отличается) будет принята с благодарностью.

Answer 1

Еще раз спасибо всем, кто прокомментировал. Смущает то, что различия, которые я обнаружил между процедурами SAS proc princomp и R princomp(), на самом деле были результатом ошибки данных, которую я сделал. Извините тем, кто нашел время, чтобы помочь ответить.

Но вместо того, чтобы оставить этот вопрос напрасным, я предложу статистически эквивалентные процедуры для SAS и R при выполнении анализа основных компонентов (PCA).

Следующие процедуры статистически эквивалентны, с данными с именем "mydata" и переменными с именем "Var1", "Var2" и "Var3".

В САС:

* Run the PCA on your data;
proc princomp data=mydata out=pc_out outstat=pc_outstat; 
var Var1 Var2 Var3; 
run;
* Use GLM on the individual components to obtain the coefficients to calculate the PCA scoring;
proc glm order=data data=pc_out;
model Prin1 = Var1 Var2 Var3;
run;

В R:

PCA.model <- princomp(mydata[,c("Var1","Var2","Var3")], scores=T, cor=T)
scored <- predict(PCA.model, mydata)
scored <- cbind(PCA.model$scores, mydata)
lm(Comp.1~Var1+Var2+Var3, data=scored)