Практическая работа №5

6. Непараметрические коэффициенты корреляции. Значимость частных коэффициентов регрессии

7. Факторный анализ

Глушков Егор Александрович, гр. 20.М04-мм
Вариант 4

---

6. Непараметрические коэффициенты корреляции. Значимость частных коэффициентов регрессии

Данные addicts.xls. Переменные (вариант 4):

• \(X_1\) – sstati – оценка уровня тревожности
• \(X_2\) – asi1_med – индекс зависимости тяжести медицинского статуса
• \(X_3\) – asi2_emp – индекс зависимости тяжести социального статуса (занятость, заработок)
• \(X_4\) – asid3_dyr – длительность героиновой зависимости
  
• \(Y\) – asi7_psy – индекс зависимости тяжести психиатрического статуса

1. Построить корреляционную матрицу по метрическим переменным \(X_1, ... , X_4\) и выделить наиболее значимую корреляцию. Для этого случая построить двумерную диаграмму.

2. Применить модель множественной регресссии для зависимой переменной \(Y\) и для независимых переменных \(X_1, ... , X_4\). Определить значимость прогноза и указать наиболее значимые переменные.

3. Вычислить частный коэффициент корреляции \(\rho_{YX_1 \cdot X_2 X_3 X_4}\).

---

Корреляционная матрица, значимые корреляции

library(readxl)
addicts <- read_excel("addicts.xlsx")
# View(addicts)

Исследуем переменные на наличие пропусков. Выделим нужные переменные

data <- na.omit(addicts[ , c("sstati", "asi1_med", "asi2_emp", "asid3_dyr", "asi7_psy")])
summary(data)

     sstati         asi1_med         asi2_emp        asid3_dyr         asi7_psy   
 Min.   :23.00   Min.   :0.0000   Min.   :0.0000   Min.   : 0.500   Min.   :0.00  
 1st Qu.:43.00   1st Qu.:0.0000   1st Qu.:0.6900   1st Qu.: 2.000   1st Qu.:0.12  
 Median :48.00   Median :0.1900   Median :0.9100   Median : 3.000   Median :0.24  
 Mean   :48.52   Mean   :0.2447   Mean   :0.7834   Mean   : 3.605   Mean   :0.26  
 3rd Qu.:54.00   3rd Qu.:0.4200   3rd Qu.:1.0000   3rd Qu.: 4.000   3rd Qu.:0.36  
 Max.   :72.00   Max.   :1.0000   Max.   :1.0000   Max.   :24.000   Max.   :0.87  

Проверяем данные на нормальность

lshap <- lapply(data, shapiro.test)
lres <- sapply(lshap, `[`, c("statistic","p.value")); lres

          sstati    asi1_med     asi2_emp     asid3_dyr   asi7_psy    
statistic 0.9929767 0.8394149    0.7918929    0.7567038   0.9517803   
p.value   0.2184984 2.858381e-16 1.639616e-18 6.13559e-20 6.385948e-08

Лишь p.value, соответствующий переменной sstati, больше уровня значимости 0.05, что говорит о том, что распределение этих данных незначимо отличается от нормального. Другими словами, можно предположить их нормальность. Для других переменных гипотеза о нормальности отвергается, вследствие чего следует предпочесть непараметрические коэффициенты корреляции (Спирмена, Кендалла), чем более привычный коэффициент корреляции Пирсона.

cor(data, method="pearson")

               sstati    asi1_med     asi2_emp  asid3_dyr  asi7_psy
sstati    1.000000000  0.21570052  0.004541223 0.12211473 0.1947717
asi1_med  0.215700516  1.00000000 -0.031071218 0.09333766 0.1346685
asi2_emp  0.004541223 -0.03107122  1.000000000 0.11027514 0.1899920
asid3_dyr 0.122114730  0.09333766  0.110275143 1.00000000 0.1402910
asi7_psy  0.194771680  0.13466846  0.189992031 0.14029098 1.0000000

ms <- cor(data, method="spearman"); ms

               sstati     asi1_med     asi2_emp  asid3_dyr   asi7_psy
sstati    1.000000000  0.160604387  0.009650113 0.09534054 0.21629588
asi1_med  0.160604387  1.000000000 -0.009739093 0.05920707 0.06179187
asi2_emp  0.009650113 -0.009739093  1.000000000 0.09440814 0.17970881
asid3_dyr 0.095340538  0.059207069  0.094408140 1.00000000 0.12806424
asi7_psy  0.216295882  0.061791868  0.179708808 0.12806424 1.00000000

cor(data, method="kendall")

               sstati     asi1_med     asi2_emp  asid3_dyr   asi7_psy
