Ein Beispiel: LASSO-Regression unter Verwendung von glmnet für binäre Ergebnisse

77

Ich beginne mit der Verwendung von dabble glmnetmit LASSO Regression , wo mein Ergebnis von Interesse dichotomous ist. Ich habe unten einen kleinen nachgebildeten Datenrahmen erstellt:

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)

Die Spalten (Variablen) im obigen Datensatz lauten wie folgt:

  • age (Alter des Kindes in Jahren) - kontinuierlich
  • gender - binär (1 = männlich; 0 = weiblich)
  • bmi_p (BMI-Perzentil) - kontinuierlich
  • m_edu (höchste Schulstufe der Mutter) - ordinär (0 = weniger als Gymnasium; 1 = Abitur; 2 = Bachelor-Abschluss; 3 = Abschluss nach dem Abitur)
  • p_edu (höchster Bildungsabschluss des Vaters) - ordinal (wie m_edu)
  • f_color (Lieblingsgrundfarbe) - nominal ("blau", "rot" oder "gelb")
  • asthma (Asthmastatus des Kindes) - binär (1 = Asthma; 0 = kein Asthma)

Das Ziel dieses Beispiels ist die Verwendung von LASSO zu machen , ein Modell der Vorhersage Kind Asthma - Status aus der Liste der sechs potentieller Einflussvariablen zu erstellen ( age, gender, bmi_p, m_edu, p_edu, und f_color). Natürlich ist die Stichprobengröße hier ein Problem, aber ich hoffe, mehr Einsicht in den Umgang mit den verschiedenen Arten von Variablen (dh kontinuierlich, ordinal, nominal und binär) innerhalb des glmnetFrameworks zu erhalten, wenn das Ergebnis binär ist (1 = Asthma) ; 0 = kein Asthma).

Wäre jemand bereit, ein Beispielskript Rzusammen mit Erklärungen für dieses Beispiel bereitzustellen , in dem LASSO mit den oben genannten Daten zur Vorhersage des Asthmastatus verwendet wird? Obwohl sehr einfach, ich weiß, ich und wahrscheinlich viele andere im Lebenslauf, würde dies sehr schätzen!

Matt Reichenbach
quelle
2
Sie erhalten möglicherweise mehr Glück, wenn Sie die Daten als dputein tatsächliches R-Objekt buchen. Lassen Sie die Leser nicht Zuckerguss darauf legen und backen Sie auch keinen Kuchen !. Wenn Sie den entsprechenden Datenrahmen in R generieren foo, bearbeiten Sie beispielsweise in der Frage die Ausgabe von dput(foo).
Gavin Simpson
Vielen Dank an GavinSimpson! Ich habe den Beitrag mit einem Datenrahmen aktualisiert, damit ich hoffentlich etwas Kuchen essen kann, ohne zu bereifen! :)
Matt Reichenbach
2
Durch die Verwendung des BMI-Perzentils widersetzen Sie sich gewissermaßen den Gesetzen der Physik. Fettleibigkeit betrifft Personen nach physischen Maßen (Länge, Volumen, Gewicht), nicht danach, wie viele Personen dem aktuellen Thema ähnlich sind, wie Perzentilierung dies tut.
Frank Harrell
3
Ich stimme zu, das BMI-Perzentil ist keine Metrik, die ich am liebsten verwende. In den CDC-Richtlinien wird jedoch empfohlen, für Kinder und Jugendliche unter 20 Jahren das BMI-Perzentil gegenüber dem BMI zu verwenden (ebenfalls eine sehr fragwürdige Messgröße!), da neben Größe und Gewicht auch Alter und Geschlecht berücksichtigt werden. Alle diese Variablen und Datenwerte wurden vollständig für dieses Beispiel ausgedacht. Dieses Beispiel spiegelt keine meiner aktuellen Arbeiten wider, da ich mit Big Data arbeite. Ich wollte nur ein Beispiel von glmnetin Aktion mit einem binären Ergebnis sehen.
Matt Reichenbach
Hier finden Sie ein Paket von Patrick Breheny mit dem Namen ncvreg, das für lineare und logistische Regressionsmodelle geeignet ist, die von MCP, SCAD oder LASSO bestraft werden. ( cran.r-project.org/web/packages/ncvreg/index.html )
bdeonovic

Antworten:

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")

Bildbeschreibung hier eingeben

Kategoriale Variablen werden normalerweise zuerst in Faktoren transformiert, dann wird eine Dummy-Variablenmatrix von Prädiktoren erstellt und zusammen mit den kontinuierlichen Prädiktoren an das Modell übergeben. Denken Sie daran, dass glmnet sowohl Ridge- als auch Lasso-Strafen verwendet, aber auch für sich alleine festgelegt werden kann.

