multifactor dimensionality reduction Опции

[_alena_]

всем привет. очень нуждаюсь в помощи и советах.... кто то знаком с multifactor dimensionality reduction??? что можно почитать о том как работать в этой среде? для каких видов данных она предназначена??? а то мне все что встречалось касается генной информации. можно ли применять другие выборки? спасибо

[_alena_]

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

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

[p2004r]

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

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

лучше больше случаев видеть конечно.

переменная Statiny : что означают её уровни?

[_alena_]

лучше больше случаев видеть конечно.

переменная Statiny : что означают её уровни?

сейчас загружу остальное... но я ж говорю мне бы понять самой что делать да как....
Статины:
1-не применяются
2-применяются (неэффективная доза)
3-применяются (эффективная доза)
4-прекратили прием (эффективная доза)
5-прекратили прием (неэффективная доза)

[p2004r]

сейчас загружу остальное... но я ж говорю мне бы понять самой что делать да как....
Статины:
1-не применяются
2-применяются (неэффективная доза)
3-применяются (эффективная доза)
4-прекратили прием (эффективная доза)
5-прекратили прием (неэффективная доза)

хорошо еще и его учтем сейчас,

по неполным данным логит регрессия полностью разделяет Progress

минимальная модель --- glm(formula = Progres ~ OP + Shunt + ChSS + AG.therapia + cholesterin,
family = binomial(logit), data = data.sh.f)


сейчас расчеты закончу и прокомментирую

[p2004r]

хорошо еще и его учтем сейчас,

сейчас расчеты закончу и прокомментирую

ха, над статинами думать надо как закодировать, так сходу не придумалось

Код

# читаем данные в таблицу data
data<-read.csv2("1.csv")

## сменим уровни шкал в таблицу data.sh

# копируем оригинальные данные
data.sh <- data

data.sh$Progres[data$Progres==1]<-0  # перекодируем все 1 по показателю Progres в
                                                          # исходной таблице data в 0 таблицы data.sh
data.sh$Progres[data$Progres==2]<-1  # тоже для 2, заменяем их на  1

data.sh$AIK[data$AIK==1]<-0
data.sh$AIK[data$AIK==2]<-1

data.sh$CA[data$CA==1]<-0
data.sh$CA[data$CA==2]<-1

data.sh$Progres[data$Progres==1]<-0
data.sh$Progres[data$Progres==2]<-1

data.sh$OP[data$OP==1]<-0
data.sh$OP[data$OP==2]<-1

data.sh$Shunt[data$Shunt==1]<-0
data.sh$Shunt[data$Shunt==2]<-1

data.sh$pol[data$pol==1]<-0
data.sh$pol[data$pol==2]<-1

data.sh$AG.therapia[data$AG.therapia==1]<-0
data.sh$AG.therapia[data$AG.therapia==2]<-1

data.sh$stenoz.a.renalis[data$stenoz.a.renalis==1]<-0
data.sh$stenoz.a.renalis[data$stenoz.a.renalis==2]<-1

data.sh$stenoz.arterij.nog[data$stenoz.arterij.nog==1]<-0
data.sh$stenoz.arterij.nog[data$stenoz.arterij.nog==2]<-1

data.sh$insult.v.anamnese[data$insult.v.anamnese==1]<-0
data.sh$insult.v.anamnese[data$insult.v.anamnese==2]<-1

data.sh$diabetus.melitus[data$diabetus.melitus==1]<-0
data.sh$diabetus.melitus[data$diabetus.melitus==2]<-1


# заменяем tip на три переменных, каждая из которых содержит признак только одного типа кровобращения
tip.1<- rep(0,length(data.sh$tip.krovoobraschenija))
tip.2<- rep(0,length(data.sh$tip.krovoobraschenija))
tip.3<- rep(0,length(data.sh$tip.krovoobraschenija))
tip.1[data.sh$tip.krovoobraschenija==1]<-1
tip.2[data.sh$tip.krovoobraschenija==2]<-1
tip.3[data.sh$tip.krovoobraschenija==3]<-1


