Ein Arbeitsblatt zu statistischen Testverfahren

Wir verwenden das frei erhältliche Statistik-Programm RStudio. Dieses HTML-Dokument kann im RMD-Format abgespeichert werden, siehe die Schaltfläche Code rechts oben. Das RMD-Dokument kann dann in RStudio geöffnet werden.

Binomialtests

Binomialtests werden verwendet, um Hypothesen über die unbekannte Wahrscheinlichkeit \(p\) einer Münze zu testen. Hierbei wirft man die Münze \(n\)-mal unabhängig und betrachtet die Zahl \(S_n\) der Erfolge. Ein kleiner Wert von \(S_n\) deutet auf ein kleines \(p\), ein großer Wert von \(S_n\) deutet auf ein großes \(p\). Die Verteilung der Zufallsgröße \(S_n\) ist die Binomialverteilung \(B_{n,p}\).

Rechtsseitiger Binomialtest

\[ H_0: p\le p_0 \qquad H_1:p>p_0 \]

Wir akzeptieren \(H_0\), falls \(S_n\) kleine Werte annimmt. Wir lehnen \(H_0\) ab, falls \(S_n\) große Werte annimmt. Wir verwerfen \(H_0\), wenn \(S_n\) einen kritischen Wert \(c_\alpha\) überschreitet. Dabei soll die Wahrscheinlichkeit für einen Fehler 1.Art höchstens \(\alpha\) betragen. Wir fordern also, dass für \(p\in H_0\) stets

\[ \mathbb P_p(H_0 \text{ verwerfen})=\mathbb P_p(S_n>c_\alpha) \le \mathbb P_{p_0}(S_n>c_\alpha) \le \alpha \]

gilt. In der ersten Ungleichung haben wir die Verwerfungs-Wahrscheinlichkeiten für \(p\in H_0\) durch ihren größten Wert ersetzt. Dies ist intuitiv verständlich: obiges Ereignis betrifft große Werte von \(S_n\). Die Wahrscheinlichkeit für große Werte von \(S_n\) steigt, falls \(p\) steigt.

Darüber hinaus soll die Güte des Tests hoch sein. Das bedeutet, dass wir \(c_\alpha\) möglichst klein wählen wollen, denn dann ist der Verwerfungsbereich möglichst groß. Wir verlangen also

\[ c_\alpha=\min\{k\in \mathbb N_0: \mathbb P_{p_0}(S_n>k)\le \alpha\} \]

Diese Zahl ist das kleinste \(\alpha\)-Fraktil bzw. \((1-\alpha)\)-Quantil der Binomialverteilung \(B_{n,p_0}\) mit Parametern \(n\) und \(p_0\).

rechtsseitiger Binomialtest in R

In R kann der rechtsseitige Binomialtest mit dem Befehl binom.test() durchgeführt werden. Wie zu allen Befehlen von R erhält man mit ?binom.test Informationen zum Befehl, und mit binom.test wird der Programmcode des Befehls angezeigt.

?binom.test

Wir testen auf \(p>1/2\) zum Niveau \(0.05\) mit den Daten der Jungengeburten in Erlangen im Jahr 2022.

k=536; n=1059 # Jungengeburten in Erlangen 2022
p_0=0.5; alpha=0.05
binom.test(k,n,p_0, alt="greater", conf.level=1-alpha)

    Exact binomial test

data:  k and n
number of successes = 536, number of trials = 1059, p-value = 0.3562
alternative hypothesis: true probability of success is greater than 0.5
95 percent confidence interval:
 0.4804054 1.0000000
sample estimates:
probability of success 
             0.5061379 

Der \(p\)-Wert der Zahl der Erfolge \(k\) ist definiert als \(p(k)=\mathbb P_{p_0}(S_n\ge k)\). Dies ist die Wahrscheinlichkeit, dass die Teststatistik einen Wert annimmt, der mindestens so extrem ist wie der beobachtete Wert. Der \(p\)-Wert einer Beobachtung ist also eine Maßzahl für die Plausibilität der Nullhypothese. Tatsächlich wird \(H_0\) abgelehnt genau dann wenn \(p(k)\le \alpha\) gilt. Das Konfidenzintervall gibt einen Bereich an, in dem der wahre unbekannte Wert von \(p\) mit Sicherheit \(1-\alpha\) erwartet werden kann.

Als zweites Beispiel führen wir den Test auf die Störfall-Wahrscheinlichkeit für ein AKW durch.

k=2; n=442*30 # 2 AKW-GAUs in 442*30 Reaktorjahren
p_0=4e-6 # berechnete Stoerfall-Wkeit
# p_0=10*4e-6 # Wkeit fuer Fragestellung "bedeutend groesser"
alpha=0.05
binom.test(k,n,p_0, alt="greater", conf.level=1-alpha)

    Exact binomial test

data:  k and n
number of successes = 2, number of trials = 13260, p-value = 0.001358
alternative hypothesis: true probability of success is greater than 4e-06
95 percent confidence interval:
 2.680016e-05 1.000000e+00
sample estimates:
probability of success 
          0.0001508296 

Gütefunktion des rechtsseitigen Binomialtests

Die Gütefunktion des rechtsseitigen Binomialtests kann man in R wie folgt berechnen.

