Макрос расчета канонических корреляций

 

Загрузить как .sps Показать простым текстом

 

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
* Макрос расчета канонических корреляций.

* Поставляется с SPSS.
* Пример вызова макроса - в конце файла.

preserve.                                                                       
set printback=off.                                                              
define cancorr (set1  =!charend('/')                                            
               /set2  =!charend('/')                                            
               /debug =!charend('/') !DEFAULT ('N')                             
               /KEEPSC=!charend('/') !DEFAULT ('Y') 
               /PRCOR =!charend('/') !DEFAULT (25 ) ).                                                                                                        
preserve.
!IF ( !DEBUG !EQ 'N') !THEN
set printback=off mprint off.                                                   
!ELSE
set printback on mprint on.
!IFEND .

*---------------------------------------------------------------------------.
* Сохранение исходного файла для возможного последующего восстановления.
*---------------------------------------------------------------------------.
!IF (!DEBUG !EQ 'N') !THEN 
SET RESULTS ON.
DO IF $CASENUM=1.
PRINT / "ВНИМАНИЕ: ВСЯ ВЫДАЧА, ВКЛЮЧАЯ СООБЩЕНИЯ ОБ ОШИБКАХ, ВРЕМЕННО ПРИОСТАНОВЛЕНА." 
      / "ЕСЛИ ВАМ КАЖЕТСЯ, ЧТО ПРОГРАММА РАБОТАЕТ НЕ ПРАВИЛЬНО, ПОВТОРНО ЗАПУСТИТЕ" 
      / "ЭТОТ МАКРОС С ДОПОЛНИТЕЛЬНЫМ АРГУМЕНТОМ  /DEBUG='Y'."  
      / "ПЕРЕД ПОВТОРНЫМ ЗАПУСКОМ НЕОБХОДИМО ЗАНОВО ЗАГРУЗИТЬ ФАЙЛ ДАННЫХ."
      / "УКАЗАННЫЙ ПАРАМЕТР ДАСТ ИНФОРМАЦИЮ ДЛЯ ДИАГНОСТИРОВАНИЯ ВОЗНИКАЮЩИХ ПРОБЛЕМ".
END IF.
!IFEND .
save outfile='cc__tmp1.sav'.

*---------------------------------------------------------------------------.
* Построение корреляционной матрицы и передача ее в команду MATRIX.
*---------------------------------------------------------------------------.

* ПО УМОЛЧАНИЮ:  ВЫДАЧА РЕЗУЛЬТАТОВ И ОШИБОК ОТКЛЮЧЕНА ДЛЯ ПРЕДОТВРАЩЕНИЯ ВЫВОДА МАТРИЦЫ КОРРЕЛЯЦИЙ В ОКНО РЕЗУЛЬТАТОВ *.
!IF (!DEBUG='N') !THEN
set results off errors off.
!IFEND 
corr variables=!set1 !set2 /missing=listwise/matrix out(*).                                                                              
set errors on results listing.                                                                      

*---------------------------------------------------------------------------.
* Вычисление канонич. корреляций и основных величин, необходимых для анализа.
*---------------------------------------------------------------------------.
* Команда SET помещена для изменения установки по умолчанию на число циклов MXLOOPS=40 * .                

SET MXLOOPS=199 MITERATE 199.
matrix.
get r /variables=!set1/file=*.
compute p1=ncol(r).
get r /file=* /names=varname/variables=!set1 !set2.
compute p2=ncol(r)-p1.
compute nx1=varname(1:p1).
compute nv=p1+p2.
compute nx2=varname((p1+1):nv).
compute rr=r(4:(nv+3),1:nv).
compute ns=r(3,1).
compute r11=rr(1:p1,1:p1).
compute r22=rr((p1+1):nv,(p1+1):nv).
compute r12=rr(1:p1,(p1+1):nv).
compute d1=r(2,1:p1).
compute d2=r(2,(p1+1):nv).
compute d1=mdiag(d1).
compute d2=mdiag(d2).
compute s1=d1*r11*d1.
compute s12=d1*r12*d2.
compute s2=d2*r22*d2.
compute d1=inv(d1).
compute d2=inv(d2).

