Ghostwriting-Service Dr. Rainer Hastedt

Erfolgsmessung im Marketing - Folge 28: Konversionsvergleiche für mehr als zwei Gruppen (logistische Regression)

Ich bin in Folge 27 für mein Zahlenbeispiel zu dem Ergebnis gekommen, dass sich die Konversionswahrscheinlichkeiten für mindestens zwei der vier Landing-Page-Varianten signifikant unterscheiden. Ich will jetzt untersuchen, für welche Varianten dies zutrifft.

Ein hierfür geeignetes Verfahren ist die logistische Regression, mit der ich die vier Landing-Page-Varianten anhand von Regressionskoeffizienten vergleichen kann.

Ich erläutere zunächst den Ansatz der logistischen Regression. Anschließend schätze ich mein logistisches Regressionsmodell und validiere es mit zwei Devianztests, einmal auf Basis der Chi-Quadrat-Verteilung und einmal als Permutationstest. Zum Schluss bespreche ich die Schätzergebnisse.

Logistische Regression (Grundidee)

Mit einer logistischen Regression schätzen Sie, bezogen auf mein Zahlenbeispiel, die folgende Gleichung:

log(p/(1-p)) = β0 + β1x1 + β2x2 + β3x3 + β4x4

Die unabhängigen Variablen x1 bis x4 stehen für die 140 Beobachtungswerte, jeweils entweder 1 für Konversion ja oder 0 für Konversion nein. x1 umfasst die 35 Beobachtungswerte für die Landing-Page-Variante H1T1, x2 die 35 Beobachtungswerte für H1T2, x3 die 35 Beobachtungswerte für H2T1 und x4 die 35 Beobachtungswerte für H2T2.

β0, β1, β2, β3, und β4 sind die zu schätzenden Regressionskoeffizienten.

Als Schätzverfahren für eine logistische Regression dient normalerweise die Maximum-Likelihood-Methode. Danach werden die Schätzwerte für β0, β1, β2, β3, und β4 so gewählt, dass die mit diesen Schätzwerten und den 140 Beobachtungswerten berechneten Konversionswahrscheinlichkeiten p = (p1, …, p140) den tatsächlichen Ergebnissen so nahe kommen wie möglich. Allen Beobachtungswerten gleich Eins (Konversion ja) soll ein möglichst hoher Schätzwert für die Konversionswahrscheinlichkeit entsprechen und allen Beobachtungswerten gleich Null (Konversion nein) ein möglichst niedriger.

Die anhand der Schätzwerte für die Regressionskoeffizienten berechneten Schätzwerte für die Konversionswahrscheinlichkeiten werden auch als gefittete Werte bezeichnet. Leads aus der gleichen Gruppe (Leads, denen beim Test die gleiche Variante der Landing-Page zugewiesen worden ist) haben den gleichen Regressionskoeffizienten und erhalten daher den gleichen gefitteten Wert. Ein gefitteter Wert ist daher ein Schätzwert für die Konversionswahrscheinlichkeit auf der betreffenden Landing-Page-Variante.

Die Transformation der abhängigen Variablen - log(p/(1-p)) statt p - verhindert unsinnige Schätzwerte (Wahrscheinlichkeiten kleiner als Null oder größer als Eins). Dieses Ziel ließe sich auch durch andere Transformationen erreichen.

Für die logistische Transformation log(p/(1-p)) spricht, dass die Regressionskoeffizienten β1, β2, β3, und β4 bei dieser Transformation Logits sind, was ich weiter unten näher erläutern werde. Eine andere Transformation wie zum Beispiel log(-log(1-p)) würde andere Regressionskoeffizienten liefern und wäre daher für meine Zwecke ungeeignet.

Für die logistische Transformation spricht außerdem, dass sich damit auch S-förmige Zusammenhänge schätzen lassen.

Betrachten Sie als Beispiel die folgende Regressionsgleichung:

log(p/(1-p)) = α + βx

Nehmen Sie an, ich hätte für die Parameter α und β die Schätzwerte a und b ermittelt. Auf dieser Basis kann ich die Gleichung

log(p/(1-p)) = a + bx

nach p auflösen und somit p in Abhängigkeit von x als Grafik darstellen.

Ich berechne p:

