Przykład: regresja LASSO z użyciem glmnet dla wyniku binarnego

77

Zaczynam bawić sięglmnet za pomocą regresji LASSO, gdzie moje wyniki zainteresowania są dychotomiczne. Poniżej utworzyłem małą próbną ramkę danych:

age     <- c(4, 8, 7, 12, 6, 9, 10, 14, 7) 
gender  <- c(1, 0, 1, 1, 1, 0, 1, 0, 0)
bmi_p   <- c(0.86, 0.45, 0.99, 0.84, 0.85, 0.67, 0.91, 0.29, 0.88)
m_edu   <- c(0, 1, 1, 2, 2, 3, 2, 0, 1)
p_edu   <- c(0, 2, 2, 2, 2, 3, 2, 0, 0)
f_color <- c("blue", "blue", "yellow", "red", "red", "yellow", "yellow", 
             "red", "yellow")
asthma  <- c(1, 1, 0, 1, 0, 0, 0, 1, 1)
# df is a data frame for further use!
df <- data.frame(age, gender, bmi_p, m_edu, p_edu, f_color, asthma)

Kolumny (zmienne) w powyższym zestawie danych są następujące:

  • age (wiek dziecka w latach) - ciągły
  • gender - binarny (1 = męski; 0 = żeński)
  • bmi_p (Percentyl BMI) - ciągły
  • m_edu (najwyższy poziom wykształcenia matki) - zwykły (0 = niższy niż liceum; 1 = dyplom ukończenia szkoły średniej; 2 = licencjat; 3 = licencjat)
  • p_edu (najwyższy poziom wykształcenia ojca) - porządkowy (taki sam jak m_edu)
  • f_color (ulubiony kolor podstawowy) - nominalny („niebieski”, „czerwony” lub „żółty”)
  • asthma (stan astmy dziecięcej) - dwójkowy (1 = astma; 0 = brak astmy)

Celem tego przykładu jest wykorzystanie Lasso do stworzenia modelu prognozowania stanu astmy dziecko z listy 6 potencjalnych zmiennych objaśniających ( age, gender, bmi_p, m_edu, p_edu, i f_color). Oczywiście wielkość próby jest tutaj problemem, ale mam nadzieję uzyskać lepszy wgląd w sposób obsługi różnych typów zmiennych (tj. Ciągły, porządkowy, nominalny i binarny) w glmnetramach, gdy wynik jest binarny (1 = astma ; 0 = bez astmy).

W związku z tym, czy ktoś byłby skłonny dostarczyć przykładowy Rskrypt wraz z objaśnieniami dla tego fałszywego przykładu, używając LASSO z powyższymi danymi do przewidzenia stanu astmy? Chociaż bardzo prosty, wiem, że i prawdopodobnie wielu innych na CV, bardzo by to docenił!

Matt Reichenbach
źródło
2
Może masz więcej szczęścia, jeśli wysłane dane jako dputo rzeczywistej obiektu R; nie każ czytelnikom nakładać szronu na wierzch, a także upiec ci ciasto! Jeśli wygenerujesz odpowiednią ramkę danych w R, powiedzmy foo, następnie edytuj w pytaniu wynik dput(foo).
Gavin Simpson
Dzięki @GavinSimpson! Zaktualizowałem post za pomocą ramki danych, więc mam nadzieję, że zjem trochę ciasta bez szronu! :)
Matt Reichenbach,
2
Używając percentyla BMI, w pewnym sensie przeczysz prawom fizyki. Otyłość wpływa na osoby zgodnie z pomiarami fizycznymi (długości, objętości, waga), a nie w zależności od liczby osób podobnych do obecnego pacjenta, co robi procentowanie.
Frank Harrell,
3
Zgadzam się, percentyl BMI nie jest metryką, którą wolę wykorzystywać; Wytyczne CDC zalecają jednak stosowanie percentyla BMI zamiast BMI (również wysoce wątpliwy wskaźnik!) dla dzieci i młodzieży w wieku poniżej 20 lat, ponieważ uwzględnia wiek i płeć oprócz wzrostu i masy ciała. Wszystkie te zmienne i wartości danych zostały całkowicie przemyślane na potrzeby tego przykładu. Ten przykład nie odzwierciedla żadnej z moich bieżących prac, ponieważ pracuję z dużymi zbiorami danych. Chciałem tylko zobaczyć przykład glmnetdziałania z podwójnym wynikiem.
Matt Reichenbach,
Podłącz tutaj pakiet Patricka Breheny'ego o nazwie ncvreg, który pasuje do modeli regresji liniowej i logistycznej karanych przez MCP, SCAD lub LASSO. ( cran.r-project.org/web/packages/ncvreg/index.html )
bdeonovic

Odpowiedzi:

100
library(glmnet)