Einige Ergebnisse:

# 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
# ...

Koeffizienten können aus dem glmmod extrahiert werden. Hier gezeigt mit 3 ausgewählten Variablen.

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 

Schließlich kann auch eine Kreuzvalidierung verwendet werden, um Lambda auszuwählen.

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

Bildbeschreibung hier eingeben

(best.lambda <- cv.glmmod$lambda.min)
# [1] 0.2732972
klopfen
quelle
4
Dies ist genau das, wonach ich gesucht habe +1, die einzigen Fragen, die ich habe, sind 1) Was können Sie mit der Kreuzvalidierung Lambda von 0,2732972 tun? und 2) Sind die ausgewählten Variablen aus dem Standard die Lieblingsfarbe (gelb), das Geschlecht und die Ausbildung des Vaters (Bachelor-Abschluss)? Vielen Dank!
Matt Reichenbach
4
1) Mit der Kreuzvalidierung werden Lambda und Koeffizienten (bei minimalem Fehler) ausgewählt. In diesem Modell gibt es kein lokales min (es gab eine Warnung, die auch auf zu wenige obs bezogen wurde); Ich würde interpretieren, dass alle Koeffizienten mit den Schrumpfungsstrafen auf Null geschrumpft wurden (das beste Modell hat nur einen Schnitt) und mit mehr (echten) Beobachtungen wiederholt wurden und möglicherweise den Lambda-Bereich vergrößerten. 2) Ja, in dem Beispiel, in dem ich coef (glmmod) [, 10] gewählt habe, wählen Sie Lambda für das Modell über den Lebenslauf oder die Interpretation der Ergebnisse. Könnten Sie als gelöst markieren, wenn Sie der Meinung wären, dass ich Ihre Frage gelöst habe? Vielen Dank.
Pat
2
Kann ich fragen, wie dies mit der f_colorVariablen umgeht? Wird Faktorstufe 1 bis 4 als eine größere Stufe angesehen als 1 bis 2, oder sind alle gleich gewichtet, ungerichtet und kategorisch? (Ich möchte es auf eine Analyse mit allen ungeordneten Prädiktoren anwenden.)
beroe
3
Die Zeile xfactors <- model.matrix(asthma ~ gender + m_edu + p_edu + f_color)[,-1]codiert die kategoriale Variable f_color (wie as.factorin den vorherigen Zeilen deklariert ). Es sollte die standardmäßige R-Dummy-Variablencodierung verwendet werden, sofern das contrasts.argArgument nicht angegeben wird. Dies bedeutet, dass alle Ebenen von f_color gleich gewichtet und ungerichtet sind, mit Ausnahme der ersten, die als Referenzklasse verwendet und in den Achsenabschnitt aufgenommen wird.
Alex
1
@Alex würde nicht model.matrix(asthma ~ gender + m_edu + p_edu + f_color + age + bmi_p)[, -1]dasselbe Ergebnis liefern wie die beiden obigen Zeilen? Warum einen zusätzlichen Schritt verwenden, um die stetigen Variablen mit zu verketten data.frame?
Jiggunjer
6

Ich werde das Paket enet verwenden, da dies meine bevorzugte Methode ist. Es ist etwas flexibler.

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
quelle
4
danke fürs teilen elasticnet; Ich weiß jedoch nicht, was ich von der Ausgabe des obigen RSkripts halten soll. Können Sie bitte klarstellen? Danke im Voraus!
Matt Reichenbach
4

Nur um das hervorragende Beispiel von pat zu erweitern. Das ursprüngliche Problem stellte Ordinalvariablen (m_edu, p_edu) mit einer inhärenten Reihenfolge zwischen Ebenen (0 <1 <2 <3). In der ursprünglichen Antwort von pat denke ich, dass diese als nominale kategoriale Variablen ohne Reihenfolge zwischen ihnen behandelt wurden. Ich kann mich irren, aber ich glaube, diese Variablen sollten so codiert werden, dass das Modell ihre inhärente Reihenfolge respektiert. Wenn diese als geordnete Faktoren codiert sind (anstatt als ungeordnete Faktoren wie in Pat's Antwort), dann liefert glmnet leicht unterschiedliche Ergebnisse ... Ich denke, der folgende Code enthält die Ordinalvariablen korrekt als geordnete Faktoren und es liefert leicht unterschiedliche Ergebnisse:

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")

Bildbeschreibung hier eingeben

times_sci
quelle
1
times_sci, good catch - dies wäre der geeignetere Weg, die Variablen des Bildungsniveaus zu modellieren. Danke für Ihren Beitrag.
Matt Reichenbach
Wie würde man eine Plotlegende für die Variablen hinzufügen? Was ist zB die rote Linie in diesem Beispiel?
Jiggunjer