# склеиваем новые переменные тип в одну таблицу с трансформированными данными
data.sh.f<- cbind(data.sh, tip.1, tip.2, tip.3)

# определяем модель логистической регрессии
model<- glm(Progres~AIK+OP+Shunt+vozrast+IMT+Stepen.Oz.+ChSS+AD+AG.therapia+cholesterin+diabetus.m
elitus+tip.1+tip.2+tip.3,
            family=binomial(logit),
            data.sh.f)

# оптимизируем модель по AIC, сокращая число предикторов до минимального
model.f<-step(model)

> model.f

Call:  glm(formula = Progres ~ Shunt + cholesterin + diabetus.melitus +
    tip.1 + tip.2, family = binomial(logit), data = data.sh.f)

Coefficients:
     (Intercept)             Shunt       cholesterin  diabetus.melitus  
         -0.9145            1.8652           -0.3532            1.5293  
           tip.1             tip.2  
          1.5598            1.4139  

Degrees of Freedom: 61 Total (i.e. Null);  56 Residual
Null Deviance:        85.69
Residual Deviance: 67.62     AIC: 79.62



library(lattice)

# как разделяет полная модель
densityplot( ~ predict(model,
                       type = "response"),
            groups = data.sh.f$Progres,
            plot.points = "rug",
            auto.key = TRUE)

# как разделяет оптимальная по AIC
densityplot( ~ predict(model.f,
                       type = "response"),
            groups = data.sh.f$Progres,
            plot.points = "rug",
            auto.key = TRUE)

# оцениваем насколько модели могут что то предсказывать


library(pROC)

roc.model <- roc(data.sh.f$Progres,               # данные из выборки
                        predict(model, type="response")) # предсказание модели

roc.model.f <- roc(data.sh.f$Progres,                 # данные из выборки
                   predict(model.f, type="response")) # предсказание модели

# Рисуем ROC кривые для моделей
> plot(roc.model)
> lines.roc(roc.model.f, col="red")
> auc(roc.model.f)
Area under the curve: 0.7931
> auc(roc.model)
Area under the curve: 0.8286

Вот как то так если логистическая регрессия используется


вот протокол исключения параметров из модели

Код

Start:  AIC=90.87
Progres ~ AIK + OP + Shunt + vozrast + IMT + Stepen.Oz. + ChSS +
    AD + AG.therapia + cholesterin + diabetus.melitus + tip.1 +
    tip.2 + tip.3


Step:  AIC=90.87
Progres ~ AIK + OP + Shunt + vozrast + IMT + Stepen.Oz. + ChSS +
    AD + AG.therapia + cholesterin + diabetus.melitus + tip.1 +
    tip.2

                   Df Deviance    AIC
- Stepen.Oz.        1   62.873 88.873
- IMT               1   62.987 88.987
- ChSS              1   63.034 89.034
- AD                1   63.407 89.407
- AIK               1   63.509 89.509
- OP                1   63.598 89.598
- vozrast           1   63.662 89.662
- tip.2             1   63.776 89.776
- diabetus.melitus  1   63.925 89.925
- AG.therapia       1   64.572 90.572
- tip.1             1   64.742 90.742
<none>                  62.869 90.869
- cholesterin       1   67.407 93.407
- Shunt             1   67.655 93.655

Step:  AIC=88.87
Progres ~ AIK + OP + Shunt + vozrast + IMT + ChSS + AD + AG.therapia +
    cholesterin + diabetus.melitus + tip.1 + tip.2

                   Df Deviance    AIC
- ChSS              1   63.037 87.037
- IMT               1   63.389 87.389
- AD                1   63.409 87.409
- AIK               1   63.589 87.589
- OP                1   63.608 87.608
- vozrast           1   63.662 87.662
- tip.2             1   63.787 87.787
- diabetus.melitus  1   63.954 87.954
- AG.therapia       1   64.659 88.659
- tip.1             1   64.859 88.859
<none>                  62.873 88.873
- cholesterin       1   67.409 91.409
- Shunt             1   67.658 91.658

