## ---
## biotools version 3.0

Данные

df <- read.table("/home/ekaterina/Downloads/newthyroid.txt", sep = ",")
colnames(df) <- c("T3resin", "Thyroxin", "Triiodothyronine", "Thyroidstimulating", "TSH_valu", "Class")
head(df)
##   T3resin Thyroxin Triiodothyronine Thyroidstimulating TSH_valu Class
## 1     107     10.1              2.2                0.9      2.7     1
## 2     113      9.9              3.1                2.0      5.9     1
## 3     127     12.9              2.4                1.4      0.6     1
## 4     109      5.3              1.6                1.4      1.5     1
## 5     105      7.3              1.5                1.5     -0.1     1
## 6     105      6.1              2.1                1.4      7.0     1

Признаки - это названия гармонов щитовидной железы. Class имеет 3 категории:

  1. Норма

  2. Гипертиреоз

  3. Гипотиреоз

Наша цель выяснить отличаются ли группы между собой и чем можно объяснить это различие.

pairs.panels(df, lm = TRUE)

Прологорифмируем некоторые признаки

df$Thyroxin <- log(df$Thyroxin)
df$Triiodothyronine <- log(df$Triiodothyronine)
df$Thyroidstimulating <- log(df$Thyroidstimulating)
pairs.panels(df, lm = TRUE)

Достаточно большие корреляции между некоторыми переменными можно объяснить с точки зрения медицины.

Проверка модели

Нормальность

#Group 1
s <- subset(df, df$Class == 1)
mshapiro.test(t(as.matrix(s[,-6])))
## 
##  Shapiro-Wilk normality test
## 
## data:  Z
## W = 0.90745, p-value = 3.562e-08
#Group 2
s <- subset(df, df$Class == 2)
mshapiro.test(t(as.matrix(s[,-6])))
## 
##  Shapiro-Wilk normality test
## 
## data:  Z
## W = 0.84216, p-value = 0.0001539
#Group 3
s <- subset(df, df$Class == 3)
mshapiro.test(t(as.matrix(s[,-6])))
## 
##  Shapiro-Wilk normality test
## 
## data:  Z
## W = 0.87098, p-value = 0.001762

Гипотеза о нормальности отвергается для всех групп

Гомоскедастичность

hom <- boxM(Y,type)
hom
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  Y
## Chi-Sq (approx.) = 704.47, df = 30, p-value < 2.2e-16

Гипотеза о равенстве ковариационных матриц отвергается

MANOVA

manova2 <- manova(as.matrix(df[-6])~df$Class, data = df)
summary(manova2, test = "Wilks")
##            Df   Wilks approx F num Df den Df    Pr(>F)    
## df$Class    2 0.13378   72.137     10    416 < 2.2e-16 ***
## Residuals 212                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Дискриминантный анализ значим

Все это верно для нормальной модели.

Расстояние Махаланобиса между группами

pooled_cov <- hom$pooled
mach <- pairwise.mahalanobis(Y, grouping = type, cov = pooled_cov)
print(mach$distance)
##          [,1]     [,2]     [,3]
## [1,]  0.00000  8.26333 25.87848
## [2,]  8.26333  0.00000 55.07655
## [3,] 25.87848 55.07655  0.00000

Наибольшее расстояние между 2 и 3 группой, наименьшее между 1 и 2.

Hotelling test

#Между 1 и 2 группой
print(hotelling.test(.~Class, data = df, pair = c(1,2)))
## Test stat:  80.698 
## Numerator df:  5 
## Denominator df:  179 
## P-value:  0
#Между 1 и 3 группой
print(hotelling.test(.~Class, data = df, pair = c(1,3)))
## Test stat:  116.35 
## Numerator df:  5 
## Denominator df:  174 
## P-value:  0
#Между 2 и 3 группой
print(hotelling.test(.~Class, data = df, pair = c(2,3)))
## Test stat:  68.808 
## Numerator df:  5 
## Denominator df:  59 
## P-value:  0

Все p-value равны 0, значит все группы между сообой значимо различаются