age     <- c(4, 8, 7, 12, 6, 9, 10, 14, 7) 
gender  <- as.factor(c(1, 0, 1, 1, 1, 0, 1, 0, 0))
bmi_p   <- c(0.86, 0.45, 0.99, 0.84, 0.85, 0.67, 0.91, 0.29, 0.88) 
m_edu   <- as.factor(c(0, 1, 1, 2, 2, 3, 2, 0, 1))
p_edu   <- as.factor(c(0, 2, 2, 2, 2, 3, 2, 0, 0))
f_color <- as.factor(c("blue", "blue", "yellow", "red", "red", "yellow", 
                       "yellow", "red", "yellow"))
asthma <- c(1, 1, 0, 1, 0, 0, 0, 1, 1)

xfactors <- model.matrix(asthma ~ gender + m_edu + p_edu + f_color)[, -1]
x        <- as.matrix(data.frame(age, bmi_p, xfactors))

# Note alpha=1 for lasso only and can blend with ridge penalty down to
# alpha=0 ridge only.
glmmod <- glmnet(x, y=as.factor(asthma), alpha=1, family="binomial")

# Plot variable coefficients vs. shrinkage parameter lambda.
plot(glmmod, xvar="lambda")

wprowadź opis zdjęcia tutaj

Zmienne jakościowe są zwykle najpierw przekształcane w czynniki, a następnie tworzona jest sztuczna macierz zmiennych predyktorów i wraz z predyktorami ciągłymi przekazywana jest do modelu. Pamiętaj, że glmnet stosuje zarówno kary kalenicowe, jak i lasso, ale można ustawić je oddzielnie.

Niektóre wyniki:

# Model shown for lambda up to first 3 selected variables.
# Lambda can have manual tuning grid for wider range.

glmmod
# Call:  glmnet(x = x, y = as.factor(asthma), family = "binomial", alpha = 1) 
# 
#        Df    %Dev   Lambda
#   [1,]  0 0.00000 0.273300
#   [2,]  1 0.01955 0.260900
#   [3,]  1 0.03737 0.249000
#   [4,]  1 0.05362 0.237700
#   [5,]  1 0.06847 0.226900
#   [6,]  1 0.08204 0.216600
#   [7,]  1 0.09445 0.206700
#   [8,]  1 0.10580 0.197300
#   [9,]  1 0.11620 0.188400
#  [10,]  3 0.13120 0.179800
#  [11,]  3 0.15390 0.171600
# ...

Współczynniki można uzyskać z glmmod. Tutaj pokazano z wybranymi 3 zmiennymi.

coef(glmmod)[, 10]
#   (Intercept)           age         bmi_p       gender1        m_edu1 
#    0.59445647    0.00000000    0.00000000   -0.01893607    0.00000000 
#        m_edu2        m_edu3        p_edu2        p_edu3    f_colorred 
#    0.00000000    0.00000000   -0.01882883    0.00000000    0.00000000 
# f_coloryellow 
#   -0.77207831 

Wreszcie, do wyboru lambda można również użyć krzyżowej weryfikacji.

cv.glmmod <- cv.glmnet(x, y=asthma, alpha=1)
plot(cv.glmmod)

wprowadź opis zdjęcia tutaj

(best.lambda <- cv.glmmod$lambda.min)
# [1] 0.2732972
poklepać
źródło
4
to jest dokładnie to, czego szukałem +1, jedyne pytania, jakie mam, to 1) co możesz zrobić z lambda walidacji krzyżowej 0,2732972? oraz 2) Czy z glmmod wybrane zmienne są ulubionymi kolorami (żółty), płcią i wykształceniem ojca (stopień licencjata)? Dzięki wielkie!
Matt Reichenbach,
4
1) Walidacja krzyżowa służy do wyboru lambda i współczynników (przy minimalnym błędzie). W tym makiecie nie ma lokalnego min (było ostrzeżenie związane również ze zbyt małą liczbą obs); Zinterpretowałbym, że wszystkie współczynniki zostały zmniejszone do zera za pomocą kar skurczowych (najlepszy model ma tylko przechwytywanie) i wznowiłbym z większą liczbą (rzeczywistych) obserwacji i może zwiększyć zakres lambda. 2) Tak, w przykładzie, w którym wybrałem coef (glmmod) [, 10] ... wybierasz lambda dla modelu poprzez CV lub interpretację wyników. Czy możesz oznaczyć jako rozwiązany, jeśli uważasz, że rozwiązałem twoje pytanie? dzięki.
klep.
2
Czy mogę zapytać, w jaki sposób obsługuje f_colorzmienną? Czy współczynnik poziomu 1–4 jest uważany za większy krok niż 1–2, czy też wszystkie są jednakowo ważone, bezkierunkowe i kategoryczne? (Chcę zastosować to do analizy ze wszystkimi nieuporządkowanymi predyktorami.)
beroe
3
Linia xfactors <- model.matrix(asthma ~ gender + m_edu + p_edu + f_color)[,-1]koduje zmienną kategorialną f_kolor (zadeklarowaną przez as.factorw poprzednich wierszach). Powinien używać domyślnego kodowania zmiennych fikcyjnych R, chyba że contrasts.argpodano argument. Oznacza to, że wszystkie poziomy f_kolor są jednakowo ważone i bezkierunkowe, z wyjątkiem pierwszego, który jest używany jako klasa referencyjna i wchłaniany do przecięcia.
Alex
1
@Alex nie model.matrix(asthma ~ gender + m_edu + p_edu + f_color + age + bmi_p)[, -1]dałby takiego samego wyniku jak dwie powyższe linie? Dlaczego warto skorzystać z dodatkowego kroku, aby połączyć zmienne ciągłe data.frame?
jiggunjer
6