Step:  AIC=87.04
Progres ~ AIK + OP + Shunt + vozrast + IMT + AD + AG.therapia +
    cholesterin + diabetus.melitus + tip.1 + tip.2

                   Df Deviance    AIC
- AD                1   63.456 85.456
- IMT               1   63.640 85.640
- OP                1   63.721 85.721
- AIK               1   63.746 85.746
- tip.2             1   63.894 85.894
- vozrast           1   63.984 85.984
- diabetus.melitus  1   64.001 86.001
- AG.therapia       1   64.764 86.764
- tip.1             1   64.882 86.882
<none>                  63.037 87.037
- cholesterin       1   67.423 89.423
- Shunt             1   68.079 90.079

Step:  AIC=85.46
Progres ~ AIK + OP + Shunt + vozrast + IMT + AG.therapia + cholesterin +
    diabetus.melitus + tip.1 + tip.2

                   Df Deviance    AIC
- IMT               1   63.885 83.885
- AIK               1   64.043 84.043
- OP                1   64.053 84.053
- diabetus.melitus  1   64.535 84.535
- vozrast           1   64.610 84.610
- AG.therapia       1   64.863 84.863
- tip.2             1   64.873 84.873
<none>                  63.456 85.456
- tip.1             1   65.866 85.866
- cholesterin       1   67.721 87.721
- Shunt             1   68.349 88.349

Step:  AIC=83.89
Progres ~ AIK + OP + Shunt + vozrast + AG.therapia + cholesterin +
    diabetus.melitus + tip.1 + tip.2

                   Df Deviance    AIC
- AIK               1   64.305 82.305
- OP                1   64.479 82.479
- vozrast           1   65.102 83.102
- AG.therapia       1   65.497 83.497
- diabetus.melitus  1   65.546 83.546
- tip.2             1   65.785 83.785
<none>                  63.885 83.885
- tip.1             1   66.413 84.413
- cholesterin       1   67.830 85.830
- Shunt             1   68.972 86.972

Step:  AIC=82.31
Progres ~ OP + Shunt + vozrast + AG.therapia + cholesterin +
    diabetus.melitus + tip.1 + tip.2

                   Df Deviance    AIC
- OP                1   64.900 80.900
- AG.therapia       1   65.680 81.680
- vozrast           1   65.972 81.972
<none>                  64.305 82.305
- diabetus.melitus  1   66.572 82.572
- tip.1             1   67.013 83.013
- tip.2             1   67.361 83.361
- cholesterin       1   68.211 84.211
- Shunt             1   70.676 86.676

Step:  AIC=80.9
Progres ~ Shunt + vozrast + AG.therapia + cholesterin + diabetus.melitus +
    tip.1 + tip.2

                   Df Deviance    AIC
- AG.therapia       1   65.998 79.998
- vozrast           1   66.623 80.623
<none>                  64.900 80.900
- diabetus.melitus  1   66.981 80.981
- tip.1             1   67.535 81.535
- tip.2             1   67.586 81.586
- cholesterin       1   68.635 82.635
- Shunt             1   70.988 84.988

Step:  AIC=80
Progres ~ Shunt + vozrast + cholesterin + diabetus.melitus +
    tip.1 + tip.2

                   Df Deviance    AIC
- vozrast           1   67.619 79.619
<none>                  65.998 79.998
- tip.2             1   68.893 80.893
- diabetus.melitus  1   69.058 81.058
- tip.1             1   69.349 81.349
- cholesterin       1   69.964 81.964
- Shunt             1   71.574 83.574

Step:  AIC=79.62
Progres ~ Shunt + cholesterin + diabetus.melitus + tip.1 + tip.2

                   Df Deviance    AIC