*---------------------------------------------------------------------------.
* Вывод всех корреляций.
*---------------------------------------------------------------------------.

do if (p1 <= !PRCOR and p2 <= !PRCOR ).
print r11 /format "f7.4"/title 'Корреляции для набора №1'                         
          /space 2 /rnames=nx1 /cnames=nx1.
print r22 /format "f7.4"/title 'Корреляции для набора №2'
          /space 2 /rnames=nx2 /cnames=nx2.
print r12 /format "f7.4"/title 'Корреляции между переменными наборов №1 и №2 '
         /space 2 /rnames=nx1 /cnames=nx2.
end if.  
                                                                             
*---------------------------------------------------------------------------.
* R1 и r2 - разложения Холецкого для матриц корреляций r11 и r22, соответственно.
*---------------------------------------------------------------------------.

compute r1=chol(r11).
compute r2=chol(r22).

*---------------------------------------------------------------------------.
* R1_inv и r2_inv - матрицы, обратные к r1 и r2, соответственно.
*---------------------------------------------------------------------------.

compute r1_inv=inv(r1).
compute r2_inv=inv(r2).

*---------------------------------------------------------------------------.
* Вычисление матрицы "омега".
*---------------------------------------------------------------------------.

do if (p1 le p2).
compute omega=t(r1_inv)*r12*r2_inv.
else.
compute omega=t(r2_inv)*t(r12)*r1_inv.
end if.

*---------------------------------------------------------------------------.
* Функция SVD производит разложение по сингулярным числам матрицы "омега".
*---------------------------------------------------------------------------.

call svd(omega,u,lambda,v).

*---------------------------------------------------------------------------.
* Создание набора имен, для пометки результатов при выводе .
*---------------------------------------------------------------------------.

!LET !@=!NULL !LET !@1=!NULL !LET !@2=!NULL                                     
  !DO !N= 1 !TO 199                                                             
      !LET !@=!CONCAT(!@,!QUOTE(!N),",")                                        
      !LET !@1=!CONCAT(!@1,!QUOTE(!CONCAT('CV1-',!N)),",")                      
      !LET !@2=!CONCAT(!@2,!QUOTE(!CONCAT('CV2-',!N)),",")                      
  !DOEND                                                                        
  !LET !@=!CONCAT(!@,!QUOTE(@@))                                                
  !LET !@1=!CONCAT(!@1,!QUOTE(@@))                                              
  !LET !@2=!CONCAT(!@2,!QUOTE(@@)).                                                                               
Compute num={!@}.

*---------------------------------------------------------------------------.
* Матрица "лямбда" содержит канонические корреляции. Вывод их на печать.
*---------------------------------------------------------------------------.

print diag(lambda)/format "f8.3"/title 'Канонические корреляции'                 
          /space 2/rnames=num.

compute dlam=diag(lambda).

*---------------------------------------------------------------------------.
* Вычисление собственных значений и проверки на наличие остальных канонических корреляций (вторых, третьих пар и т.д.).
*---------------------------------------------------------------------------.