n=23; p_0=0.7; alpha=0.05
# Definition der Gütefunktion
G=function(p) {
  # bestimme kleinstes (1-alpha)-Quantil von B(n,p0)
  # mit R-Befehl qbinom
  c_al=qbinom(1-alpha, n,p_0) 
  # alternativ Berechnung von Hand
  # c_al=min(which(pbinom(0:n,n,p_0)>=1-alpha))-1
  1-pbinom(c_al, n,p) # P_p(H0 ablehnen)
}
# Zeichne Gütefunktion im Intervall [0,1] und beschrifte die Achsen 
plot(G, main=paste("Gütefunktion rechtsseitiger Binomialtest, n=",n), xlab="p", ylab="G(p)", xlim=c(0,1), lwd=2)
# Trage Zusatzinformationen ein
abline(h=alpha, col="blue"); abline(v=p_0, col="red")

text(p_0, 0.7, pos=2, expression(p[0]), col="red")
text(0, alpha, pos=3, expression(alpha), col="blue")

Falls man die Stichprobengröße erhöht, nähert sich die stetige Gütekurve einer Sprungfunktion an, die bei \(p_0\) von \(0\) auf \(1\) springt. Dies kann man so interpretieren, dass der Test um so genauer wird, je mehr Daten zur Verfügung stehen. Tatsächlich kann man so für beliebiges \(p\in H_1\) eine Stichprobengröße bestimmen, so dass die Wahrscheinlichkeit für einen Fehler 2.Art bei \(p\) kleiner als ein vorgegebener Wert \(\beta\) ist. Damit hat man eine explizite Kontrolle über die Qualität des Tests.

effektives Niveau des rechtsseitigen Binomialtests

Da die Teststatistik diskret ist, kann man nicht erwarten, dass die Wahrscheinlichkeit für einen Fehler 1.Art ihren maximal zulässigen Wert \(\alpha\) erreicht. Der Test ist also konservativer als gefordert. Dies führt auch dazu, dass die Güte in der Nähe von \(p_0\) geringer ausfällt als vielleicht möglich.

p_0=0.7; alpha=0.05
Gn=function(n) {
  c_al=qbinom(1-alpha, n,p_0) # kleinstes (1-alpha)-Quantil von B(n,p0)
  1-pbinom(c_al, n,p_0) # P_p(H0 ablehnen)
}
# fuer plot-Befehl muss Funktion "vektorisiert" werden
Gnvec=Vectorize(Gn)
xn=10:150; yn=Gnvec(xn)
plot(xn,yn, type="h", main=paste("rechtsseitiger Binomialtest zum Niveau", alpha), xlab="Stichprobengröße n", ylab="effektives Niveau", ylim=c(0,2*alpha))
abline(h=alpha, col="blue")

text(xn[1],alpha, pos=3, expression(alpha), col="blue")

Dieses unerwünschte Phänomen tritt bei kontinuierlichen Zufallsgrößen nicht auf. Um das Phänomen zu umgehen, kann man bei diskreten Zufallsgrößen sogenannte randomisierte Tests betrachten.

randomisierter rechtsseitiger Binomialtest

Bei diesem Test entscheidet man per Münzwurf, ob schon bei \(c_\alpha\) die Nullhypothese verworfen wird. Die Wahrscheinlichkeit \(\gamma\) dieser Münze ist dabei so eingestellt, dass der randomisierte Test genau Wahrscheinlichkeit \(\alpha\) für einen Fehler 1. Art besitzt. Dies führt auf

\[ \gamma=\frac{\alpha-\mathbb P_{p_0}(S_n>c_\alpha)}{\mathbb P_{p_0}(S_n=c_\alpha)} \]

Wir vergleichen die Gütefunktionen beider Tests.

Gran=function(p){
  # berechne kritischen Wert
  c_al=qbinom(1-alpha, n,p_0) # kleinstes (1-alpha)-Quantil von B(n,p0)
  # berechne gamma
  gam=(alpha-G(p_0))/dbinom(c_al,n,p_0)
  # definiere Wkeit H0 abzulehnen
  G(p)+gam*dbinom(c_al,n,p)
}
plot(G, main=paste("Gütefunktion rechtsseitiger Binomialtest, n=",n), xlab="p", ylab="G(p)", xlim=c(0,1), lwd=2)
abline(h=alpha, col="blue"); abline(v=p_0, col="red")

text(p_0, 0.7, pos=2, expression(p[0]), col="red")
text(0, alpha, pos=3, expression(alpha), col="blue")

text(0,0.9, pos=4, "randomisiert", col="brown")
text(0,0.8, pos=4, "nicht randomisiert", col="black")

curve(Gran(x), add=TRUE, col="brown", lwd=2)

Offenbar hat der randomisierte Test gleichmäßig in \(p\) eine höhere Güte als der obige nicht randomisierte Test.

Tatsächlich ist dieser randomisierte Test innerhalb der Klasse aller randomisierten Tests für das obige Testproblem optimal. Dies bedeutet, dass es keinen anderen Niveau-\(\alpha\)-Test für das obige Problem gibt, welche eine gleichmäßig höhere Güte besitzt. Tatsächlich gilt sogar, dass jeder gleichmäßig beste Niveau-\(\alpha\)-Test die obige Form besitzt.

Linksseitiger Binomialtest

Das Hypothesenpaar für den linksseitigen Binomialtest ist

\[ H_0: p\ge p_0 \qquad H_1:p<p_0 \]

Wir verwerfen \(H_0\), falls \(S_n\) einen kritischen Wert \(c_\alpha\) unterschreitet. Dabei soll die Wahrscheinlichkeit für einen Fehler 1.Art höchstens \(\alpha\) betragen. Wir fordern also

\[ \mathbb P_p(H_0 \text{ verwerfen})=\mathbb P_p(S_n<c_\alpha) \le \mathbb P_{p_0}(S_n<c_\alpha) \le \alpha \]Darüber hinaus soll die Güte des Tests hoch sein. Das bedeutet, dass wir \(c_\alpha\) möglichst groß wählen. Wir haben also