sstati    1.000000000  0.116505322  0.006477093 0.06777330 0.15059703
asi1_med  0.116505322  1.000000000 -0.005463709 0.04544845 0.04726184
asi2_emp  0.006477093 -0.005463709  1.000000000 0.07372656 0.13120600
asid3_dyr 0.067773302  0.045448452  0.073726560 1.00000000 0.09631333
asi7_psy  0.150597034  0.047261842  0.131206003 0.09631333 1.00000000

show_cor_test <- function(x1, x2) {
  sp <- cor.test(x1, x2, method="spearman", exact=FALSE)
  k <- cor.test(x1, x2, method="kendall", exact=FALSE)
  df <- data.frame(cbind(c(sp$estimate, sp$p.value), c(k$estimate, k$p.value)), 
                   row.names=c("corr coeff (rho | tau resp)", "p.value"))
  colnames(df)<-c('spearman', 'kendall')
  df
}

Выделим наиболее значимые корреляции: уровень тревожности и психиатрический статус

show_cor_test(data$sstati, data$asi7_psy)

Как видно, \(\rho\) и \(\tau\) не так высоки (0.216 и 0.15), однако с учетом соответствующих p.value можно принять гипотезы о том, что данные коэффициенты корреляции Спирмена и Кендалла ненулевые.

Аналогично с уровнем тревожности и медицинским статусом

show_cor_test(data$sstati, data$asi1_med)

Коэффициенты корреляции не высоки, однако p.value говорят о том, что коэффициенты не равны нулю

Такая же ситуация с психиатрическим и социальным статусами

show_cor_test(data$asi7_psy, data$asi2_emp)

Построим heatmap с соответствующими коэффициентами корреляции Спирмена:

library(GGally)
ggcorr(data, method=c("pairwise", "spearman"), label=TRUE)

Двумерная диаграмма для уровня тревожности и психиатрического статуса:

str <- paste("r", round(cor(data$sstati, data$asi7_psy, method='spearman'), 3), sep="=")
qplot(sstati, asi7_psy, data = data, geom = c("point", "smooth"), method = "lm") + labs(title=str) + theme(plot.title = element_text(hjust = 0.5))

Двумерная диаграмма для уровня тревожности и медицинского статуса:

str <- paste("r", round(cor(data$sstati, data$asi1_med, method='spearman'), 3), sep="=")
qplot(sstati, asi1_med, data = data, geom = c("point", "smooth"), method = "lm") + labs(title=str) + theme(plot.title = element_text(hjust = 0.5))

Двумерная диаграмма для социального и психиатрического статусов:

str <- paste("r", round(cor(data$asi2_emp, data$asi7_psy, method='spearman'), 3), sep="=")
qplot(asi2_emp, asi7_psy, data = data, geom = c("point", "smooth"), method = "lm") + labs(title=str) + theme(plot.title = element_text(hjust = 0.5))

Можно заметить, что имеющаяся линейная связь между этими парами переменных весьма слабая (тем не менее, значима)

Модель множественной регресссии

Применим модель множественной регрессии с asi7_psy как зависимой переменной

mreg <- lm(asi7_psy ~ sstati + asi1_med + asi2_emp + asid3_dyr, data=data)
summary(mreg)

Call:
lm(formula = asi7_psy ~ sstati + asi1_med + asi2_emp + asid3_dyr, 
    data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.36262 -0.13077 -0.02225  0.10122  0.62899 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept) -0.043488   0.067113  -0.648  0.51754   
sstati       0.003371   0.001239   2.721  0.00693 **
asi1_med     0.066214   0.040571   1.632  0.10383   
asi2_emp     0.124408   0.039730   3.131  0.00193 **
asid3_dyr    0.007287   0.004683   1.556  0.12085   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1764 on 272 degrees of freedom
Multiple R-squared:  0.09201,   Adjusted R-squared:  0.07866 
F-statistic: 6.891 on 4 and 272 DF,  p-value: 2.689e-05

Согласно \(p.value < 0.05\) (а именно согласно значимости прогноза 2.689e-05 < 0.05) отвергается гипотеза о том, что данная модель не зависит от предикторов. Значимыми (ненулевыми) коэффициентами с уровнем 0.05 являются значения при переменных sstati и asi2_emp, чего и следовало ожидать из проведенного выше корреляционного анализа (именно коэффициенты корреляции этих пар переменных были наибольними).

Проверим некоррелированность остатков [практически равны 0]

c(cor(mreg$residuals, data$sstati), cor(mreg$residuals, data$asi1_med ), cor(mreg$residuals, data$asi2_emp ), cor(mreg$residuals, data$asid3_dyr))

[1] -1.888976e-16  9.887577e-20 -1.548208e-16  2.494561e-17

Частный коэффициент корреляции \(\rho_{YX_1 \cdot X_2 X_3 X_4}\)

library(matlib)
ms

               sstati     asi1_med     asi2_emp  asid3_dyr   asi7_psy