pwert <- function(a,b,x) {
 exp(a+b*x)/(1+exp(a+b*x))
 }

Ich zeichne p in Abhängigkeit von x für a = -12 und b = 0,2:

a <- -12
b <- 0.2
x <- seq(20,100,by=1)
y <- pwert(a,b,x)
dfr <- data.frame(x,y)
library(ggplot2)
svg("emim-28-1.svg", width=4.3, height=3.1, bg="transparent")
ggplot(data=dfr,aes(x=x,y=y)) + geom_line() + theme_bw() +
labs(title="b > 0, S-förmig",
x="Unabhängige Variable", y="Wahrscheinlichkeit") +
theme(plot.title=element_text(size=rel(1)),
axis.text=element_text(size=rel(0.7)),
axis.title=element_text(size=rel(0.9)),
plot.background=element_rect(colour=NA,
fill="transparent"))
dev.off()

Schätzung und Validierung des logistischen Regressionsmodells

Ich lese zunächst die für mein Zahlenbeispiel verwendeten Daten von meiner Website (kurze Erläuterungen dazu in Folge 27):

url <- "http://www.ghostwriting-service.de/erfolgsmessung-im-marketing-folge-24-konversionsvergleiche-fuer-mehr-als-zwei-gruppen-einfuehrung.php"
library(XML)
daten <- readHTMLTable(url, which=4, stringsAsFactors=FALSE)
daten$Konversion <- as.numeric(daten$Konversion)

Ich ergänze die Spalte mit der Gruppeneinteilung:

library(dplyr)
daten2 <- mutate(daten,
 Variante=paste(Headline, Beschreibung, sep=""))

Die ersten fünf der insgesamt 140 Zeilen (Kennzahl-Spalte weggelassen):

Headline	Beschreibung	Konversion	Variante
H1	T1	0	H1T1
H2	T1	1	H2T1
H1	T1	1	H1T1
H2	T1	1	H2T1
H1	T1	0	H1T1

Ich schätze jetzt meine Regressionsgleichung:

reg <- glm(Konversion~Variante, daten2, family=binomial(link="logit"))

Mit link="logit" ist die logistische Transformation log(p/(1-p)) gemeint. Ich habe demnach eine logistische Regression gewählt. Die von mir erwähnte Transformation log(-log(1-p)) ließe sich durch link="cloglog" wählen.

Die durch meine logistische Regression geschätzten 140 Konversionswahrscheinlichkeiten sind im Vektor reg$fitted gespeichert, für den ich eine Häufigkeitstabelle erstelle (auf vier Dezimalstellen gerundet):

table(round(reg$fitted,4))

1	2	3	4
0.4571	0.5143	0.5714	0.7714
35	35	35	35

Es gibt für meine Daten vier verschiedene Schätzwerte für die Konversionswahrscheinlichkeit, die jeweils 35mal auftreten. Ich habe vier Gruppen mit jeweils 35 Leads.

Ich berechne jetzt die vier Gruppenmittelwerte, ebenfalls auf vier Dezimalstellen gerundet:

summarise(group_by(daten2,Variante), MW=round(mean(Konversion),4))

Ich ergänze die auf diese Weise berechnete Tabelle um eine Spalte mit den Schätzwerten für die Konversionswahrscheinlichkeiten:

Variante	MW	reg$fitted
H1T1	0.7714	0.7714
H1T2	0.5143	0.5143
H2T1	0.5714	0.5714
H2T2	0.4571	0.4571

Die anhand der logistischen Regression geschätzten Konversionswahrscheinlichkeiten (Spalte reg$fitted) entsprechen demnach den anhand der Daten berechneten Gruppenmittelwerten (Spalte MW, durchschnittliche Konversionsraten).

Der nächste Schritt besteht darin, das geschätzte logistische Regressionsmodell auf seine Tauglichkeit zu untersuchen, was ich mit einem Devianztest machen will.

Bei einem Devianztest beschäftige ich mich mit der Frage, ob mein Regressionsmodell besser ist als das so genannte Nullmodell. In meinem Regressionsmodell habe ich vier Gruppen berücksichtigt und für jede dieser Gruppen einen Regressionskoeffizienten geschätzt. Im Nullmodell verschmelze ich die vier Gruppen zu einer einzigen und schätze für diese Gruppe einen einheitlichen Regressionskoeffizienten.