Użyję pakietu enet, ponieważ jest to moja preferowana metoda. Jest trochę bardziej elastyczny.

install.packages('elasticnet')
library(elasticnet)

age <- c(4,8,7,12,6,9,10,14,7) 
gender <- c(1,0,1,1,1,0,1,0,0)
bmi_p <- c(0.86,0.45,0.99,0.84,0.85,0.67,0.91,0.29,0.88)
m_edu <- c(0,1,1,2,2,3,2,0,1)
p_edu <- c(0,2,2,2,2,3,2,0,0)
#f_color <- c("blue", "blue", "yellow", "red", "red", "yellow", "yellow", "red", "yellow")
f_color <- c(0, 0, 1, 2, 2, 1, 1, 2, 1)
asthma <- c(1,1,0,1,0,0,0,1,1)
pred <- cbind(age, gender, bmi_p, m_edu, p_edu, f_color)



enet(x=pred, y=asthma, lambda=0)
bdeonovic
źródło
4
dzięki za udostępnienie elasticnet; nie wiem jednak, co sądzić o wynikach powyższego Rskryptu. Czy możesz wyjaśnić? Z góry dziękuję!
Matt Reichenbach,
4

Aby rozwinąć doskonały przykład dostarczony przez pat. Pierwotny problem stanowiły zmienne porządkowe (m_edu, p_edu) z nieodłączną kolejnością między poziomami (0 <1 <2 <3). W pierwotnej odpowiedzi Pata myślę, że zostały one potraktowane jako nominalne zmienne kategorialne bez kolejności między nimi. Mogę się mylić, ale uważam, że te zmienne powinny być kodowane w taki sposób, aby model szanował ich naturalną kolejność. Jeśli są one zakodowane jako czynniki uporządkowane (a nie jako czynniki nieuporządkowane jak w odpowiedzi Pata), to glmnet daje nieco inne wyniki ... Myślę, że poniższy kod poprawnie zawiera zmienne porządkowe jako czynniki uporządkowane i daje nieco inne wyniki:

library(glmnet)

age     <- c(4, 8, 7, 12, 6, 9, 10, 14, 7) 
gender  <- as.factor(c(1, 0, 1, 1, 1, 0, 1, 0, 0))
bmi_p   <- c(0.86, 0.45, 0.99, 0.84, 0.85, 0.67, 0.91, 0.29, 0.88) 
m_edu   <- factor(c(0, 1, 1, 2, 2, 3, 2, 0, 1), 
                  ordered = TRUE)
p_edu   <- factor(c(0, 2, 2, 2, 2, 3, 2, 0, 0), 
                  levels = c(0, 1, 2, 3), 
                  ordered = TRUE)
f_color <- as.factor(c("blue", "blue", "yellow", "red", "red", 
                       "yellow", "yellow", "red", "yellow"))
asthma <- c(1, 1, 0, 1, 0, 0, 0, 1, 1)

xfactors <- model.matrix(asthma ~ gender + m_edu + p_edu + f_color)[, -1]
x        <- as.matrix(data.frame(age, bmi_p, xfactors))

# Note alpha=1 for lasso only and can blend with ridge penalty down to
# alpha=0 ridge only.
glmmod <- glmnet(x, y=as.factor(asthma), alpha=1, family="binomial")

# Plot variable coefficients vs. shrinkage parameter lambda.
plot(glmmod, xvar="lambda")

wprowadź opis zdjęcia tutaj

czasami_sci
źródło
1
czasami_sci, dobry haczyk - byłby to bardziej odpowiedni sposób na modelowanie zmiennych poziomu wykształcenia. Dziękuję za twój wkład.
Matt Reichenbach,
jak dodać legendę fabuły dla zmiennych? Np. Jaka jest czerwona linia w tym przykładzie?
jiggunjer