<none>                  67.619 79.619
- cholesterin       1   70.782 80.782
- tip.2             1   70.954 80.954
- diabetus.melitus  1   71.087 81.087
- tip.1             1   71.458 81.458
- Shunt             1   73.848 83.848
> model.f

Call:  glm(formula = Progres ~ Shunt + cholesterin + diabetus.melitus +
    tip.1 + tip.2, family = binomial(logit), data = data.sh.f)

Coefficients:
     (Intercept)             Shunt       cholesterin  diabetus.melitus  
         -0.9145            1.8652           -0.3532            1.5293  
           tip.1             tip.2  
          1.5598            1.4139  

Degrees of Freedom: 61 Total (i.e. Null);  56 Residual
Null Deviance:        85.69
Residual Deviance: 67.62     AIC: 79.62

Эскизы прикрепленных изображений

 

[YVR]

вот протокол исключения параметров из модели

Код

Progres ~ Shunt + cholesterin + diabetus.melitus + tip.1 + tip.2

Что интересно у меня нелинейная модель, но попробовал два варианта:

1. Взял 9 первых в таблице параметров (кроме CA)
2. Оставшиеся после исключения 5 параметров, как получилось у Вас

По п. 1 модель дала сепарабельность 96.77% т.е. в 60 из 62 наличие или отсутствие Прогресса было вычислено правильно. Теоретически подогнать модель можно и под 100% сепарабельность, но практически нельзя, т.к. модель со 100% сепарабельностью, полученная на обучающей выборке (31 пример) нестационарна. Т.е. на контрольной выборке (тоже 31 пример) имеет значительно отличающиеся матожидание и дисперсию, вычисленные по остаткам на выборке обучающей. А нестационарная модель нерепродуктивна, т.е. может быть опровергнута независимыми экспериментами.

По п. 2 модель дала сепарабельность 75.81% т.е. в 47 из 62 случаев наличие или отсутствие Прогресса было вычислено правильно. Но модель нестационарна, даже по количеству распознанных Прогрессов - на обучающей выборке их 32, а на контрольной только 15. Да впрочем почему модель нестационарна понять не сложно: все признаки кроме холестерина бинарны и распределены неравномерно, т.е. добиться более или менее приемлемой стационарности при таких данных можно только на очень больших выборках.

Сообщение отредактировал YVR - 21.03.2012 - 20:51

[p2004r]

Что интересно у меня нелинейная модель, но попробовал два варианта:

1. Взял 9 первых в таблице параметров (кроме CA)
2. Оставшиеся после исключения 5 параметров, как получилось у Вас

По п. 1 модель дала сепарабельность 96.77% т.е. в 60 из 62 наличие или отсутствие Прогресса было вычислено правильно. Теоретически подогнать модель можно и под 100% сепарабельность, но практически нельзя, т.к. модель со 100% сепарабельностью, полученная на обучающей выборке (31 пример) нестационарна. Т.е. на контрольной выборке (тоже 31 пример) имеет значительно отличающиеся матожидание и дисперсию, вычисленные по остаткам на выборке обучающей. А нестационарная модель нерепродуктивна, т.е. может быть опровергнута независимыми экспериментами.

По п. 2 модель дала сепарабельность 75.81% т.е. в 47 из 62 случаев наличие или отсутствие Прогресса было вычислено правильно. Но модель нестационарна, даже по количеству распознанных Прогрессов - на обучающей выборке их 32, а на контрольной только 15. Да впрочем почему модель нестационарна понять не сложно: все признаки кроме холестерина бинарны и распределены неравномерно, т.е. добиться более или менее приемлемой стационарности при таких данных можно только на очень больших выборках.

Все можно проверить бутстрепом

Код

> ci.auc(roc.model.f)
95% CI: 0.6836-0.9026 (DeLong)
> ci.auc(roc.model.f,method="bootstrap")
95% CI: 0.6677-0.8945 (2000 stratified bootstrap replicates)
> ci.auc(roc.model)
95% CI: 0.7209-0.926 (DeLong)
> ci.auc(roc.model,method="bootstrap")
95% CI: 0.7126-0.9164 (2000 stratified bootstrap replicates)
>