Ich beginne mit der als Devianz bezeichneten Kennzahl, die sich am einfachsten erläutern lässt, wenn ich vorher die Devianzfunktion zeichne.

Die Devianzfunktion hat als unabhängige Variable die innerhalb der Funktion mit x bezeichnete Konversionswahrscheinlichkeit:

devianz <- function(x) {
 -2*(x*log(x) + (1-x)*log(1-x))
 }

Ich bilde für diese Funktion eine Tabelle mit Wertepaaren und zeichne diese Tabelle in Form einer Linie:

p <- seq(0,1,by=0.001)
y <- devianz(p)
dfr <- data.frame(p,y)
library(ggplot2)
svg("emim-28-2.svg", width=4.3, height=3.1, bg="transparent")
ggplot(data=dfr,aes(x=p,y=y)) + geom_line() +
theme_bw() +
labs(title="Devianzfunktion",
x="Konversionswahrscheinlichkeit",
y="Devianz") +
theme(plot.title=element_text(size=rel(1)),
axis.text=element_text(size=rel(0.7)),
axis.title=element_text(size=rel(0.9)),
plot.background=element_rect(colour=NA,
fill="transparent"))
dev.off()

Bei einer Konversionswahrscheinlichkeit von 0,51 würde meine Prognose lauten, dass die Konversion erfolgreich ist. Die Gefahr, mit dieser Prognose falsch zu liegen wäre sehr groß, weil die Wahrscheinlichkeit für »Konversion nein« 0,49 beträgt.

Bei einer Konversionswahrscheinlichkeit von 0,01 würde meine Prognose lauten, dass die Konversion fehlschlägt. Die Gefahr, mit dieser Prognose falsch zu liegen wäre sehr gering.

Je niedriger die Devianz, desto höher die Wahrscheinlichkeit, anhand der betreffenden Konversionswahrscheinlichkeit eine korrekte Prognose zu formulieren (Konversion ja oder nein).

Die Devianz ist an der Stelle 0,5 am höchsten. Weil die Kurve im mittleren Bereich eine relativ geringe Steigung hat, macht es für die Höhe der Devianz nur einen relativ geringen Unterschied, ob die Konversionswahrscheinlichkeit gleich 0,5 oder 0,56 ist. Weiter weg von der Mitte, wo die Kurve wesentlich steiler ist, würde sich ein gleichgroßer Unterschied zwischen zwei Konversionswahrscheinlichkeiten viel stärker auf die Höhe der Devianz auswirken.

Eine niedrige Devianz steht für eindeutige Verhältnisse, die Sie nur haben, wenn die Konversionswahrscheinlichkeit entweder sehr klein oder sehr groß ist.

Sie erhalten die Devianz für das Regressionsmodell, indem Sie für alle 140 geschätzten Konversionswahrscheinlichkeiten die Devianz berechnen und diese Werte addieren:

sum(devianz(reg$fitted))

Das Ergebnis - 182,1862 - erhalten Sie auch folgendermaßen:

reg$deviance

Dieser Wert entspricht einem Durchschnitt von ca. 1,3 pro Konversionswahrscheinlichkeit (182,1862 geteilt durch 140). Ein Blick auf die Grafik zur Devianzfunktion zeigt, dass dies sehr viel ist.

Ich erstelle eine Tabelle, in der ich die 140 tatsächlichen Werte (Konversion ja oder nein) den geschätzten Konversionswahrscheinlichkeiten gegenüberstelle.

data.frame(daten2$Konversion,reg$fitted)

Die ersten fünf Zeilen:

daten2$Konversion	reg$fitted
0	0.7714286
1	0.5714286
1	0.7714286
1	0.5714286
0	0.7714286

Das Lead aus der ersten Zeile hat das Angebot auf der Landing-Page ausgeschlagen (Konversion nein = 0). Zu diesem Lead würde daher eine geschätzte Konversionswahrscheinlichkeit unterhalb von 0,5 passen. Der Schätzwert für die Konversionswahrscheinlichkeit liegt in Wirklichkeit oberhalb von 0,5 (0,7714286), was zu einem Lead mit dem Wert 1 in der linken Spalte passen würde.