\[ c_\alpha=\max\{k\in \mathbb N_0: \mathbb P_{p_0}(S_n<k)\le \alpha\} \]

Diese Zahl ist das größte \(\alpha\)-Quantil der Binomialverteilung \(B_{n,p_0}\) mit Parametern \(n\) und \(p_0\). Die Gütefunktion des linksseitigen Binomialtests kann man in R wie folgt bestimmen.

n=23; p_0=0.7; alpha=0.05
G=function(p) {
  # berechne groesstes alpha-Quantil von B(n,p0)
  # hierfuer gibt es keinen R-Befehl
  c_al=max(which(pbinom((-1):(n-1),n,p_0)<=alpha))-1
  pbinom(c_al-1, n,p) # P_p(H0 ablehnen)
}
plot(G, main=paste("Gütefunktion linksseitiger Binomialtest, n=",n), xlab="p", ylab="G(p)", xlim=c(0,1), lwd=2)
abline(h=alpha, col="blue"); abline(v=p_0, col="red")

text(p_0, 0.7, pos=2, expression(p[0]), col="red")
text(0, alpha, pos=3, expression(alpha), col="blue")

Zweiseitige Binomialtests

Wir konstruieren zweiseitige Binomialtests der Form

\[ H_0:p=p_0 \qquad H_1:p\ne p_0 \]

Konstruktion mit dem TOST-Verfahren

Beim TOST-Verfahren (Two One-Sided Tests) testen wir \(H_0^{(1)}:p\ge p_0\) gegen \(H_1^{(1)}: p<p_0\) und dann \(H_0^{(2)}:p\le p_0\) gegen \(H_1^{(2)}: p>p_0\), jeweils zum Niveau \(\alpha/2\). Wir lehnen \(H_0\) ab genau dann wenn \(H_0^{(1)}\) abgelehnt werden kann oder wenn \(H_0^{(2)}\) abgelehnt werden kann. Dies ergibt einen Niveau-\(\alpha\)-Test für das obige zweiseitige Testproblem.

Die Gütefunktion dieses Tests kann aus den obigen Tests wie folgt erhalten werden.

n=16; p_0=0.7; alpha=0.05
# c1 groesstes alpha/2-Quantil von B(n,p0)
# c2 kleinstes (1-alpha/2)-Quantil von B(n,p0)
# es gilt nach Konstruktion c1<=c2
c1_t=max(which(pbinom((-1):(n-1),n,p_0)<=alpha/2))-1
c2_t=min(which(pbinom(0:n,n,p_0)>=1-alpha/2))-1
# lehne H_0: p=p_0 ab, falls S_n<c1 oder S_n>c2
G_tost=function(p) pbinom(c1_t-1, n,p) + (1-pbinom(c2_t,n,p))
plot(G_tost, main=paste("Gütefunktion zweiseitiger TOST-Binomialtest, n=",n), xlab="p", ylab="G(p)", xlim=c(0,1), lwd=2)
abline(h=alpha, col="blue"); abline(v=p_0, col="red")

text(p_0, 0.7, pos=2, expression(p[0]), col="red")
text(0, alpha, pos=3, expression(alpha), col="blue")

Konstruktion mit dem MPVT-Verfahren

Beim MPVT-Verfahren (Most Probable Values Test) wählen wir als Akzeptanzbereich für die Nullhypothese diejenigen Werte aus, die für \(H_0\) am wahrscheinlichsten sind.

Dies kann man in R wie folgt implementieren.

n=16; p_0=0.7; alpha=0.05
# bestimme kritischen Wert
# Wkeiten der Werte unter H0
probs=dbinom(0:n,n,p_0)
# nach Groesse sortieren
cvals=sort(dbinom(0:n,n,p_0))
# fuege einen grosseren Wert ein
cvals=append(cvals,1.1)
# Funktion: Wert k klein genug fuer Niveau alpha?
adm=function(k) sum(probs*(probs<k))<=alpha
admvec=Vectorize(adm)
# waehle maximalen zulaessigen Wert in cvals
c_al=cvals[max(which(sapply(cvals, admvec)))]
# Definition der Guetefunktion
# lehne H_0: p=p_0 ab, falls Wkeit fuer Wert kleiner als c_al ist
G_m=function(p){
  probsp=dbinom(0:n,n,p)
  sum(probsp*(probs<c_al))
}
G_mvec=Vectorize(G_m)
# zeichne die Guetefunktion
plot(G_mvec, main=paste("Gütefunktion zweiseitiger MPVW-Binomialtest, n=",n), xlab="p", ylab="G(p)", xlim=c(0,1), ylim=c(0,1),lwd=2)
abline(h=alpha, col="blue"); abline(v=p_0, col="red")

text(p_0, 0.7, pos=2, expression(p[0]), col="red")
text(0, alpha, pos=3, expression(alpha), col="blue")

Konstruktion des optimalen zweiseitigen Tests

Die obigen Tests haben zwei kritische Werte \(c_1<c_2\). Wir lehnen \(H_0\) ab, falls \(S_n<c_1\) oder \(S_n>c_2\) gilt. Wir können solche Tests randomisieren, indem wir im Fall \(S_n=c_1\) die Nullhypothese \(H_0\) mit Wahrscheinlichkeit \(\gamma_1\) ablehnen und im Fall \(S_n=c_2\) die Nullhypothese \(H_0\) mit Wahrscheinlichkeit \(\gamma_2\) ablehnen. Solche Tests kann man mit einer Funktion \(\phi:\{0, \ldots, n\}\to [0,1]\) beschreiben, die jeder möglichen Beobachtung die Wahrscheinlichkeit zuordnet, mit der man die Gegenhypothese akzeptiert.