compute eign=(1 &/ (1-dlam &**2)) - 1.
compute wlam=1 &/ (1+eign).
compute n=nrow(wlam).
compute wilk=wlam.
compute df=wlam.
compute sig=wlam.
compute bart2=wlam.
compute tem=1.
loop  #l=1 to n.
+  compute tem=tem*wlam(n-#l+1).
+  compute df(n-#l+1)=(p1-n+#l)*(p2-n+#l).
+  compute dof=df(n-#l+1).
+  compute bart2(n-#l+1)=-(ns-0.5*(p1+p2+3))*ln(tem).
+  compute chi=bart2(n-#l+1).
+  compute sig(n-#l+1)=1-chicdf(chi,dof).
+  compute wilk(n-#l+1)=tem.
end loop.
compute test={wilk,bart2,df,sig}.
print test /format "f8.3"/title 'Проверка, что остальные канонические корреляции равны нулю:'
  /space 2/rnames=num
  /cnames={"Стат. Уилка","Хи-квадрат","  Ст. св.  ","  Знач."}.

*---------------------------------------------------------------------------.
* Вычисление и вывод стандартизированных канонических коэффициентов для набора №1.
*---------------------------------------------------------------------------.

do if (p1 le p2).
compute a=r1_inv*u.
else.
compute a=r1_inv*v.
end if.
do if (p2 lt p1).
compute a=a(:,1:p2).
end if.
print a/format "f8.3"/title 'Стандарт. канонич. коэфф. для набора №1'     
      /space 2/rnames=nx1/cnames=num.

*---------------------------------------------------------------------------.
* Вычисление и вывод нестандартизированных канонических коэффициентов для набора №1.
*---------------------------------------------------------------------------.

compute a1=d1*a.
print a1/format "f8.3"/title 'Нестандарт. канонич. коэфф. для набора №1'             
      /space 2/rnames=nx1/cnames=num.

*---------------------------------------------------------------------------.
* Вычисление и вывод стандартизированных канонических коэффициентов для набора №2.
*---------------------------------------------------------------------------.

do if (p1 le p2).
compute b=r2_inv*v.
else.
compute b=r2_inv*u.
end if.
do if (p1 le p2).
compute b=b(:,1:p1).
end if.
print b /format "f8.3"/title 'Стандарт. канонич. коэфф. для набора №1'    
      /space 2/rnames=nx2/cnames=num.

*---------------------------------------------------------------------------.
* Вычисление и вывод нестандартизированных канонических коэффициентов для набора №2.
*---------------------------------------------------------------------------.

compute b1=d2*b.
print b1/format "f8.3"/title 'Нестандарт. канонич. коэфф. для набора №1'             
      /space 2/rnames=nx2/cnames=num.

*---------------------------------------------------------------------------.
* Вычисление и вывод на печать канонических нагрузок для набора №1.
*---------------------------------------------------------------------------.

compute tem=d1*s1*a1.
print tem /format "f8.3"/title 'Канонические нагрузки для набора №1'
   /space 2/rnames=nx1/cnames=num.

*---------------------------------------------------------------------------.
* Вычисление индекса избыточности (redundancy index) как пропорции дисперсии
* в наборе №1, объясненной его собственными каноническими величинами.
*---------------------------------------------------------------------------.

compute f1=cssq(tem)/p1.
compute f1=t(f1).

*---------------------------------------------------------------------------.
* Вычисление и печать перекрестных нагрузок для набора №1.
*---------------------------------------------------------------------------.

compute tem=d1*s12*b1.
print tem /format "f8.3"/title 'Перекрестные нагрузки для набора №1'                       
   /space 2/rnames=nx1/cnames=num.
                                                                               
*---------------------------------------------------------------------------.
* Вычисление индекса избыточности как пропорции дисперсии в наборе №1, 
* объясненной каноническими величинами набора №2.
*---------------------------------------------------------------------------.                                                                               

compute cs3=cssq(tem)/p1.
compute cs3=t(cs3).

*---------------------------------------------------------------------------.
* Вычисление и печать канонических нагрузок для набора №2.
*---------------------------------------------------------------------------.

compute tem=d2*s2*b1.
print tem /format "f8.3"/title 'Канонические нагрузки для набора №2'                   
   /space 2/rnames=nx2/cnames=num.

*---------------------------------------------------------------------------.
* Вычисление индекса избыточности как пропорции дисперсии в наборе №2, 
* объясненной его собственными каноническими величинами.
*---------------------------------------------------------------------------.

compute f2=cssq(tem)/p2.
compute f2=t(f2).

*---------------------------------------------------------------------------.
* Вычисление и печать перекрестных нагрузок для набора №2.
*---------------------------------------------------------------------------.

compute tem=d2*t(s12)*a1.
print tem /format "f8.3"/title 'Перекрестные нагрузки для набора №2'                       
   /space 2/rnames=nx2/cnames=num.

*---------------------------------------------------------------------------.
* Вычисление индекса избыточности как пропорции дисперсии в наборе №2, 
* объясненной каноническими величинами из набора №1.
*---------------------------------------------------------------------------.

compute cs4=cssq(tem)/p2.
compute cs4=t(cs4).

*---------------------------------------------------------------------------.
* Печать результатов анализа избыточности.
*---------------------------------------------------------------------------.

compute c1={!@1}.
compute c2={!@2}.

print /title '            Анализ избыточности:' /space 2.
print f1/format "f15.3"
  /title 'Доля дисперсии набора №1, объясненная его собственными канонич. велич.'
    /space 2/rnames=c1/cnames= {"Доля дисп."}.
print cs3/format "f15.3"
  /title 'Доля дисперсии набора №2, объясненная канонич. велич. набора №2'
    /space 2/rnames=c2/cnames= {"Доля дисп."}.
print f2/format "f15.3"
  /title 'Доля дисперсии набора №2, объясненная его собственными канонич. велич.'
    /space 2/rnames=c2/cnames= {"Доля дисп."}.
print cs4/format "f15.3"
  /title 'Доля дисперсии набора №2, объясненная канонич. велич. набора №1'
    /space 2/rnames=c1/cnames= {"Доля дисп."}.

*---------------------------------------------------------------------------.
* Создание файлов для использования при вычислении значений канонических переменных (canonical scores).
*---------------------------------------------------------------------------.

SAVE {P1,P2} / OUTFILE 'CC__SIZE.SAV'.
SAVE {T(A1),T(B1)} / OUTFILE 'CC__AB.SAV' .
END MATRIX.
*---------------------------------------------------------------------------.
* Создание файла с именами переменных и переменной номера набора            .
*---------------------------------------------------------------------------.

SET MESSAGES OFF RESULTS OFF.
SELECT IF $CASENUM=1.
DO REPEAT V=!SET1.
COMPUTE V=1.
END REPEAT.
DO REPEAT V=!SET2.
COMPUTE V=2.
END REPEAT.
STRING VARNAME (A8).
COMPUTE VARNAME='SET_NUM'.
FLIP VARIABLES !SET1 !SET2 / NEWNAMES=VARNAME .
COMPUTE VARSEQ=1.
SPLIT FILE BY SET_NUM.
CREATE VARSEQ=CSUM(VARSEQ).
SAVE OUTFILE 'CC_NAMES.SAV'.
GET FILE 'CC__SIZE.SAV' .

*---------------------------------------------------------------------------.
* Подготовка необходимой информации для записи инструкций COMPUTE для скоринга.
*---------------------------------------------------------------------------.

WRITE OUTFILE 'CC__AB.INC'
   /'STRING @NMA001 TO @NMA',COL1 (N3), ' (A8)'
   /'       @NMB001 TO @NMB',COL2 (N3), ' (A8)'
   /'VECTOR @NMA= @NMA001 TO @NMA',COL1 (N3)
   /'       @NMB= @NMB001 TO @NMB',COL2 (N3)
   /'COMPUTE N_A='COL1 (N3)
   /'COMPUTE N_B='COL2 (N3)
   /'IF (SET_NUM=1) @NMA(VARSEQ)=CASE_LBL'
   /'IF (SET_NUM=2) @NMB(VARSEQ)=CASE_LBL'
   /'COMPUTE @=1'
   /'AGGREGATE OUTFILE "CC__SPRD.SAV" / BREAK @'
   / ' / N_A=MAX(N_A) / N_B=MAX(N_B)'
   / ' / @NMA001 TO @NMA',COL1 (N3) '=MAX (@NMA001 TO @NMA',COL1 (N3),')'
   / ' / @NMB001 TO @NMB',COL2 (N3) '=MAX (@NMB001 TO @NMB',COL2 (N3),')'
   / 'GET FILE "CC__AB.SAV"'
   / 'COMPUTE @=1'
   / 'MATCH FILES FILE * / TABLE "CC__SPRD.SAV"/BY @'
   / 'VECTOR @NMA= @NMA001 TO @NMA',COL1 (N3)
   / '       @NMB= @NMB001 TO @NMB',COL2 (N3)
   / '       COEF= COL1 TO @'.
EXECUTE.

GET FILE 'CC_NAMES.SAV'.
INCLUDE FILE 'CC__AB.INC'.
SET PRINTBACK OFF.

*---------------------------------------------------------------------------.
* Запись инструкций COMPUTE для скоринга.
*---------------------------------------------------------------------------.

STRING @SCNM@ (A8).
STRING @OLDNM@ (A8).
COMPUTE @SCNM@=CONCAT('S1_CV',STRING($CASENUM,N3)).
WRITE OUTFILE 'CC__.INC' /'COMPUTE ',@SCNM@ ,'= 0'.
LOOP CC@@@ = 1 TO N_A.
COMPUTE @OLDNM@=@NMA(CC@@@).
COMPUTE @COEF@ =COEF(CC@@@).
WRITE OUTFILE 'CC__.INC' / ' +',@COEF@ (comma18.16),' * ',@OLDNM@ .
END LOOP.
COMPUTE @SCNM@=CONCAT('S2_CV',STRING($CASENUM,N3)).
WRITE OUTFILE 'CC__.INC' /'COMPUTE ',@SCNM@ ,'= 0'.
LOOP CC@@@=1 TO N_B.
COMPUTE @OLDNM@=@NMB(CC@@@).
COMPUTE @COEF@ =COEF(CC@@@+N_A).
WRITE OUTFILE 'CC__.INC' / ' +',@COEF@ (comma18.16),' * ',@OLDNM@ .
END LOOP.
EXECUTE.

*---------------------------------------------------------------------------.
* Загрузка исходного файла данных и запуск скоринговой программы.
*---------------------------------------------------------------------------.

GET FILE 'cc__tmp1.sav'.
INCLUDE FILE 'CC__.INC' .
ERASE FILE 'CC__SIZE.SAV' .
ERASE FILE 'CC__AB.INC'.
ERASE FILE 'CC_NAMES.SAV'.
ERASE FILE 'CC__AB.SAV'.
ERASE FILE "CC__SPRD.SAV" 
!IF (!KEEPSC ='N') !THEN
ERASE FILE 'CC__.INC'.
!ELSE

DO IF ($CASENUM=1).
SET RESULTS ON.
PRINT /'Оценки значений канонических переменных записаны в активный файл данных.'
      /'Был также создан файл со скоринговой программой для SPSS.'
      /'Чтобы использовать эту программу, откройте другой файл данных'
      /'С ТЕМИ ЖЕ переменными, что были использованы в данном анализе. '
      /'Затем примените команду INCLUDE для запуска скоринговой программы.'
      /'Например :' /
      /'GET FILE какой-тодругойфайлданных'
      /'INCLUDE FILE "CC__.INC".'
      /'EXECUTE.'.
END IF.
EXECUTE.

!IFEND.
RESTORE.
!enddefine.
RESTORE.

* Примеч. перев:.
* Пример вызова макроса (для вызова - убрать "*" перед командой cancorr):.
*cancorr set1=x1 x2 x3 /set2=y1 y2 y3 y4.
* Изучение канонических корреляций для наборов переменных x1 - x3 и y1 - y4.
* Возможные дополнительные параметры:
* debug - запрос дополнительной информации о работе макроса ("N" или "Y").
* PRCOR - максимальная размерность выводимой на печать матрицы корреляций.
* KEEPSC - оценка значений канонических переменных ("N" или "Y").
Related pages

...

Navigate from here