Die erste Zeile der obigen Tabelle ist so gesehen ein Beispiel für eine Fehlklassifikation. Richtig klassifiziert sind die Leads in den Zeilen 2, 3 und 4 (jeweils Konversion = 1 und Konversionswahrscheinlichkeit > 0,5). Die fünfte Zeile zeigt eine weitere Fehlklassifikation.

Durch meine Regressionsgleichung sind 40 Prozent der insgesamt 140 Leads falsch klassifiziert.

mean(ifelse(reg$fitted<0.5,0,1) != daten2$Konversion)

Die hohe Devianz und der hohe Anteil der falsch klassifizierten Leads zeigen, dass meine logistische Regression von mäßiger Qualität ist.

Ich komme jetzt zum Devianztest und damit zu der Frage, ob es besser wäre, das Nullmodell zu verwenden.

In R geht dies wie folgt:

anova(reg, test="Chisq")

Das ausgegebene Signifikanzniveau in Höhe von 0,038 bedeutet, dass mein Regressionsmodell, bei dem ich zwischen den vier Varianten der Landing-Page differenziere laut Devianztest besser ist als das Nullmodell, bei dem ich die Unterscheidung zwischen den vier Varianten weglasse.

Der Devianztest setzt voraus, dass ich meine Daten als Stichprobe deuten kann und dass meine Daten zur Chi-Quadrat-Verteilung passen. Beides ist unzutreffend. Ich werde den Devianztest daher als Permutationstest wiederholen. Dies ist sachgemäß, weil die vier Gruppen in meinem Zahlenbeispiel durch Zufallsauswahl gebildet worden sind (das Zahlenbeispiel ist in Folge 24 erläutert, das Grundprinzip von Permutationstests in Folge 21 und Folge 22; Hinweise zu Chi-Quadrat-Tests finden sich in Folge 27).

Ich schätze zunächst das Nullmodell und damit die folgende Regressionsgleichung:

log(p/1-p) = β0

reg0 <- glm(Konversion~1, daten2, family=binomial(link="logit"))

Alle 140 Leads erhalten jetzt den gleichen Schätzwert für die Konversionswahrscheinlichkeit (0,5785714). Dieser Schätzwert entspricht der durchschnittlichen Konversionsrate aller 140 Leads.

table(reg0$fitted)
mean(daten2$Konversion)

Die Devianz des Nullmodells kann ich auf dreierlei Weise ermitteln:

reg0$deviance
reg$null.deviance
sum(devianz(reg0$fitted))

Als Devianz des Nullmodells erhalte ich 190,6097 und damit einen höheren Wert als für mein umfangreicheres Modell (182,1862).

Die Differenz zwischen den beiden Devianzwerten beträgt 8,42 (gerundet).

diff <- reg0$deviance-reg$deviance

Ich gehe von der Nullhypothese H0 aus, dass zwischen den vier Varianten der Landing-Page keine bedeutenden Unterschiede bestehen, dass es daher bedeutungslos ist, welche zufallsbedingte Gruppeneinteilung sich ergeben hat und dass mein umfangreicheres Regressionsmodell somit nicht besser ist als das Nullmodell.

Gesetzt den Fall, die Nullhypothese ist richtig. Wie groß ist dann die Wahrscheinlichkeit für eine Differenz zwischen den Devianzwerten größer oder gleich 8,42?

Ich beantworte die Frage mit dem folgenden Permutations-Devianztest:

nsim <- 20000
diffneu <- vector(length=nsim)
for (i in 1:nsim) {
 neueinteilung <- sample(daten2$Variante)
 regneu <- glm(daten2$Konversion~neueinteilung,
 family=binomial(link="logit"))
 diffneu[i] <- reg0$deviance-regneu$deviance
 }
bedingung <- diffneu >= diff mean(bedingung)

Wenn ich die Simulation mit nsim = 20.000 Iterationen mehrmals durchführe, dann erhalte ich im Allgemeinen unterschiedliche Werte, die jeweils bei 0,04 liegen. Ich lehne die Nullhypothese daher ab und verwende das erweiterte Modell.