Der optimale Test für das zweiseitige Testproblem wird im Buch von Lehmann und Romano im Beispiel 4.2.1 beschrieben. Er ist von der obigen Form. Der Test werde durch die Funktion \(\phi\) beschrieben. Die Konstanten \(c_1,c_2, \gamma_1,\gamma_2\) sind eindeutig bestimmt über die zwei Gleichungen

\[ \mathbb E_{p_0}[\phi]=\alpha \ , \qquad \mathbb E_{p_0}[\phi \cdot S_n]=\alpha\cdot n p_0 \ . \]

Die erste Gleichung fordert, dass der Test das Niveau \(\alpha\) ausschöpft. Es soll also \(G_\phi(p_0)=\alpha\) gelten. Die zweite Gleichung legt den Erwartungswert von \(S_n\) auf dem Ablehnungsbereich fest. Es sollte sich der natürliche Werte \(\alpha\cdot n p_0\) ergeben.

Für festes \(c_1,c_2\) sind die obigen Gleichungen lineare Gleichungen für \(\gamma_1,\gamma_2\). Es gilt nämlich

\[ \begin{split} \gamma_1 \mathbb P_{p_0}[S_n=c_1]+\gamma_2 \mathbb P_{p_0}[S_n=c_2]&=\alpha-\mathbb P_{p_0}[S_n<c_1 \text{ oder } S_n>c_2]\\ \gamma_1 c_1\mathbb P_{p_0}[S_n=c_1]+\gamma_2 c_2\mathbb P_{p_0}[S_n=c_2]&= \alpha n p_0-\mathbb E_{p_0}[S_n \cdot 1_{\{S_n<c_1 \text{ oder } S_n>c_2\}}] \end{split} \]

Da die Konstanten \(0\le c_1<c_2\le n\) und \(0\le \gamma_1,\gamma_2<1\) durch obige Gleichungen eindeutig bestimmt sind, kann man sie durch “Ausprobieren” erhalten. Hierzu dient das folgende Programm. Es kann für große Stichproben noch optimiert werden.

n=16; p_0=0.7; alpha=0.05
# probiere alle (!) Wahlen c1<c2
# einfach zu programmieren
# aber fuer hohes n sehr ineffizient
for (c1 in 0:n){
  for (c2 in (c1+1):n){
    # Wkeit auf Ablehnungsbereich
    p=(pbinom(c1-1,n,p_0) + 1-pbinom(c2,n,p_0))
    # Erwartungswert auf Ablehnungsbereich
    probs=dbinom(0:n, n, p_0)
    e=sum((0:n)*probs*(((0:n)<c1)+((0:n)>c2)))
    if ((p>alpha) || (e>alpha*n*p_0)) next
    # bestimme Wkeiten fuer Muenzwurf bei kritischen Werten
    p1=dbinom(c1,n,p_0); p2=dbinom(c2,n,p_0)
    # loese die Matrixgleichung A*gam=b
    b=c(alpha-p, alpha*n*p_0-e)
    A=matrix(c(p1,c1*p1,p2,c2*p2), nrow = 2)
    # ignoriere zu große Werte gam_1, gam_2
    # diese erkennt man an sehr kleiner Determinante
    prec=.Machine$double.eps
    if (abs(det(A))<10*prec) next
    gam=solve(A,b)
    # Falls gam_1, gam_2 Wkeiten sind, haben wir die
    # korrekten Parameter gefunden
    # Ausgabe der korrekten Werte
    if (sum(gam>=0 & gam <1)==2) {
      print(c(c1,c2))
      print(gam)
    }
  }
}

[1]  7 14
[1] 0.85148648 0.01322789

Mit den oben bestimmten Parametern des optimalen Tests kann die Gütefunktion gezeichnet werden.

# lehne H_0: p=p_0 ab, falls S_n<c1 oder S_n>c2
# randomisiere Test mit Wkeiten gamma1, gamma2
# die folgenden Parameter wurden mit dem obigen Programm bestimmt
n=16; c1=7; c2=14; g1=0.85148648; g2=0.01322789
# Definition der Guetefunktion
G_opt=function(p) pbinom(c1-1, n,p) + (1-pbinom(c2,n,p)) + g1*dbinom(c1,n,p)+g2*dbinom(c2,n,p)
# Guetefunktion zeichnen
plot(G_opt, main=paste("Gütefunktion zweiseitiger Binomialtest, n=",n), xlab="p", ylab="G(p)", xlim=c(0,1), lwd=2, col="gray")
abline(h=alpha, col="blue"); abline(v=p_0, col="red")

text(p_0, 0.7, pos=2, expression(p[0]), col="red")
text(0, alpha, pos=3, expression(alpha), col="blue")

text(0,0.7, pos=4, "OPTI", col="gray")

Vergleich der Gütefunktionen

Wir vergleichen die Gütefunktion der drei Tests. Hierzu verwenden wir die nicht randomisierte Variante des optimalen Tests.

# lehne H_0: p=p_0 ab, falls S_n<c1 oder S_n>c2
c1=7; c2=14; g1=0.85148648; g2=0.01322789
G_opt=function(p) pbinom(c1-1, n,p) + (1-pbinom(c2,n,p)) 
#+ g1*dbinom(c1,n,p)+g2*dbinom(c2,n,p)
# (Zusatzterm fuer randomisierte Version)
# zeichne die drei nicht randomisierten Guetefunktionen
plot(G_mvec, main=paste("Gütefunktion zweiseitiger Binomialtest, n=",n), xlab="p", ylab="G(p)", xlim=c(0,1), lwd=2, col="black")
abline(h=alpha, col="blue"); abline(v=p_0, col="red")