sstati    1.000000000  0.160604387  0.009650113 0.09534054 0.21629588
asi1_med  0.160604387  1.000000000 -0.009739093 0.05920707 0.06179187
asi2_emp  0.009650113 -0.009739093  1.000000000 0.09440814 0.17970881
asid3_dyr 0.095340538  0.059207069  0.094408140 1.00000000 0.12806424
asi7_psy  0.216295882  0.061791868  0.179708808 0.12806424 1.00000000

\[ \rho_{YX_1 \cdot X_2 X_3 X_4} = - \frac{\Lambda_{Y X_1}}{\sqrt{\Lambda_{Y Y} \Lambda_{1 1}}}\] где \(\Lambda_{ij}\) – алгебраическое дополнение элемента \(\lambda_{ij}\) корреляционной матрицы \(\Lambda\)

rho <- - cofactor(ms, 5, 1) / sqrt(cofactor(ms, 5, 5) * cofactor(ms, 1, 1))
c(r.yx1=ms[1, 5], rho=rho)

    r.yx1       rho 
0.2162959 0.2034030 

\(r_{Y X_1} > \rho_{YX_1 \cdot X_2 X_3 X_4}\), то факторы \(X_2, X_3, X4\) [немного] искажают взаимосвязь между \(Y\) и \(X_1\) в сторону её увеличения

---

7. Факторный анализ

Задание.
Построить матрицу факторных нагрузок. Проинтерпретировать главные компоненты, определить вклад первых двух в общую дисперсию и построить двумерную диаграмму первых двух факторов.

---

Вклад главных компонент в общую дисперсию

pc <- princomp(scale(data))
summary(pc)

Importance of components:
                          Comp.1    Comp.2    Comp.3    Comp.4    Comp.5
Standard deviation     1.2207259 1.0501082 0.9385849 0.8843622 0.8520637
Proportion of Variance 0.2991142 0.2213445 0.1768267 0.1569860 0.1457286
Cumulative Proportion  0.2991142 0.5204587 0.6972854 0.8542714 1.0000000

Первые две компоненты дают 52.04587%

Факторные нагрузки

pc$loadings[, seq(5)]

             Comp.1     Comp.2      Comp.3      Comp.4      Comp.5
sstati    0.5088674  0.3666091  0.13377642  0.58353704  0.49823700
asi1_med  0.4316568  0.5052637  0.07655432 -0.74242863  0.03633432
asi2_emp  0.2806453 -0.7246625  0.22612083 -0.28145720  0.51551306
asid3_dyr 0.4233503 -0.1845726 -0.88210384  0.02520452 -0.08924714
asi7_psy  0.5447377 -0.2260610  0.38341295  0.16861444 -0.69044981

Первая компонента отвечает за “медицинские” факторы: тревожность, медицинский и психиатрический статус, длительность героиновой зависимости, притом везде связь прямо пропорциональна.
Вторая компонента связана с понижением социального статуса и повышением медицинского статуса.
Третья – с понижением длительности героиновой зависимости.
Четвертая – с понижением медицинского статуса и повышением уровня тревожности.
Пятая компонента говорит о повышении уровня тревожности и социального статуса вместе с понижением психиатрического статуса.

Проверка [совпадают с точностью до перестановки знаков для, например, первого фактора]:

ei <- eigen(cor(data))

A <- ei$vectors
Y <- as.matrix(scale(data)) %*% A

cbind(Y[seq(10), seq(3)], pc$scores[seq(10),seq(3)])

                                             Comp.1      Comp.2      Comp.3
 [1,] -0.1983072  0.72089292  0.19465100  0.1983072  0.72089292  0.19465100
 [2,]  1.4273690  2.44695870  0.06755617 -1.4273690  2.44695870  0.06755617
 [3,] -1.0878609  0.12771542 -0.30632264  1.0878609  0.12771542 -0.30632264
 [4,]  0.2421745  0.06698835  0.61873096 -0.2421745  0.06698835  0.61873096
 [5,] -1.2291057 -0.43829516  0.62489850  1.2291057 -0.43829516  0.62489850
 [6,] -1.8404684  1.85516675  0.84688457  1.8404684  1.85516675  0.84688457
 [7,] -0.5164741 -0.90396856  0.22615760  0.5164741 -0.90396856  0.22615760
 [8,]  0.6092661  1.22103938  0.14656096 -0.6092661  1.22103938  0.14656096
 [9,]  0.5496977 -0.84859723  0.39636680 -0.5496977 -0.84859723  0.39636680
[10,] -1.4000370 -0.80176615  0.17534824  1.4000370 -0.80176615  0.17534824

plot(pc$scores, type="n")
text(pc$scores, row.names(data), cex=0.6)

Заметим выброс в виде пациента №166: для него характерно большое значение Comp.1 – повышенные “медицинские” факторы – при достаточно низком Comp.2