Избыточность признаков

Выберем признаки, которые влияют на качество разделения. Forward Greedy Wilks (Stepwise Analysis).

type <- as.numeric(type)
greedy.wilks(df[-6], type)
## Formula containing included variables: 
## 
## type ~ Thyroxin + TSH_valu + Thyroidstimulating + Triiodothyronine + 
##     T3resin
## <environment: 0x68b9080>
## 
## 
## Values calculated in each step of the selection procedure: 
## 
##                 vars Wilks.lambda F.statistics.overall p.value.overall
## 1           Thyroxin    0.2999984            247.33526    3.754935e-56
## 2           TSH_valu    0.1981744            131.48925    8.140747e-73
## 3 Thyroidstimulating    0.1707704             99.39156    1.962108e-77
## 4   Triiodothyronine    0.1473669             83.85882    4.627563e-82
## 5            T3resin    0.1337772             72.13712    1.863610e-84
##   F.statistics.diff p.value.diff
## 1         247.33526 3.754935e-56
## 2          54.20691 0.000000e+00
## 3          16.84967 1.625689e-07
## 4          16.59572 2.034494e-07
## 5          10.56480 4.257579e-05

Получили, что все признаки оказывают значимое влияние на результат.

Canonical Discriminant Analysis

cd <- candisc(manova2, data = df)
summary(cd)
## 
## Canonical Discriminant Analysis for df$Class:
## 
##    CanRsq Eigenvalue Difference Percent Cumulative
## 1 0.81441    4.38815     4.0008 91.8893     91.889
## 2 0.27919    0.38733     4.0008  8.1107    100.000
## 
## Class means:
## 
##       Can1     Can2
## 1 -0.29954  0.39696
## 2 -2.72641 -1.14370
## 3  4.67854 -0.65049
## 
##  std coefficients:
##                        Can1     Can2
## T3resin             0.23510  0.41904
## Thyroxin           -0.67705 -0.01736
## Triiodothyronine   -0.14526 -0.66284
## Thyroidstimulating  0.25937 -0.63111
## TSH_valu            0.52457 -0.14436
cd$coeffs.std
##                          Can1       Can2
## T3resin             0.2350990  0.4190447
## Thyroxin           -0.6770482 -0.0173597
## Triiodothyronine   -0.1452639 -0.6628392
## Thyroidstimulating  0.2593749 -0.6311149
## TSH_valu            0.5245687 -0.1443624
cd$structure
##                          Can1       Can2
## T3resin             0.5711575  0.4249287
## Thyroxin           -0.9236947 -0.1356796
## Triiodothyronine   -0.7224633 -0.5375347
## Thyroidstimulating  0.8416865 -0.3776452
## TSH_valu            0.7230735 -0.3387810

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

Значимое число канонических переменных

Lambda Prime Test

print(cd)
## 
## Canonical Discriminant Analysis for df$Class:
## 
##    CanRsq Eigenvalue Difference Percent Cumulative
## 1 0.81441    4.38815     4.0008 91.8893     91.889
## 2 0.27919    0.38733     4.0008  8.1107    100.000
## 
## Test of H0: The canonical correlations in the 
## current row and all that follow are zero
## 
##   LR test stat approx F numDF denDF   Pr(> F)    
## 1      0.13378   72.137    10   416 < 2.2e-16 ***
## 2      0.72081   20.238     4   209 4.188e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Все канонические переменные значимы

plot(cd)
## Vector scale factor set to 4.81

По графику видно, что первую каноническую переменную можно инерпретировать как различие между группами в целом. У людей из 3 группы (Гипотиреоз, зеленый цвет) повышены Thyroidstimulating, T3resin и Thyroidstimulating. У людей и 2 группы (гипертиреоз, красный цвет) повышены Thyroxin и Triiodothyronine. Люди без выраженных отколнений в функциях щитовидной железы скучковались у нуля, то есть у них значения гармонов не привышвет норму.

По второй канонической переменной описывается разница между гармоном T3resin и всеми остальными. И можно предположить, что это различие между больными и здоровыми.