text(p_0, 0.7, pos=2, expression(p[0]), col="red")
text(0, alpha, pos=3, expression(alpha), col="blue")

text(0,0.9, pos=4, "TOST", col="brown")
text(0,0.8, pos=4, "MPVT", col="black")

text(0,0.7, pos=4, "OPTI", col="gray")
curve(G_tost(x), add=TRUE, col="brown", lwd=2)

curve(G_opt(x), add=TRUE, col="gray", lwd=2)

Es ist auch instruktiv, die Gütefunktion des randomisierten optimalen Tests zu zeichnen. Dann kann man erkennen, dass der MPVT-Test teilweise eine höhere Güte als der optimale Test haben kann. Dies ist kein Widerspruch, weil die Optimalitätsaussage nur für unverfälschte Tests gilt. Der MPVT-Test ist aber verfälscht.

zweiseitiger Binomialtest in R

Der R-Befehl binom.test() erlaubt auch zweiseitige Binomialtests.

n=510400; k=334 # Geburtstagkinder 29. Februar in Nürnberg
p_0=(4*365+1)^(-1); alpha=0.05
binom.test(k,n,p_0, alt="two.sided", conf.level=1-alpha)

    Exact binomial test

data:  k and n
number of successes = 334, number of trials = 510400, p-value = 0.4376
alternative hypothesis: true probability of success is not equal to 0.0006844627
95 percent confidence interval:
 0.0005861023 0.0007284424
sample estimates:
probability of success 
          0.0006543887 

In binom.test() sind das TOST-Verfahren und das MPVT-Verfahren implementiert. Das TOST-Verfahren erhält man, in dem man das Konfidenzintervall für die Testentscheidung verwendet. Man akzeptiert \(H_0\), falls \(p_0\) im Konfidenzintervall für \(p\) liegt. Das MPVT-Verfahren erhält man, in dem für die Testentscheidung den \(p\)-Wert verwendet. Man lehnt \(H_0\) ab, falls \(p(k)\le \alpha\) gilt.

In unserem Beispiel führen beide Tests auf dieselbe Entscheidung. Dies ist auch der Fall für den optimalen Test, welcher in binom.test() nicht implementiert ist. Tatsächlich könnte man in unserem Beispiel auch einen \(t\)-Test oder einen \(z\)-Test verwenden, weil wegen der hohen Stichprobengröße die Normalapproximation eine gute Näherung ist (beachte \(n p_0(1-p_0)>10\)). Dies Verfahren ist weiter unten am Beispiel eines Äquivalenztests beschrieben.

Äquivalenztests

Bei Äquivalenztests untersucht man das Hypothesenpaar.

\[ H_0 : p\le p_1 \text{ oder } p\ge p_2 \qquad H_1:p_1<p<p_2 \]

Wir besprechen, wie solche Tests durchgeführt werden können. Es ist nicht schwer, hieraus Tests auf Inäquivalenz zu entwickeln.

TOST-Äquivalenztest

Beim TOST-Verfahren (Two One-Sided Tests) testen wir \(H_0^{(1)}:p\le p_1\) gegen \(H_1^{(1)}: p>p_1\) und dann \(H_0^{(2)}:p\ge p_2\) gegen \(H_1^{(2)}: p<p_2\), jeweils zum Niveau \(\alpha\). Wir lehnen \(H_0\) ab genau dann wenn \(H_0^{(1)}\) abgelehnt werden kann und wenn \(H_0^{(2)}\) abgelehnt werden kann. Dies ergibt einen Niveau-\(\alpha\)-Test für das obige zweiseitige Testproblem. Dieser Test ist asymptotisch optimal.

Der Test wird ähnlich wie beim TOST-Verfahren für den zweiseitigen Binomialtest konstruiert. Wir lehnen \(H_0\) ab, falls \(c_1<S_n<c_2\) gilt. Wir berechnen die kritischen Grenzen \(c_1\) und \(c_2\).

n=100; p_1=0.6; p_2=0.8; alpha=0.05
# c1 kleinstes (1-alpha)-Quantil von B(n,p_1)
# c2 groesstes alpha-Quantil von B(n,p_2)
c_1=min(which(pbinom(0:n,n,p_1)>=1-alpha))-1
c_2=max(which(pbinom((-1):(n-1),n,p_2)<=alpha))-1
c(c_1,c_2)

[1] 68 73

Falls für die kritischen Werte \(c_1<c_2\) gilt, ist der TOST-Test anwendbar. Dies ist für hinreichend große Stichproben der Fall.

n=100; p_1=0.6; p_2=0.8; alpha=0.05
# bestimme kritische Werte wie oben
c_1=68; c_2=73
# lehne H0 ab, falls c_1<S_n<c_2 gilt
G_tost=function(p) {
  vals=dbinom(0:n, n, p)
  sum(vals*(((0:n)>c_1) & ((0:n)<c_2)))
}
G_tostv=Vectorize(G_tost)
# zeichne Guetefunktion
plot(G_tostv, main=paste("Gütefunktion TOST-Äquivalenztest, n=",n), xlab="p", ylab="G(p)", xlim=c(0,1), lwd=2)
# zeichne Niveau-alpha-Linie
abline(h=alpha, col="blue")