Mein Permutations-Devianztest führt für meine Daten zur gleichen Entscheidung wie der Devianztest anhand der Chi-Quadrat-Verteilung. Dies muss nicht immer so sein. Der Devianztest anhand der Chi-Quadrat-Verteilung liefert für meine Daten ein deutlich geringeres Signifikanzniveau.

Die Schätzergebnisse der logistischen Regression

Als Ausgangspunkt dient meine Regressionsgleichung

log(p/(1-p)) = β0 + β1x1 + β2x2 + β3x3 + β4x4

x1 entspricht den Werten aus der Spalte Konversion für Variante H1T1, x2 den Werten aus der Spalte Konversion für H1T2, x3 den Werten aus der Spalte Konversion für H2T1 und x4 den Werten aus der Spalte Konversion für H2T2.

p ist der Vektor mit den Konversionswahrscheinlichkeiten, p/(1-p) der Vektor mit den Odds für eine erfolgreiche Konversion und log(p/(1-p)) der Vektor mit den logarithmierten Odds für eine erfolgreiche Konversion (Logits). Wahrscheinlichkeiten, Odds und Logits - p, p/(1-p) und log(p/(1-p)) - messen das Gleiche, nur auf unterschiedlichen Skalen.

Eine grafische Darstellung der Odds und Logits in Abhängigkeit von der Wahrscheinlichkeit verdeutlicht dies:

Odds <- function(x) {x/(1-x)}
x <- seq(0,0.85,by=0.001)
y <- Odds(x)
y2 <- log(y)
dfr <- data.frame(x,y,y2)
library(ggplot2)
svg("emim-28-3.svg", width=4.3, height=3.1, bg="transparent")
ggplot(data=dfr,aes(x=x,y=y)) + geom_line() +
geom_line(data=dfr,(aes(x=x,y=y2))) +
theme_bw() +
labs(title="Odds und Logits",
x="Wahrscheinlichkeit p",
y="Odds / Logits") +
annotate("text", x=0.1, y=1, label="Odds = p/(1-p)", size=4, hjust=0, vjust=0) +
annotate("text", x=0.3, y=-1.2, label="Logit = log(p/(1-p))", size=4, hjust=0, vjust=1) +
theme(plot.title=element_text(size=rel(1)),
axis.text=element_text(size=rel(0.7)),
axis.title=element_text(size=rel(0.9)),
plot.background=element_rect(colour=NA,
fill="transparent"))
dev.off()

Beide Linien sind streng monoton steigend. Ist Wahrscheinlichkeit p1 größer als Wahrscheinlichkeit p2, so gilt dies auch für die entsprechenden Odds und die entsprechenden Logits. Umgekehrt: Ist Logit l3 größer als Logit l4, so gilt dies auch für die entsprechenden Odds und die entsprechenden Wahrscheinlichkeiten.

Im letzten Abschnitt hatte ich erläutert, dass die durchschnittlichen Konversionsraten der vier Gruppen mit den Schätzwerten der vier Konversionswahrscheinlichkeiten übereinstimmen. Ich berechne jetzt die entsprechenden Odds und Logits.

summarise(group_by(daten2,Variante),
p=round(mean(Konversion),4),
Odds=round(mean(Konversion)/(1-mean(Konversion)),4),
Logit=round(log(Odds),4))

Variante	p	Odds	Logit
H1T1	0,7714	3,375	1,2164
H1T2	0,5143	1,0588	0,0571
H2T1	0,5714	1,3333	0,2877
H2T2	0,4571	0,8421	-0,1719

Die Werte in der Spalte p sind die geschätzten Konversionswahrscheinlichkeiten, die ich in der Spalte Odds als Odds angegeben habe und in der Spalte Logit als Logits (logarithmierte Odds).

Für meine Regressionsgleichung

log(p/(1-p)) = β0 + β1x1 + β2x2 + β3x3 + β4x4

sind die in der Spalte Logit angegebenen Werte die Schätzwerte für die Parameter β1, β2, β3, und β4. Der Schätzwert für die Niveaukonstante β0 ist Null.

Signifikanztests über die Parameter β1, β2, β3, und β4 oder Linearkombinationen dieser Parameter sind daher Signifikanztests über Konversions-Logits (logarithmierte Odds für eine erfolgreiche Konversion).