если коэффициенты модели, то они тоже могут быть оценены

Код

r.boot <- t(replicate(200,
                      as.numeric((glm(Progres ~ Shunt + cholesterin + diabetus.melitus + tip.1 + tip.2,
                                      family=binomial(logit),
                                      data.sh.f[sample.int(length(data.sh.f[,1]),
                                                           replace = TRUE),]))$coefficients)))

parallel(r.boot, horizontal.axis=FALSE)

безусловно случаев крайне немного даже для сколько нибудь уверенных выводах даже для линейной модели

Эскизы прикрепленных изображений

 

[p2004r]

Ну нелинейно так нелинейно

Итак randomForest, он сразу пишет бутстреп оценку насколько можно доверять его результатам (естественно обучающую выборку пилит со 100% точностью)

Код

> names(data.sh.f)
[1] "X..bolnogo"            "AIK"                   "CA"                  
[4] "Progres"               "OP"                    "Shunt"                
[7] "Srok.gruppy"           "vozrast"               "IMT"                  
[10] "Stepen.Oz."            "pol"                   "tip.krovoobraschenija"
[13] "ChSS"                  "AD"                    "AG.therapia"          
[16] "Statiny"               "cholesterin"           "stenoz.a.renalis"    
[19] "stenoz.arterij.nog"    "insult.v.anamnese"     "diabetus.melitus"    
[22] "tip.1"                 "tip.2"                 "tip.3"                
[25] "statin.ter"            "statin.break"          "statin.doz"    

> randomForest(data.sh.f[,-c(1,4,16,12)], y=factor(data.sh.f$Progres))

Call:
randomForest(x = data.sh.f[, -c(1, 4, 16, 12)], y = factor(data.sh.f$Progres))
               Type of random forest: classification
                     Number of trees: 500
No. of variables tried at each split: 4

        OOB estimate of  error rate: 40.32%
Confusion matrix:
   0  1 class.error
0 15 14   0.4827586
1 11 22   0.3333333
>

важность переменных

> randomForest(data.sh.f[,c(7,8,9,17)], y=factor(data.sh.f$Progres))

Call:
randomForest(x = data.sh.f[, c(7, 8, 9, 17)], y = factor(data.sh.f$Progres))
               Type of random forest: classification
                     Number of trees: 500
No. of variables tried at each split: 2

        OOB estimate of  error rate: 40.32%
Confusion matrix:
   0  1 class.error
0 18 11   0.3793103
1 14 19   0.4242424

varImpPlot(randomForest(data.sh.f[,-c(1,4,16,12)], y=factor(data.sh.f$Progres)), type=2)

ну типа вот нелинейная важность данных для классификации.

присоединим картинку как он классифицирует по первым самым важным.

а не такая плохая модель

Код

> m.rf <-glm(Progres ~ Srok.gruppy + cholesterin + IMT + vozrast ,
+            family=binomial(logit),
+            data.sh.f)


> roc.m.rf.f <- roc(data.sh.f$Progres,               # данные из выборки
+                   predict(m.rf, type="response")) # предсказание модели


> densityplot( ~ predict(m.rf,
+                        type = "response"),
+             groups = data.sh.f$Progres,
+             plot.points = "rug",
+             auto.key = TRUE)


> plot(roc.m.rf.f)

Call:
roc.default(response = data.sh.f$Progres, predictor = predict(m.rf,     type = "response"))

Data: predict(m.rf, type = "response") in 29 controls (data.sh.f$Progres 0) < 33 cases (data.sh.f$Progres 1).
Area under the curve: 0.7827

> ci.auc(roc.m.rf.f)
95% CI: 0.6579-0.9074 (DeLong)

Сообщение отредактировал p2004r - 23.03.2012 - 17:57

Эскизы прикрепленных изображений