text(0, alpha, pos=3, expression(alpha), col="blue")
# zeichne p1, p2
abline(v=p_1, col="red"); abline(v=p_2, col="red")

text(p_1, 0.7, pos=2, expression(p[1]), col="red")
text(p_2, 0.7, pos=2, expression(p[1]), col="red")

Optimaler Äquivalenztest

Für das obige Testproblem gibt es einen gleichmäßig besten Niveau-\(\alpha\)-Test. Wir lehnen \(H_0\) ab, falls \(c_1<S_n<c_2\) gilt. Wir akzeptieren \(H_0\), falls \(S_n<c_1\) oder \(S_n>c_2\) gilt. Im Fall \(S_n=c_1\) lehnen wir \(H_0\) mit Wahrscheinlichkeit \(\gamma_1\) ab, im Fall \(S_n=c_2\) lehnen wir \(H_0\) mit Wahrscheinlichkeit \(\gamma_2\) ab.

Es bezeichne \(p\mapsto G(p)\) die Gütefunktion des optimalen Tests. Dann gelten \(G(p_1)=\alpha\) und \(G(p_2)=\alpha\). Diese beiden Gleichungen legen die Konstanten \(0\le c_1<c_2\le n\) und \(0\le \gamma_1,\gamma_2<1\) eindeutig fest. Sie können durch “Ausprobieren” bestimmt werden.

Das R-Paket EQUIVNONINF implementiert die Algorithmen aus dem Buch von Wellek für Äquivalenztests. Das Paket kann über Packages/Install im rechten unteren Fenster installiert werden. Der Befehl für den Binomial-Äquivalenztest ist bi1st(). Der Algorithmus funktioniert nur für hinreichend kleine Stichproben.

library(EQUIVNONINF)
bi1st(0.05, 100, .6, .8)

 alpha = 0.05    n = 100    P1 = 0.6    P2 = 0.8    C1 = 68    C2 = 73 
 gam1 = 0.5965293    gam2 = 0.8318769    POWNONRD = 0.3367418    POW = 0.4429064

Es ergeben sich dieselben kritischen Werte \(c_1,c_2\) wie oben. Das Programm gibt auch die Güte des Tests bei \(p=(p_1+p_2)/2\) an. Für POWNONRD wird die nicht randomisierte Version des Tests verwendet, bei POW wird die randomisierte Version des Tests verwendet.

Anstatt das obige Programmpaket zu verwenden, kann man diese Funktion auch selbst programmieren, indem man den Code für den optimalen zweiseitigen Binomialtest anpasst.

Optimaler Äquivalenz-Test für normalverteilte Zufallsgrößen

Wir beschreiben den optimalen Äquivalenztest für den Erwartungswert einer normalverteilten Zufallsgröße bei bekannter Varianz. Der optimale Inäquivalenztest kann analog entwickelt werden. Dasselbe gilt für die klassischen einseitigen und zweiseitigen Tests auf Erwartungswert einer Normalverteilung bei bekannter Varianz (\(z\)-Tests).

Seien \(X_1, \ldots, X_n\) unabhängige wie \(N(\mu,\sigma^2)\) verteilte Zufallsvariablen. Wir nehmen an, dass \(\mu\) unbekannt ist und dass \(\sigma^2\) bekannt ist. Wir wollen auf das Hypothesenpaar

\[ H_0:|\mu-\mu_0|\ge \delta \qquad H_1: |\mu-\mu_0|<\delta \] testen. Hierzu verwenden wir die Test-Statistik

\[ T_n=(Y_n)^2 \ , \qquad Y_n=\frac{\overline{X}_n-\mu_0}{\sigma/\sqrt{n}} \ , \] mit \(\overline{X}_n=(X_1+\ldots+X_n)/n\) das arithmetische Mittel der Stichprobe. Wir lehnen \(H_0\) ab, falls \(T_n\) zu klein ist. Es handelt sich also um einen linksseitigen Test für den Parameter \(|\mu-\mu_0|\).

Die Zufallsgröße \(T_n\) ist das Quadrat der Zufallsgröße \(Y_n\). Die Zufallsgröße \(Y_n\) ist normalverteilt mit Erwartungswert \((\mu-\mu_0)/(\sigma/\sqrt{n})\) und Varianz \(1\). Die Verteilung von \(T_n\) heißt nichtzentrale \(\chi^2\)-Verteilung mit einem Freiheitsgrad und Nichtzentralitätsparameter

\[ \lambda=\left(\frac{\mu-\mu_0}{\sigma/\sqrt{n}}\right)^2 \] Wir lehnen also \(H_0\) ab, falls \(T_n< q\) gilt, mit \(q\) das \(\alpha\)-Quantil dieser Verteilung. Die Verteilung ist in R implementiert.

Wir wenden diesen Äquivalenz-Test auf die Frage nach der Wahrscheinlichkeit für den 29. Februar als Geburtstag an. Dies ist zulässig, weil beim Münzwurf die relative Häufigkeit der Erfolge bei großer Stichprobe in guter Näherung normalverteilt ist. Dies führt auf folgenden Code.

# Aequivalenztest fuer normalverteilte Zufallsgroessen
# Beispiel: p_0 fuer Geburtstage am 29. Februar?
p_0=(4*365+1)^(-1)
# Daten aus Nuernberg 2012
n=510400; k=334
# Normalapproximation zulaessig, falls n*p_0*(1-p_0)>10
n*p_0*(1-p_0)

[1] 349.1106