R liefert die Schätzwerte standardmäßig in der für meine Zwecke passenden Form:

round(coefficients(summary(reg))[,c(1,3,4)],digits=4)

	Estimate	z value	Pr(>|z|)
(Intercept)	1,2164	3,0218	0,0025
VarianteH1T2	-1,1592	-2,2049	0,0275
VarianteH2T1	-0,9287	-1,7592	0,0785
VarianteH2T2	-1,3882	-2,6369	0,0084

Die rechte Spalte der Tabelle enthält die Ergebnisse von z-Tests, die auf der Normalverteilungsannahme beruhen.

»(Intercept)« = 1,2164 ist der Schätzwert für den Parameter β1 und daher für die logarithmierten Odds einer erfolgreichen Konversion auf Variante H1T1 der Landing-Page. Die Nullhypothese des entsprechenden Signifikanztests lautet H0: β1 = 0.

H1T1 fungiert hier als Referenzkategorie. Dies ist sinnvoll, weil Variante H1T1 die höchste durchschnittliche Konversionsrate aufweist und ich prüfen möchte, ob diese Variante besser ist als die anderen Varianten.

»VarianteH1T2« = -1,1592 ist der Schätzwert für die Differenz β2 - β1 (H1T2 im Vergleich zu H1T1). Die Nullhypothese des entsprechenden Signifikanztests lautet H0: β2 - β1 = 0. Der z-Test legt nahe, H0 abzulehnen (Signifikanzniveau kleiner als 0,05).

»VarianteH2T1« = -0,9287 ist der Schätzwert für die Differenz β3 - β1 (H2T1 im Vergleich zu H1T1). Die Nullhypothese des entsprechenden Signifikanztests lautet H0: β3 - β1 = 0. Der z-Test legt nahe, H0 beizubehalten (Signifikanzniveau größer als 0,05).

»VarianteH2T2« = -1,3882 ist der Schätzwert für die Differenz β4 - β1 (H2T2 im Vergleich zu H1T1). Die Nullhypothese des entsprechenden Signifikanztests lautet H0: β4 - β1 = 0. Der z-Test legt nahe, H0 abzulehnen (Signifikanzniveau kleiner als 0,05).

Ich wiederhole die Signifikanztests als Permutationstests:

zwerte <- coefficients(summary(reg))[,"z value"]
nsim <- 20000
zwerteneu <- matrix(nrow=nsim, ncol=length(zwerte),byrow=TRUE)
colnames(zwerteneu) <- names(zwerte)
for (i in 1:nsim) {
 neueinteilung <- sample(daten2$Variante)
 regneu <- glm(daten2$Konversion~neueinteilung,
 family="binomial")
 zwerteneu[i,] <-
 coefficients(summary(regneu))[,3]
 }
zwerteneu
pwerte <- vector(length=4)
for (i in 1:4) {
 pwerte[i] <- mean(abs(zwerteneu[,i]) >=
 abs(zwerte[i]))
 }
pwerte

Ich vergleiche meine Ergebnisse:

Nullhypothese	z-Test	Permutations-z-Test
H0: β1 = 0	0,003	0,005
H0: β2 - β1 = 0	0,027	0,022
H0: β3 - β1 = 0	0,079	0,068
H0: β4 - β1 = 0	0,008	0,006

Meine Permutations-z-Tests liefern demnach im Wesentlichen das Gleiche wie die herkömmlichen z-Tests. Die Ergebnisse der Permutationstests variieren nur geringfügig, wenn ich die Simulation mehrmals wiederhole.

Somit unterscheiden sich die Logits für eine erfolgreiche Konversion zwischen Variante H1T1 und Variante H1T2 sowie zwischen Variante H1T1 und Variante H2T2 signifikant.

Es besteht demnach Kausalität:

Die im Vergleich zu H1T2 und H2T2 höhere durchschnittliche Konversionsrate für Variante H1T1 ist nicht zufallsbedingt, sondern systematisch, bedingt durch die gewählte Kombination aus Headline und Produktbeschreibung.

Meine Daten legen nahe, Produktbeschreibung T1 mit Headline H1 für eindeutig besser zu halten als Produktbeschreibung T2 mit Headline H1 oder H2.