# waehle fuer delta relativen Fehler von 10 Prozent
del=p_0*0.1
# alpha-Quantil der nichtzentralen chi2-Verteilung
# 1 Freiheitsgrad (df), Nichtzentralitaetsparameter lam
alpha=0.05
sig=sqrt(p_0*(1-p_0)); lam=(del*sqrt(n)/sig)^2
q=qchisq(alpha, df=1, ncp=lam)
# Wert der Teststatistik
xbar=k/n; sig=sqrt(p_0*(1-p_0))
T=((xbar-p_0)*sqrt(n)/sig)^2
# Lehne H0 ab, falls T<q 
T; q

[1] 0.6749011
[1] 0.1178393

Im Beispiel können wir \(H_0\) nicht ablehnen. Dies könnte daran liegen, dass die Stichprobe nicht groß genug ist. Es wäre nun instruktiv, beispielsweise die Zahlen mehrerer Großstädte zusammenzufassen.

Jenseits des Binomialtests

Tests auf Erwartungswert in R

Der Einfachheit soll die hypothetische Fragestellung untersucht werden, ob der mittlere Wert einer Messgröße aus einem wiederholten Experiment größer als \(2.1\) ist.

Anstatt reale Daten zu verwenden, erzeugen wir mit dem Computer Daten aus einer Normalverteilung. Wir vergessen danach, wie die Daten erzeugt wurden. Mit diesem Verfahren kann man simulieren, wie sich der Test verhält. Damit kann man ein intuitives Verständnis für die Anwendung des Tests entwickeln.

mu0=2.2; sig=2; n=200
# n Zufallszahlen aus N(mu0, sig^2)
x=rnorm(n, mu0, sig)

Wir nehmen an, dass die Daten in guter Näherung einer Normalverteilung folgen. Dies kann zum Beispiel mit einem qq-Plot oder einem Test auf Normalverteilung untersucht werden.

Für den \(z\)-Test bei bekannter Varianz gehen wir wie folgt vor. Aus Bequemlichkeit verwenden wir das Paket BSDA, welches über Packages/Install im rechten unteren Fenster installiert werden kann.

# z-Test
# lade Paket mit z-Test
library(BSDA)
# rechtsseitiger z-Test bei bekannter Varianz auf mu>2.1 zum Niveau alpha
alpha=0.05
z.test(x,mu=2.1, sigma.x=sig, alt="greater", conf.level = 1-alpha)

    One-sample z-Test

data:  x
z = 0.10213, p-value = 0.4593
alternative hypothesis: true mean is greater than 2.1
95 percent confidence interval:
 1.881825       NA
sample estimates:
mean of x 
 2.114443 

Für den \(t\)-Test bei unbekannter Varianz gehen wir wie folgt vor.

# t-Test
# rechtsseitiger t-Test mu>2.1 bei unbekannter Varianz zum Niveau alpha
alpha=0.05
t.test(x,mu=2.1, alt="greater", conf.level = 1-alpha)

    One Sample t-test

data:  x
t = 0.1058, df = 199, p-value = 0.4579
alternative hypothesis: true mean is greater than 2.1
95 percent confidence interval:
 1.888846      Inf
sample estimates:
mean of x 
 2.114443 

Anpassungstests in R für diskrete Verteilungen

Als Standard-Beispiel untersuchen wir die Ziehung der Lottozahlen. Daten hierzu finden sich unter dem folgenden Link.

# Haeufigkeiten der Lottozahlen 1 bis 49 (Stand 4.3.2024)
x=c(607,587,607,593,584,649,594,544,585,587,614,557,520,558,560,581,591,582,602,558,547,615,565,586,606,618,589,550,589,559,617,621,618,570,563,591,581,606,575,579,599,603,611,569,533,561,583,592,636)
# Liste der Haeufigkeiten
x

 [1] 607 587 607 593 584 649 594 544 585 587 614 557 520 558 560 581 591 582 602 558
[21] 547 615 565 586 606 618 589 550 589 559 617 621 618 570 563 591 581 606 575 579
[41] 599 603 611 569 533 561 583 592 636

Wir verwenden das statistische Modell, dass sich die 49 Lottozahlen wie Zufallszahlen aus einer unbekannten gegebenen Verteilung mit 49 Werten verhalten. Die Nullhypothese des \(\chi^2\)-Anpassungstests lautet, dass die unbekannten Wahrscheinlichkeiten für die Lottozahlen den vermuteten Wahrscheinlichkeiten \(1/49\) gleichen. Als Teststatistik wird eine Größe verwendet, welche die Abweichung der beobachteten Häufigkeiten von den erwarteteten Häufigkeiten quantifiziert.

Die Durchführung des Tests in R ist sehr einfach. Wir wählen das Niveau \(\alpha=0.05\).

# Referenz-Wahrscheinlichkeiten
p0=rep(1/49,49)
# chi2-Anpassungstest
chisq.test(x,p=p0)

    Chi-squared test for given probabilities

data:  x
X-squared = 57.225, df = 48, p-value = 0.17

Weil der \(p\)-Wert größer als \(\alpha\) ist, müseen wir die Nullhypothese akzeptieren. Die Daten liefern keinen Anlass daran zu zweifeln, dass jede Zahl mit der gleichen Wahrscheinlichkeit auftritt.

Für die Testentscheidung wurde die asymptotisch exakte \(\chi^2\)-Verteilung der Teststatistik verwendet, welche auf einer Normalapproximation beruht. Als Faustregel gilt, dass die Normalapproximation zulässig ist, falls für jede vermutete Wahrscheinlichkeit \(p_i\) die Ungleichung \(n p_i\ge 5\) gilt, mit \(n\) die Stichprobengröße. In unserem Beispiel ist die Normalapproximation also unproblematisch.

Man kann die exakte Verteilung der Teststatistik ohne großen Aufwand mit dem Computer simulieren. Damit kann man zum einen überprüfen, wie gut die Normalapproximation ist. Zum anderen kann man den Anpassungstest mit der simulierten exakten Verteilung durchführen.

Man kann die \(\chi^2\)-Approximation für die asymptotische Verteilung der Teststatistik verwenden, um einen asymptotisch exakten Äquivalenztest auf Anpassung zu entwickeln.

Unabhängigkeitstests in R für diskrete Verteilungen

Wir beschränken uns auf den einfachsten Fall von Zufallsgrößen mit nur zwei Merkmalen und besprechen ein fiktives Beispiel.

Es soll zum Niveau \(\alpha=0.05\) untersucht werden, ob Smartphone-Nutzung einen signifikanten Einfluss auf schulische Leistungen hat. Hierzu wurde in der 12. Jahrgangsstufe eine Mathematik-Klausur gestellt, an der 87 Schüler:innen teilgenommen haben. Die Schüler:innen haben über einen Zeitraum von einem Monat vor der Klausur ihre Smartphone-Nutzungsdauer gemessen. Ihr Klausurergebnis war entweder besser (B) oder nicht besser (S) als das durchschnittliche Ergebnis. Ihre Smartphone-Nutzungsdauer war entweder niedrig (N, weniger als 2h/Tag) oder hoch (H, mindestens 2h/Tag).

Das Ergebnis ist in der folgenden Tabelle angegeben.

M=cbind(c(19,19),c(16,33))
rownames(M)=c("B", "S")
colnames(M)=c("N", "H")
M

   N  H
B 19 16
S 19 33

Wir verwenden das statistische Modell, dass die beiden Merkmale sich wie Zufallsgrößen mit fester unbekannter gemeinsamer Verteilung verhalten. Beim Unabhängigkeitstest lautet die Nullhypothese, dass die Merkmale unabhängig sind. Falls wir \(H_0\) ablehnen, haben wir eine signifikante Verletzung von Unabhängigkeit erkannt. Falls wir \(H_0\) akzeptieren, liefern die Daten keinen Hinweis darauf, dass die Annahme der Unabhängigkeit der Merkmale verletzt sein könnte.

Aufgrund der relativ großen Stichprobe sollten wir erwarten, dass eine Normalapproximation zulässig ist. Wir verwenden hierzu die Faustregel, dass jeder Eintrag in der Matrix größer als 5 ist. Also gehen wir davon aus, dass ein \(\chi^2\)-Anpassungstest zulässig ist.

# chi2-Unabhaengigkeitstest
chisq.test(M)

    Pearson's Chi-squared test with Yates' continuity correction

data:  M
X-squared = 2.0055, df = 1, p-value = 0.1567

Es ist vielleicht überraschend, dass die Unabhängigkeitshypothese akzeptiert werden muss. Die Stichprobe ist nicht groß genug, um eine Abweichung von Unabhängigkeit zu erkennen.

Der sogenannte exakte Test von Fisher hätte dasselbe Ergebnis geliefert. Seine Teststatistik ist der linke obere Eintrag in der Matrix. Die Verteilung der Teststatistik ist bei unabhängigen Merkmalen und festen Randhäufigkeiten hypergeometrisch. (Letztere Annahme ist für die Anwendung unproblematisch, weil man den Ergebnisraum in Teilräume mit festen Randhäufigkeiten zerlegen kann.) Es wird der \(p\)-Wert für das MPVT-Verfahren ausgegeben. Also ist nicht der optimale Test implementiert!

# exakter Test von Fisher
fisher.test(M)

    Fisher's Exact Test for Count Data

data:  M
p-value = 0.1252
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.7895282 5.4037752
sample estimates:
odds ratio 
  2.044939 

Tatsächlich implementiert der R-Befehl ein verwandtes Testproblem. Man kann sich dafür interessieren, ob zwei unabhängige Merkmale mit denselben Wahrscheinlichkeiten auftreten. Diese Frage nach Homogenität unabhängiger Merkmale unterscheidet sich von der Frage nach Unabhängigkeit der beiden Merkmale. Als Teststatistik wird das Verhältnis der geschätzten Chancen (odds) für überdurchschnittliches Abschneiden bei niedriger und hoher Nutzungsdauer verwendet. Der R-Befehl gibt auch ein Konfidenzintervall für das unbekannte Chancenverhältnis \(\varrho\) an. Ein Niveau-\(\alpha\)-Test auf \(H_0:\varrho=\varrho_0\) gegen \(H_1:\varrho\neq\varrho_0\) kann hieraus auch mit folgender Entscheidungsregel erhalten werden: Akzeptiere \(H_0\) zum Niveau \(\alpha\), falls \(\varrho_0\) im \((1-\alpha)\)-Konfidenzintervall für \(\varrho\) liegt. Tatsächlich gilt, dass bei festen Randhäufigkeiten die Bedingung der Unabhängigkeit gerade durch \(\varrho=1\) beschrieben wird. Für Details siehe das Buch von Lehmann und Romano, Kapitel 4.5 und 4.6.

# h_N Haeufigkeit fuer besseres Abschneiden bei niedriger Dauer
h_N=M["B","N"]/(M["B","N"]+M["S","N"])
# h_M Haeufigkeit fuer besseres Abschneiden bei hoher Dauer
h_H=M["B","H"]/(M["B","H"]+M["S","H"])
# Wert des Chancen-Verhaeltnisses (odds ratio)
(h_N/(1-h_N))/(h_H/(1-h_H))

[1] 2.0625