09_FairerVergleich__20220531: Kapitel 9: Fairer Vergleich in der Rueckmeldung
In LSAmitR: Daten, Beispiele und Funktionen zu 'Large-Scale Assessment mit R'
Kapitel 9 R Documentation
Kapitel 9: Fairer Vergleich in der Rueckmeldung
Description
Das ist die Nutzerseite zum Kapitel 9, Fairer Vergleich in der
Rückmeldung, im Herausgeberband Large-Scale Assessment mit R: Methodische
Grundlagen der österreichischen Bildungsstandardüberprüfung.
Im Abschnitt Details werden die im Kapitel verwendeten R-Syntaxen zur
Unterstützung für Leser/innen kommentiert und dokumentiert.
Im Abschnitt Examples werden die R-Syntaxen des Kapitels vollständig
wiedergegeben und gegebenenfalls erweitert.
Details
Vorbereitungen
Der zur Illustration verwendete Datensatz dat beinhaltet (imputierte)
aggregierte Leistungs- und Hintergrunddaten von 244 Schulen, bestehend aus 74
ländlichen allgemeinbildenden Pflichtschulen (APS, Stratum 1), 69 städtischen
APS (Stratum 2), 52 ländlichen allgemeinbildenden höheren Schulen (AHS, Stratum
3) und 49 städtischen AHS (Stratum 4).
Im Kapitel wird zur Bildung von Interaktionseffekten und quadratischen Termen
der Kovariaten eine neue Funktion covainteraction verwendet.
data(datenKapitel09)
dat <- datenKapitel09
covainteraction <- function(dat,covas,nchar){
for(ii in 1:(length(covas))){
vv1 <- covas[ii]
# Interaktion von vv1 mit sich selbst
subname1 <- substr(vv1,1,nchar)
intvar <- paste0(subname1, subname1)
if(vv1 == covas[1]){
dat.int <- dat[,vv1]*dat[,vv1];
newvars <- intvar } else {
dat.int <- cbind(dat.int,dat[,vv1]*dat[,vv1]);
newvars <- c(newvars,intvar)
}
# Interaktion von vv1 mit restlichen Variablen
if(ii < length(covas)){
for(jj in ((ii+1):length(covas))){
vv2 <- covas[jj]
subname2 <- substr(vv2,1,nchar)
intvar <- paste0(subname1, subname2)
newvars <- c(newvars, intvar)
dat.int <- cbind(dat.int,dat[,vv1]*dat[,vv2])
}
}
}
dat.int <- data.frame(dat.int)
names(dat.int) <- newvars
return(dat.int)
}
Abschnitt 9.2.5.1: Kovariaten und Interaktionsterme
Listing 1: Kovariatenauswahl und z-Standardisierung
Als Variablen zur Beschreibung von Kontext und Schülerzusammensetzung in den
Schulen werden in diesem Beispiel die logarithmierte Schulgröße groesse,
der Anteil an Mädchen female, der Anteil an Schülerinnen und Schülern mit
Migrationshintergrund mig und der mittlere sozioökonomische Status (SES)
sozstat eingeführt.
Die abhängige Variable ist die aggregierte
Schülerleistung der Schule TWLE. Alle Kovariaten (vars) werden
zunächst z-standardisiert (zvars).
vars <- c("groesse","female","mig","sozstat")
zvars <- paste0("z",vars)
dat[,zvars] <- scale(dat[,vars],scale = TRUE)
Listing 2: Interaktionen bilden, z-standardisieren
Zur Optimierung der Modellspezifikation werden Interaktionseffekte und
quadratische Terme der Kovariaten gebildet, dann z-standardisiert und in den
Gesamtdatensatz hinzugefügt.
Die neuen Variablennamen sind in der Liste intvars aufgelistet.
dat1 <- LSAmitR::covainteraction(dat = dat,covas = zvars,nchar = 4)
intvars <- names(dat1) # Interaktionsvariablen
dat1[,intvars] <- scale(dat1[,intvars],scale = TRUE)
dat <- cbind(dat,dat1)
Listing 3: Haupt- und Interaktionseffekte
maineff <- zvars # Haupteffekte
alleff <- c(zvars,intvars) # Haupt- und Interaktionseffekte
Abschnitt 9.2.5.2: OLS-Regression
Listing 4: OLS-Regression mit Haupteffekten
Das OLS-Regressionsmodell mit den Haupteffekten als Modellprädiktoren
(ols.mod1) für Schulen im Stratum (st) 4
fm.ols1 <- paste0("TWLE ~ ",paste(maineff,collapse=" + "))
fm.ols1 <- as.formula(fm.ols1) # Modellgleichung
st <- 4
pos <- which(dat$stratum == st) # Schulen im Stratum st
ols.mod1 <- lm(formula = fm.ols1,data = dat[pos,]) # Regression
Abschnitt 9.2.5.3: Lasso-Regression
Listing 5: Datenaufbereitung
Für die Durchführung der Lasso-Regression wird das R-Paket glmnet
(Friedman et al., 2010) eingesetzt. Die Kovariatenmatrix (Z) sowie der
Vektor der abhängigen Leistungswerte (Y) müssen definiert werden.
library(glmnet)
Z <- as.matrix(dat[pos,alleff]) # Kovariatenmatrix
Y <- dat$TWLE[pos] # Abhängige Variable
Listing 6: Bestimmung Teilmengen für Kreuzvalidierung,
Lasso-Regression
Das Lasso-Verfahren wird mit der Funktion cv.glmnet() durchgeführt.
Zur Auswahl eines optimalen shrinkage λ wird das Verfahren der
K-fachen Kreuzvalidierung verwendet.
Die Schulstichprobe wird durch ID-Zuweisung (foldid) verschiedenen
Teilmengen zugewiesen.
nid <- floor(length(pos)/3) # Teilmengen definieren
foldid <- rep(c(1:nid),3,length.out=length(pos)) # Zuweisung
lasso.mod2 <- glmnet::cv.glmnet(x=Z,y=Y,alpha = 1, foldid = foldid)
Listing 7: Erwartungswerte der Schulen
Entsprechend lambda.min werden die Regressionskoeffizienten und die
entsprechenden Erwartungswerte der Schulen bestimmt.
lasso.pred2 <- predict(lasso.mod2,newx = Z,s="lambda.min")
dat$expTWLE.Lasso2[pos] <- as.vector(lasso.pred2)
Listing 8: Bestimmung R^2
Bestimmung des aufgeklärten Varianzanteils der Schulleistung R^2.
varY <- var(dat$TWLE[pos])
varY.lasso.mod2 <- var(dat$expTWLE.Lasso2[pos])
R2.lasso.mod2 <- varY.lasso.mod2/varY
Abschnitt 9.2.5.4: Nichtparametrische Regression
Listing 9: Distanzberechnung
Der erste Schritt zur Durchführung einer nichtparametrischen Regression ist die
Erstellung der Distanzmatrix zwischen Schulen. In diesem Beispiel wird die
euklidische Distanz als Distanzmaß berechnet, alle standardisierten Haupteffekte
sind eingeschlossen. Außerdem setzen wir die Gewichte von allen Kovariaten
(gi) auf 1. dfr.i in diesem Beispiel ist die Distanzmatrix für
erste Schule im Stratum.
N <- length(pos) # Anzahl Schulen im Stratum
schools <- dat$idschool[pos] # Schulen-ID
i <- 1
# Teildatensatz von Schule i
dat.i <- dat[pos[i],c("idschool","TWLE",maineff)]
names(dat.i) <- paste0(names(dat.i),".i")
# Daten der Vergleichsschulen
dat.vgl <- dat[pos[-i],c("idschool","TWLE",maineff)]
index.vgl <- match(dat.vgl$idschool,schools)
# Daten zusammenfügen
dfr.i <- data.frame("index.i"=i,dat.i,"index.vgl"=index.vgl,
dat.vgl, row.names=NULL)
# Distanz zur Schule i
dfr.i$dist <- 0
gi <- c(1,1,1,1)
for(ii in 1:length(maineff)){
vv <- maineff[ii]
pair.vv <- grep(vv, names(dfr.i), value=T)
dist.vv <- gi[ii]*((dfr.i[,pair.vv[1]]-dfr.i[,pair.vv[2]])^2)
dfr.i$dist <- dfr.i$dist + dist.vv }
Listing 10: H initiieren
p(x) = λ^x exp(-λ)/x!.
Die Gewichte w_{ik} für jedes Paar (i, k) von Schulen werden mithilfe der
Distanz, der Gauß’schen Kernfunktion (dnorm) als Transformationsfunktion
sowie einer schulspezifischen Bandweite h_i berechnet. Die Auswahl des optimalen
Werts \hat{h_i} für jede Schule i erfolgt nach Vieu (1991). Zunächst wird
ein Vektor H so gewählt, dass der optimale Wert \hat{h_i} größer
als der kleinste und kleiner als der größte Wert in H ausfällt.
Je kleiner das Intervall zwischen den Werten in H ist, desto
wahrscheinlicher ist, dass ein Listenelement den optimalen Wert erlangt.
Auf der anderen Seite korrespondiert die Rechenzeit mit der Länge von H.
Gemäß der Größe der Vergleichsgruppe wählen wir eine Länge von 30 für H,
zusätzlich wird ein sehr großer Wert (100000) für die Fälle hinzugefügt, bei
denen alle Gewichte beinahe gleich sind.
d.dist <- max(dfr.i$dist)-min(dfr.i$dist)
H <- c(seq(d.dist/100,d.dist,length=30),100000)
V1 <- length(H)
# Anzahl Vergleichsschulen
n <- nrow(dfr.i)
Listing 11: Leave-One-Out-Schätzer der jeweiligen Vergleichsschule k
nach h in H
Auf Basis aller Werte in H und dem jeweils entsprechenden Gewicht w_{ik}
(wgt.h) werden die Leave-One-Out-Schätzer der jeweiligen Vergleichsschule
(pred.k) berechnet.
sumw <- 0*H # Vektor w_{ik} initiieren, h in H
av <- "TWLE"
dfr0.i <- dfr.i[,c("idschool",av)]
# Schleife über alle h-Werte
for (ll in 1:V1 ){
h <- H[ll]
# Gewicht w_{ik} bei h
dfr.i$wgt.h <- dnorm(sqrt(dfr.i$dist), mean=0, sd=sqrt(h))
# Summe von w_{ik} bei h
sumw[which(H==h)] <- sum(dfr.i$wgt.h)
# Leave-one-out-Schätzer von Y_k
for (k in 1:n){
# Regressionsformel
fm <- paste0(av,"~",paste0(maineff,collapse="+"))
fm <- as.formula(fm)
# Regressionsanalyse ohne Beitrag von Schule k
dfr.i0 <- dfr.i[-k,]
mod.k <- lm(formula=fm,data=dfr.i0,weights=dfr.i0$wgt.h)
# Erwartungswert anhand Kovariaten der Schule k berechnen
pred.k <- predict(mod.k, dfr.i)[k]
dfr0.i[k,paste0( "h_",h) ] <- pred.k
}}
# Erwartungswerte auf Basis verschiedener h-Werte
dfr1 <- data.frame("idschool.i"=dfr.i$idschool.i[1],"h"=H )
Listing 12: Kreuzvalidierungskriterium nach h in H
Zur Berechnung der Kreuzvalidierungskriterien (CVh, je kleiner, desto
reliabler sind die Schätzwerte) für alle Werte h in H verwenden wir in
diesem Beispiel die Plug-in-Bandweite nach Altman und Leger (1995) (hAL),
die mit der Funktion ALbw() des R-Pakets kerdiest aufrufbar ist.
library(kerdiest)
hAL <- kerdiest::ALbw("n",dfr.i$dist) # Plug-in Bandweite
dfr.i$cross.h <- hAL
dfr.i$crosswgt <- dnorm( sqrt(dfr.i$dist), mean=0, sd = sqrt(hAL) )
# Kreuzvalidierungskriterium CVh
vh <- grep("h_",colnames(dfr0.i),value=TRUE)
for (ll in 1:V1){
dfr1[ll,"CVh"] <- sum( (dfr0.i[,av] - dfr0.i[,vh[ll]])^2 *
dfr.i$crosswgt) / n}
Listing 13: Bestimmung des optimalen Wertes von h
Der optimale Wert von h in H (h.min) entspricht dem mit dem
kleinsten resultierenden CVh.
dfr1$min.h.index <- 0
ind <- which.min( dfr1$CVh )
dfr1$min.h.index[ind] <- 1
dfr1$h.min <- dfr1$h[ind]
Listing 14: Kleinste Quadratsumme der Schätzfehler
Kleinste Quadratsumme der Schätzfehler der nichtparametrischen Regression mit
h=h.min.
dfr1$CVhmin <- dfr1[ ind , "CVh" ]
Listing 15: Effizienzsteigerung
Die Effizienz (Steigerung der Schätzungsreliabilität) der nichtparametrischen
Regression gegenüber der linearen Regression (äquivalent zu CVh bei
h=100000).
dfr1$eff_gain <- 100 * ( dfr1[V1,"CVh"] / dfr1$CVhmin[1] - 1 )
Listing 16: Durchführung der nichtparametrischen Regression
h <- dfr1$h.min[1] # h.min
dfr.i$wgt.h <- dnorm(sqrt(dfr.i$dist),sd=sqrt(h))/
dnorm(0,sd= sqrt(h)) # w_{ik} bei h.min
dfr.i0 <- dfr.i
# Lokale Regression
mod.ii <- lm(formula=fm,data=dfr.i0,weights=dfr.i0$wgt.h)
# Kovariaten Schule i
predM <- data.frame(dfr.i[1,paste0(maineff,".i")])
names(predM) <- maineff
pred.ii <- predict(mod.ii, predM) # Schätzwert Schule i
dat$expTWLE.np[match(dfr1$idschool.i[1],dat$idschool)] <- pred.ii
Abschnitt 9.2.7, Berücksichtigung der Schätzfehler
Der Erwartungsbereich wird nach der im Buch beschriebenen Vorgehensweise
bestimmt.
Listing 17: Bestimmung des Erwartungsbereichs
Bestimmung der Breite des Erwartungsbereichs aller Schulen auf Basis der
Ergebnisse der OLS-Regression mit Haupteffekten.
vv <- "expTWLE.OLS1" # Variablenname
mm <- "OLS1" # Kurzname des Modells
dfr <- NULL # Ergebnistabelle
# Schleife über alle möglichen Breite von 10 bis 60
for(w in 10:60){
# Variablen für Ergebnisse pro w
var <- paste0(mm,".pos.eb",w) # Position der Schule
var.low <- paste0(mm,".eblow",w) # Untere Grenze des EBs
var.upp <- paste0(mm,".ebupp",w) # Obere Grenze des EBs
# Berechnen
dat[,var.low] <- dat[,vv]-w/2 # Untere Grenze des EBs
dat[,var.upp] <- dat[,vv]+w/2 # Obere Grenze des EBs
# Position: -1=unterhalb, 0=innerhalb, 1=oberhalb des EBs
dat[,var] <- -1*(dat$TWLE < dat[,var.low]) + 1*(dat$TWLE > dat[,var.upp])
# Verteilung der Schulpositionen
tmp <- data.frame(t(matrix(prop.table(table(dat[,var])))))
names(tmp) <- c("unterhalb","innerhalb","oberhalb")
tmp <- data.frame("ModellxBereich"=var,tmp); dfr <- rbind(dfr,tmp) }
# Abweichung zur Wunschverteilung 25-50-25
dfr1 <- dfr
dfr1[,c(2,4)] <- (dfr1[,c(2,4)] - .25)^2
dfr1[,3] <- (dfr1[,3] - .5)^2
dfr1$sumquare <- rowSums(dfr1[,-1])
# Auswahl markieren
dfr$Auswahl <- 1*(dfr1$sumquare == min(dfr1$sumquare) )
Author(s)
Giang Pham, Alexander Robitzsch, Ann Cathrice George, Roman Freunberger
References
Pham, G., Robitzsch, A., George, A. C. & Freunberger, R. (2016).
Fairer Vergleich in der Rückmeldung.
In S. Breit & C. Schreiner (Hrsg.), Large-Scale Assessment mit R:
Methodische Grundlagen der österreichischen Bildungsstandardüberprüfung
(pp. 295–332). Wien: facultas.
See Also
Zu datenKapitel09, den im Kapitel verwendeten Daten.
Zurück zu Kapitel 8, Fehlende Daten und Plausible Values.
Zu Kapitel 10, Reporting und Analysen.
Zur Übersicht.
Examples
## Not run:
library(miceadds)
library(glmnet)
library(kerdiest)
covainteraction <- function(dat,covas,nchar){
for(ii in 1:(length(covas))){
vv1 <- covas[ii]
# Interaktion von vv1 mit sich selbst
subname1 <- substr(vv1,1,nchar)
intvar <- paste0(subname1, subname1)
if(vv1 == covas[1]){
dat.int <- dat[,vv1]*dat[,vv1];
newvars <- intvar } else {
dat.int <- cbind(dat.int,dat[,vv1]*dat[,vv1]);
newvars <- c(newvars,intvar)
}
# Interaktion von vv1 mit restlichen Variablen
if(ii < length(covas)){
for(jj in ((ii+1):length(covas))){
vv2 <- covas[jj]
subname2 <- substr(vv2,1,nchar)
intvar <- paste0(subname1, subname2)
newvars <- c(newvars, intvar)
dat.int <- cbind(dat.int,dat[,vv1]*dat[,vv2])
}
}
}
dat.int <- data.frame(dat.int)
names(dat.int) <- newvars
return(dat.int)
}
data(datenKapitel09)
dat <- datenKapitel09
# Platzhalter für Leistungsschätzwerte aller Modelle
dat$expTWLE.OLS1 <- NA
dat$expTWLE.OLS2 <- NA
dat$expTWLE.Lasso1 <- NA
dat$expTWLE.Lasso2 <- NA
dat$expTWLE.np <- NA
## -------------------------------------------------------------
## Abschnitt 9.2.5, Umsetzung in R
## -------------------------------------------------------------
# -------------------------------------------------------------
# Abschnitt 9.2.5.1, Listing 1: Kovariatenauswahl und
# z-Standardisierung
#
vars <- c("groesse","female","mig","sozstat")
zvars <- paste0("z",vars)
dat[,zvars] <- scale(dat[,vars],scale = TRUE)
# -------------------------------------------------------------
# Abschnitt 9.2.5.1, Listing 2:
#
# Interaktionen bilden, z-standardisieren
dat1 <- LSAmitR::covainteraction(dat = dat,covas = zvars,nchar = 4)
intvars <- names(dat1) # Interaktionsvariablen
dat1[,intvars] <- scale(dat1[,intvars],scale = TRUE)
dat <- cbind(dat,dat1)
# -------------------------------------------------------------
# Abschnitt 9.2.5.1, Listing 3: Modellprädiktoren: Haupt- und
# Interaktionseffekte
#
maineff <- zvars # Haupteffekte
alleff <- c(zvars,intvars) # Haupt- und Interaktionseffekte
# -------------------------------------------------------------
# Abschnitt 9.2.5.2, Listing 4: OLS-Regression mit Haupteffekten
#
fm.ols1 <- paste0("TWLE ~ ",paste(maineff,collapse=" + "))
fm.ols1 <- as.formula(fm.ols1) # Modellgleichung
st <- 4
pos <- which(dat$stratum == st) # Schulen im Stratum st
ols.mod1 <- lm(formula = fm.ols1,data = dat[pos,]) # Regression
# -------------------------------------------------------------
# Abschnitt 9.2.5.3, Listing 5: Lasso-Regression
# Datenaufbereitung
#
library(glmnet)
Z <- as.matrix(dat[pos,alleff]) # Kovariatenmatrix
Y <- dat$TWLE[pos] # Abhängige Variable
# -------------------------------------------------------------
# Abschnitt 9.2.5.3, Listing 6: Lasso-Regression
# Bestimmung Teilmengen für Kreuzvalidierung, Lasso-Regression
#
nid <- floor(length(pos)/3) # Teilmengen definieren
foldid <- rep(c(1:nid),3,length.out=length(pos)) # Zuweisung
lasso.mod2 <- glmnet::cv.glmnet(x=Z,y=Y,alpha = 1, foldid = foldid)
# -------------------------------------------------------------
# Abschnitt 9.2.5.3, Listing 7: Lasso-Regression
# Erwartungswerte der Schulen
#
lasso.pred2 <- predict(lasso.mod2,newx = Z,s="lambda.min")
dat$expTWLE.Lasso2[pos] <- as.vector(lasso.pred2)
# -------------------------------------------------------------
# Abschnitt 9.2.5.3, Listing 8: Lasso-Regression
# Bestimmung R^2
#
varY <- var(dat$TWLE[pos])
varY.lasso.mod2 <- var(dat$expTWLE.Lasso2[pos])
R2.lasso.mod2 <- varY.lasso.mod2/varY
# -------------------------------------------------------------
# Abschnitt 9.2.5.4, Listing 9: Nichtparametrische Regression
# Distanzberechnung zur Schule i (Stratum st)
#
N <- length(pos) # Anzahl Schulen im Stratum
schools <- dat$idschool[pos] # Schulen-ID
i <- 1
# Teildatensatz von Schule i
dat.i <- dat[pos[i],c("idschool","TWLE",maineff)]
names(dat.i) <- paste0(names(dat.i),".i")
# Daten der Vergleichsschulen
dat.vgl <- dat[pos[-i],c("idschool","TWLE",maineff)]
index.vgl <- match(dat.vgl$idschool,schools)
# Daten zusammenfügen
dfr.i <- data.frame("index.i"=i,dat.i,"index.vgl"=index.vgl,
dat.vgl, row.names=NULL)
# Distanz zur Schule i
dfr.i$dist <- 0
gi <- c(1,1,1,1)
for(ii in 1:length(maineff)){
vv <- maineff[ii]
pair.vv <- grep(vv, names(dfr.i), value=T)
dist.vv <- gi[ii]*((dfr.i[,pair.vv[1]]-dfr.i[,pair.vv[2]])^2)
dfr.i$dist <- dfr.i$dist + dist.vv }
# -------------------------------------------------------------
# Abschnitt 9.2.5.4, Listing 10: Nichtparametrische Regression
#
# H initiieren
d.dist <- max(dfr.i$dist)-min(dfr.i$dist)
H <- c(seq(d.dist/100,d.dist,length=30),100000)
V1 <- length(H)
# Anzahl Vergleichsschulen
n <- nrow(dfr.i)
# -------------------------------------------------------------
# Abschnitt 9.2.5.4, Listing 11: Nichtparametrische Regression
# Berechnung der Leave-One-Out-Schätzer der jeweiligen
# Vergleichsschule k nach h in H
#
sumw <- 0*H # Vektor w_{ik} initiieren, h in H
av <- "TWLE"
dfr0.i <- dfr.i[,c("idschool",av)]
# Schleife über alle h-Werte
for (ll in 1:V1 ){
h <- H[ll]
# Gewicht w_{ik} bei h
dfr.i$wgt.h <- dnorm(sqrt(dfr.i$dist), mean=0, sd=sqrt(h))
# Summe von w_{ik} bei h
sumw[which(H==h)] <- sum(dfr.i$wgt.h)
# Leave-one-out-Schätzer von Y_k
for (k in 1:n){
# Regressionsformel
fm <- paste0(av,"~",paste0(maineff,collapse="+"))
fm <- as.formula(fm)
# Regressionsanalyse ohne Beitrag von Schule k
dfr.i0 <- dfr.i[-k,]
mod.k <- lm(formula=fm,data=dfr.i0,weights=dfr.i0$wgt.h)
# Erwartungswert anhand Kovariaten der Schule k berechnen
pred.k <- predict(mod.k, dfr.i)[k]
dfr0.i[k,paste0( "h_",h) ] <- pred.k
}}
# Erwartungswerte auf Basis verschiedener h-Werte
dfr1 <- data.frame("idschool.i"=dfr.i$idschool.i[1],"h"=H )
# -------------------------------------------------------------
# Abschnitt 9.2.5.4, Listing 12: Nichtparametrische Regression
# Berechnung des Kreuzvalidierungskriteriums nach h in H
#
library(kerdiest)
hAL <- kerdiest::ALbw("n",dfr.i$dist) # Plug-in Bandweite
dfr.i$cross.h <- hAL
dfr.i$crosswgt <- dnorm( sqrt(dfr.i$dist), mean=0, sd = sqrt(hAL) )
# Kreuzvalidierungskriterium CVh
vh <- grep("h_",colnames(dfr0.i),value=TRUE)
for (ll in 1:V1){
dfr1[ll,"CVh"] <- sum( (dfr0.i[,av] - dfr0.i[,vh[ll]])^2 *
dfr.i$crosswgt) / n}
# -------------------------------------------------------------
# Abschnitt 9.2.5.4, Listing 13: Nichtparametrische Regression
# Bestimmung optimales Wertes von h (h.min)
#
dfr1$min.h.index <- 0
ind <- which.min( dfr1$CVh )
dfr1$min.h.index[ind ] <- 1
dfr1$h.min <- dfr1$h[ind]
# -------------------------------------------------------------
# Abschnitt 9.2.5.4, Listing 14: Nichtparametrische Regression
# Kleinste Quadratsumme der Schätzfehler
#
dfr1$CVhmin <- dfr1[ ind , "CVh" ]
# -------------------------------------------------------------
# Abschnitt 9.2.5.4, Listing 15: Nichtparametrische Regression
# Effizienzsteigerung berechnen
#
dfr1$eff_gain <- 100 * ( dfr1[V1,"CVh"] / dfr1$CVhmin[1] - 1 )
# -------------------------------------------------------------
# Abschnitt 9.2.5.4, Listing 16: Nichtparametrische Regression
# Durchführung der nichtparametrischen Regression bei h=h.min
#
h <- dfr1$h.min[1] # h.min
dfr.i$wgt.h <- dnorm(sqrt(dfr.i$dist),sd=sqrt(h))/
dnorm(0,sd= sqrt(h)) # w_{ik} bei h.min
dfr.i0 <- dfr.i
# Lokale Regression
mod.ii <- lm(formula=fm,data=dfr.i0,weights=dfr.i0$wgt.h)
# Kovariaten Schule i
predM <- data.frame(dfr.i[1,paste0(maineff,".i")])
names(predM) <- maineff
pred.ii <- predict(mod.ii, predM) # Schätzwert Schule i
dat[match(dfr1$idschool.i[1],dat$idschool), "expTWLE.np"] <- pred.ii
## -------------------------------------------------------------
## Abschnitt 9.2.5, Umsetzung in R, Ergänzung zum Buch
## -------------------------------------------------------------
# Korrelationen zwischen Haupteffekten
cor(dat[,maineff]) # gesamt
# Pro Stratum
for(s in 1:4) print(cor(dat[which(dat$stratum == s),maineff]))
# -------------------------------------------------------------
# Abschnitt 9.2.5.2, Ergänzung zum Buch
# OLS-Regression
#
# Modellgleichung nur mit Haupteffekten
fm.ols1 <- paste0("TWLE ~ ",paste(maineff,collapse=" + "))
fm.ols1 <- as.formula(fm.ols1)
# Modellgleichung mit Haupteffekten ohne zgroesse
fm.ols1a <- paste0("TWLE ~ ",paste(setdiff(maineff,c("zgroesse")),
collapse=" + "))
fm.ols1a <- as.formula(fm.ols1a)
# Modellgleichung mit Haupt- und Interaktionseffekten
fm.ols2 <- paste0("TWLE ~ ",paste(alleff,collapse=" + "))
fm.ols2 <- as.formula(fm.ols2)
# Ergebnistabelle über 4 Strata hinweg vorbereiten
tab1 <- data.frame("Variable"=c("(Intercept)",maineff))
tab2 <- data.frame("Variable"=c("(Intercept)",alleff))
# Durchführung: Schleife über vier Strata
for(st in 1:4){
# st <- 4
# Position Schulen des Stratums st im Datensatz
pos <- which(dat$stratum == st)
#---------------------------------
# OLS-Modell 1
# Durchführung
ols.mod1 <- lm(formula = fm.ols1,data = dat[pos,])
ols.mod1a <- lm(formula = fm.ols1a,data = dat[pos,])
# Modellergebnisse anzeigen
summary(ols.mod1)
summary(ols.mod1a)
# Erwartungswerte der Schulen
dat$expTWLE.OLS1[pos] <- fitted(ols.mod1)
# Ergebnisse in Tabelle speichern
par <- summary(ols.mod1)
tab.s <- data.frame(par$coef,R2=par$r.squared,R2.adj=par$adj.r.squared)
names(tab.s) <- paste0("stratum",st,
c("_coef","_SE","_t","_p","_R2","_R2.adj"))
tab1 <- cbind(tab1, tab.s)
# Durchführung OLS-Modell 2
ols.mod2 <- lm(formula = fm.ols2,data = dat[pos,])
# Modellergebnisse anzeigen
summary(ols.mod2)
# Erwartungswerte der Schulen
dat$expTWLE.OLS2[pos] <- fitted(ols.mod2)
# Ergebnisse in Tabelle speichern
par <- summary(ols.mod2)
tab.s <- data.frame(par$coef,R2=par$r.squared,R2.adj=par$adj.r.squared)
names(tab.s) <- paste0("stratum",st,
c("_coef","_SE","_t","_p","_R2","_R2.adj"))
tab2 <- cbind(tab2, tab.s)
}
# Daten Schule 1196 ansehen
dat[which(dat$idschool == 1196),]
# Schätzwerte nach ols.mod1 und ols.mod2 vergleichen
summary(abs(dat$expTWLE.OLS1 - dat$expTWLE.OLS2))
cor.test(dat$expTWLE.OLS1,dat$expTWLE.OLS2)
# Grafische Darstellung des Vergleich (Schule 1196 rot markiert)
plot(dat$expTWLE.OLS1,dat$expTWLE.OLS2,xlim=c(380,650),ylim=c(380,650),
col=1*(dat$idschool == 1196)+1,pch=15*(dat$idschool == 1196)+1)
abline(a=0,b=1)
# -------------------------------------------------------------
# Abschnitt 9.2.5.3, Ergänzung zum Buch
# Lasso-Regression
#
library(glmnet)
# Variablen für Erwartungswerte
dat$expTWLE.Lasso2 <- dat$expTWLE.Lasso1 <- NA
# Tabelle für Modellergebnisse
tab3 <- data.frame("Variable"=c("(Intercept)",maineff))
tab4 <- data.frame("Variable"=c("(Intercept)",alleff))
for(st in 1:4){
# st <- 4
# Position Schulen des Stratums st im Datensatz
pos <- which(dat$stratum == st)
#------------------------------------------------------------#
# Lasso-Regression mit den Haupteffekten
# Kovariatenmatrix
Z <- as.matrix(dat[pos,maineff])
# Abhängige Variable
Y <- dat$TWLE[pos]
# Kreuzvalidierung: Teilmengen definieren
nid <- floor(length(pos)/3)
# Schulen zu Teilmengen zuordnen
foldid <- rep(c(1:nid),3,length.out=length(pos))
# Regression
lasso.mod1 <- cv.glmnet(x=Z,y=Y,alpha = 1, foldid = foldid)
# Ergebnisse ansehen
print(lasso.mod1)
# Lasso-Koeffizienten bei lambda.min
print(lasso.beta <- coef(lasso.mod1,s="lambda.min"))
# Erwartungswerte der Schulen
lasso.pred1 <- predict(lasso.mod1,newx = Z,s="lambda.min")
dat$expTWLE.Lasso1[pos] <- as.vector(lasso.pred1)
# R2 bestimmen
varY <- var(dat$TWLE[pos])
varY.lasso.mod1 <- var(dat$expTWLE.Lasso1[pos])
print(R2.lasso.mod1 <- varY.lasso.mod1/varY)
# Ergebnistabelle
vv <- paste0("coef.stratum",st); tab3[,vv] <- NA
tab3[lasso.beta@i+1,vv] <- lasso.beta@x
vv <- paste0("lambda.stratum",st); tab3[,vv] <- lasso.mod1$lambda.min
vv <- paste0("R2.stratum",st); tab3[,vv] <- R2.lasso.mod1
#------------------------------------------------------------#
# Lasso-Regression mit Haupt- und Interaktionseffekten
# Kovariatenmatrix
Z <- as.matrix(dat[pos,alleff])
# Regression
lasso.mod2 <- cv.glmnet(x=Z,y=Y,alpha = 1, foldid = foldid)
# Ergebnisausdruck
print(lasso.mod2)
# Lasso-Koeffizienten bei lambda.min
print(lasso.beta <- coef(lasso.mod2,s="lambda.min"))
# Erwartungswerte der Schulen
lasso.pred2 <- predict(lasso.mod2,newx = Z,s="lambda.min")
dat$expTWLE.Lasso2[pos] <- as.vector(lasso.pred2)
# R2 bestimmen
varY.lasso.mod2 <- var(dat$expTWLE.Lasso2[pos])
R2.lasso.mod2 <- varY.lasso.mod2/varY
R2.lasso.mod2
# Ergebnistabelle
vv <- paste0("coef.stratum",st); tab4[,vv] <- NA
tab4[lasso.beta@i+1,vv] <- lasso.beta@x
vv <- paste0("lambda.stratum",st); tab4[,vv] <- lasso.mod2$lambda.min
vv <- paste0("R2.stratum",st); tab4[,vv] <- R2.lasso.mod2
}
# Regressionresiduen = Schätzung von SChul- und Unterrichtseffekt
dat$resTWLE.Lasso1 <- dat$TWLE - dat$expTWLE.Lasso1
dat$resTWLE.Lasso2 <- dat$TWLE - dat$expTWLE.Lasso2
# -------------------------------------------------------------
# Abschnitt 9.2.5.4, Ergänzung zum Buch
# Nichtparametrische Regression
#
#
# Achtung: Der nachfolgende Algorithmus benötigt viel Zeit!
#
av <- "TWLE" # Abhängige Variable
dfr3 <- NULL # Ergebnistabelle
# Variable für Leistungsschätzwerte
# Schleife über 4 Strata
for(st in 1:4){
# st <- 1
pos <- which(dat$stratum == st)
N <- length(pos)
schools <- dat$idschool[pos]
###
# Distanzmatrix dfr für alle Schulen im Stratum erstellen
dfr <- NULL
for (i in 1:N){
# i <- 1
# Teildatensatz von Schule i
dat.i <- dat[pos[i],c("idschool","TWLE",maineff)]
# Daten der Vergleichsgruppe
dat.vgl <- dat[pos[-i],c("idschool","TWLE",maineff)]
# Variablennamen von dat.vgl umbenennen
# names(dat.vgl) <- paste0("vgl.",names(dat.vgl))
# Variablennamen von dat.i umbenennen
names(dat.i) <- paste0(names(dat.i),".i")
# Daten zusammenfügen
index.vgl <- match(dat.vgl$idschool,schools)
dfr.i <- data.frame("index.i"=i,dat.i,
"index.vgl"=index.vgl,dat.vgl,
row.names=NULL)
# Distanz zur i
dfr.i$dist <- 0
gi <- c(1,1,1,1)
for(ii in 1:length(maineff)){
vv <- maineff[ii]
pair.vv <- grep(vv, names(dfr.i), value=T)
dist.vv <- gi[ii]*((dfr.i[,pair.vv[1]]-dfr.i[,pair.vv[2]])^2)
dfr.i$dist <- dfr.i$dist + dist.vv
}
print(i) ; flush.console()
dfr <- rbind( dfr , dfr.i )
}
dfr1 <- index.dataframe( dfr , systime=TRUE )
###
# h-Auswahl und Nichtparametrische Regression pro Schule i
dfr1.list <- list()
for (i in 1:N){
# i <- 1
dfr.i <- dfr[ dfr$index.i == i , ]
n <- nrow(dfr.i)
# Startwertliste für h initiieren
d.dist <- max(dfr.i$dist)-min(dfr.i$dist)
H <- c(seq(d.dist/100,d.dist,length=30),100000)
V1 <- length(H) # Anzahl der Startwerte in H
# Startwerte: Summe von w_ik
sumw <- 0*H
dfr0.i <- dfr.i[,c("idschool",av)]
# Schleife über alle h-Werte
for (ll in 1:V1 ){
h <- H[ll]
# Gewicht w_ik bei h
dfr.i$wgt.h <- dnorm(sqrt(dfr.i$dist), mean=0, sd=sqrt(h))
# Summe von w_ik bei h
sumw[which(H==h)] <- sum(dfr.i$wgt.h)
# Leave-one-out-Schätzer von Y_k
for (k in 1:n){
# Regressionsformel
fm <- paste0(av,"~",paste0(maineff,collapse="+"))
fm <- as.formula(fm)
# Regressionsanalyse ohne Beitrag von Schule k
dfr.i0 <- dfr.i[-k,]
mod.k <- lm(formula=fm,data=dfr.i0,weights=dfr.i0$wgt.h)
# Erwartungswert anhand Kovariaten der Schule k berechnen
pred.k <- predict(mod.k, dfr.i)[k]
dfr0.i[k,paste0( "h_",h) ] <- pred.k
}
print(paste0("i=",i,", h_",ll))
}
# Erwartungswerte auf Basis verschiedener h-Werte
dfr1 <- data.frame("idschool.i"=dfr.i$idschool.i[1],"h"=H )
# Berechnung des Kreuzvalidierungskriteriums
library(kerdiest)
hAL <- kerdiest::ALbw("n",dfr.i$dist) # Plug-in Bandbreite nach Altman und
# Leger
name <- paste0( "bandwidth_choice_school" , dfr.i$idschool.i[1] ,
"_cross.h_" , round2(hAL,1) )
# Regressionsgewichte auf Basis cross.h
dfr.i$cross.h <- hAL
dfr.i$crosswgt <- dnorm( sqrt(dfr.i$dist), mean=0, sd = sqrt(hAL) )
dfr.i <- index.dataframe( dfr.i , systime=TRUE )
# Kreuzvalidierungskriterium CVh
vh <- grep("h_",colnames(dfr0.i),value=TRUE)
for (ll in 1:V1){
# ll <- 5
dfr1[ll,"CVh"] <- sum( (dfr0.i[,av] - dfr0.i[,vh[ll]])^2 *
dfr.i$crosswgt) / n
print(ll)
}
# Bestimmung h.min
dfr1$min.h.index <- 0
ind <- which.min( dfr1$CVh )
dfr1$min.h.index[ind ] <- 1
dfr1$h.min <- dfr1$h[ind]
# Kleinste Quadratsumme der Schätzfehler
dfr1$CVhmin <- dfr1[ ind , "CVh" ]
# Effizienzsteigerung berechnen
dfr1$eff_gain <- 100 * ( dfr1[V1,"CVh"] / dfr1$CVhmin[1] - 1 )
# h auswählen
h <- dfr1$h.min[1]
# Gewichte anhand h berechnen
dfr.i$wgt.h <- dnorm( sqrt( dfr.i$dist ) , sd = sqrt( h) ) /
dnorm( 0 , sd = sqrt( h) )
dfr.i0 <- dfr.i
mod.ii <- lm(formula = fm,data = dfr.i0,weights = dfr.i0$wgt.h)
# Leistungsschätzwerte berechnen
predM <- data.frame(dfr.i[1,paste0(maineff,".i")])
names(predM) <- maineff
pred.ii <- predict( mod.ii , predM )
dfr1$fitted_np <- pred.ii
dfr1$h.min_sumwgt <- sum( dfr.i0$wgt.h )
dfr1$h_sumwgt <- sumw
# Leistungsschätzwerte zum Datensatz hinzufügen
dat$expTWLE.np[match(dfr1$idschool.i[1],dat$idschool)] <- pred.ii
dfr1.list[[i]] <- dfr1
}
###
# Ergebnisse im Stratum st zusammenfassen
dfr2 <- NULL
for(i in 1:length(dfr1.list)){
dat.ff <- dfr1.list[[i]]
dfr2.ff <- dat.ff[1,c("idschool.i","h.min","fitted_np","h.min_sumwgt",
"CVhmin","eff_gain")]
dfr2.ff$CVlinreg <- dat.ff[V1,"CVh"]
names(dfr2.ff) <- c("idschool","h.min","fitted_np","h.min_sumwgt",
"CVhmin","eff_gain","CVlinreg")
dfr2 <- rbind(dfr2, dfr2.ff)
print(i)
}
#---------------------------------------------------##
# R2 berechnen
varY <- var(dat$TWLE[pos])
varY.np <- var(dat$expTWLE.np[pos])
dfr2$R2.np <- varY.np/varY
#---------------------------------------------------##
# Zur Gesamtergebnistabelle
dfr3 <- rbind(dfr3,cbind("Stratum"=st,dfr2))
}
# Effizienz der NP-Regression gegenüber OLS-Regression
summary(dfr3$eff_gain)
table(dfr3$eff_gain > 5)
table(dfr3$eff_gain > 10)
table(dfr3$eff_gain > 20)
# Regressionsresiduen
dat$resTWLE.np <- dat$TWLE - dat$expTWLE.np
## -------------------------------------------------------------
## Abschnitt 9.2.6, Ergänzung zum Buch
## Ergebnisse im Vergleich
## -------------------------------------------------------------
# Output-Variablen
out <- grep("expTWLE",names(dat),value=T)
lt <- length(out)
# Korrelationsmatrix
tab <- tab1 <- as.matrix(round2(cor(dat[,out]),3))
# Varianzmatrix
tab2 <- as.matrix(round2(sqrt(var(dat[,out])),1))
tab3 <- matrix(NA,lt,lt)
# Differenzmatrix
for(ii in 1:(lt-1))
for(jj in (ii+1):lt) tab3[ii,jj] <- round2(mean(abs(dat[,out[jj]] -
dat[,out[ii]])),1)
tab4 <- matrix(NA,lt,lt)
# Differenzmatrix
for(ii in 1:(lt-1))
for(jj in (ii+1):lt) tab4[ii,jj] <- round2(sd(abs(dat[,out[jj]] -
dat[,out[ii]])),1)
# Ergebnistabelle
diag(tab) <- diag(tab2)
tab[upper.tri(tab)] <- tab3[upper.tri(tab3)]
# R2 Gesamt
varY <- var(dat$TWLE)
varexp.OLS1 <- var(dat$expTWLE.OLS1); R2.OLS1 <- varexp.OLS1/varY
varexp.OLS2 <- var(dat$expTWLE.OLS2); R2.OLS2 <- varexp.OLS2/varY
varexp.Lasso1 <- var(dat$expTWLE.Lasso1); R2.Lasso1 <- varexp.Lasso1/varY
varexp.Lasso2 <- var(dat$expTWLE.Lasso2); R2.Lasso2 <- varexp.Lasso2/varY
varexp.np <- var(dat$expTWLE.np); R2.np <- varexp.np/varY
R2 <- c(R2.OLS1,R2.OLS2,R2.Lasso1,R2.Lasso2,R2.np)
tab <- cbind(tab,R2)
# R2 pro Stratum
dat0 <- dat
for(st in 1:4){
# st <- 1
dat <- dat0[which(dat0$stratum == st),]
varY <- var(dat$TWLE)
varexp.OLS1 <- var(dat$expTWLE.OLS1); R2.OLS1 <- varexp.OLS1/varY
varexp.OLS2 <- var(dat$expTWLE.OLS2); R2.OLS2 <- varexp.OLS2/varY
varexp.Lasso1 <- var(dat$expTWLE.Lasso1); R2.Lasso1 <- varexp.Lasso1/varY
varexp.Lasso2 <- var(dat$expTWLE.Lasso2); R2.Lasso2 <- varexp.Lasso2/varY
varexp.np <- var(dat$expTWLE.np); R2.np <- varexp.np/varY
R2 <- c(R2.OLS1,R2.OLS2,R2.Lasso1,R2.Lasso2,R2.np)
tab <- cbind(tab,R2)
}
colnames(tab)[7:10] <- paste0("R2_stratum",1:4)
## -------------------------------------------------------------
## Abschnitt 9.2.7, Berücksichtigung der Schätzfehler
## -------------------------------------------------------------
# -------------------------------------------------------------
# Abschnitt 9.2.7, Listing 17: Bestimmung des Erwartungsbereichs
#
vv <- "expTWLE.OLS1" # Variablenname
mm <- "OLS1" # Kurzname des Modells
dfr <- NULL # Ergebnistabelle
# Schleife über alle möglichen Breite von 10 bis 60
for(w in 10:60){
# Variablen für Ergebnisse pro w
var <- paste0(mm,".pos.eb",w) # Position der Schule
var.low <- paste0(mm,".eblow",w) # Untere Grenze des EBs
var.upp <- paste0(mm,".ebupp",w) # Obere Grenze des EBs
# Berechnen
dat[,var.low] <- dat[,vv]-w/2 # Untere Grenze des EBs
dat[,var.upp] <- dat[,vv]+w/2 # Obere Grenze des EBs
# Position: -1=unterhalb, 0=innerhalb, 1=oberhalb des EBs
dat[,var] <- -1*(dat$TWLE < dat[,var.low]) + 1*(dat$TWLE > dat[,var.upp])
# Verteilung der Schulpositionen
tmp <- data.frame(t(matrix(prop.table(table(dat[,var])))))
names(tmp) <- c("unterhalb","innerhalb","oberhalb")
tmp <- data.frame("ModellxBereich"=var,tmp); dfr <- rbind(dfr,tmp) }
# Abweichung zur Wunschverteilung 25-50-25
dfr1 <- dfr
dfr1[,c(2,4)] <- (dfr1[,c(2,4)] - .25)^2
dfr1[,3] <- (dfr1[,3] - .5)^2
dfr1$sumquare <- rowSums(dfr1[,-1])
# Auswahl markieren
dfr$Auswahl <- 1*(dfr1$sumquare == min(dfr1$sumquare) )
# -------------------------------------------------------------
# Abschnitt 9.2.7, Ergänzung zum Buch
# Bestimmung des Erwartungsbereichs
#
# Ergebnisse aller Schulen werden aus Ursprungsdatensatz geladen.
dat <- datenKapitel09
# Liste der Erwartungswerte-Variablen
exp.vars <- grep("expTWLE",names(dat),value=T)
# Modellnamen
m.vars <- gsub("expTWLE.","",exp.vars, fixed = TRUE)
# Liste der Ergebnistabelle
list0 <- list()
# Ergebnisse
tab.erg <- NULL
# Schleife über alle Erwartungswerte aller Modelle
for(ii in 1:length(exp.vars)){
# ii <- 1
vv <- exp.vars[ii]
mm <- m.vars[ii]
# Ergebnistabelle
dfr <- NULL
# Schleife über alle möglichen Breite von 10 bis 60
for(w in 10:60){
# eb <- 10
var <- paste0(mm,".pos.eb",w) # Position der Schule
var.low <- paste0(mm,".eblow",w) # Untere Grenze des EBs
var.upp <- paste0(mm,".ebupp",w) # Obere Grenze des EBs
# Untere Grenze des EBs = Erwartungswert - w/2
dat[,var.low] <- dat[,vv]-w/2
# Obere Grenze des EBs = Erwartungswert + w/2
dat[,var.upp] <- dat[,vv]+w/2
# Position der Schule bestimmen
# -1 = unterhalb, 0 = innterhalb, 1 = oberhalb des EBs
dat[,var] <- -1*(dat$TWLE < dat[,var.low]) + 1*(dat$TWLE > dat[,var.upp])
# Verteilung der Positionen
tmp <- data.frame(t(matrix(prop.table(table(dat[,var])))))
names(tmp) <- c("unterhalb","innerhalb","oberhalb")
tmp <- data.frame("ModellxBereich"=var,tmp)
dfr <- rbind(dfr,tmp)
}
# Vergleich mit Wunschverteilung 25-50-25
dfr1 <- dfr
dfr1[,c(2,4)] <- (dfr1[,c(2,4)] - .25)^2
dfr1[,3] <- (dfr1[,3] - .5)^2
dfr1$sumquare <- rowSums(dfr1[,-1])
# Auswahl markieren
dfr$Auswahl <- 1*(dfr1$sumquare == min(dfr1$sumquare) )
# Zum Liste hinzufügen
list0[[ii]] <- dfr
print(dfr[which(dfr$Auswahl == 1),])
tab.erg <- rbind(tab.erg, dfr[which(dfr$Auswahl == 1),])
}
# Nur gewählte Ergebnisse im Datensatz beibehalten
all.vars <- grep("eb",names(dat),value=T)
# Untere und Obere Grenze mit speichern
eb.vars <- tab.erg[,1]
low.vars <- gsub("pos.eb","eblow",eb.vars)
upp.vars <- gsub("pos.eb","ebupp",eb.vars)
del.vars <- setdiff(all.vars, c(eb.vars,low.vars,upp.vars))
dat <- dat[,-match(del.vars,names(dat))]
## -------------------------------------------------------------
## Appendix: Abbildungen
## -------------------------------------------------------------
# -------------------------------------------------------------
# Abbildung 9.4
#
# Koeffizienten bei der ersten 50 lambdas ausdrucken
# Stratum 4
lambda <- lasso.mod2$lambda[1:50]
a <- round2(lambda,2)
a1 <- a[order(a)]
L <- length(a)
dfr <- NULL
for(ll in 1:L){
dfr.ll <- as.matrix(coef(lasso.mod2,newx = Z,s=a[ll] ))
colnames(dfr.ll) <- paste0("a_",ll)
dfr.ll <- data.frame("coef"=rownames(dfr.ll),dfr.ll)
rownames(dfr.ll) <- NULL
if(ll == 1) dfr <- dfr.ll else dfr <- merge(dfr, dfr.ll)
}
# Ohne Intercept
dfr <- dfr[-1,]
rownames(dfr) <- 1:nrow(dfr)
cl <- colors()
cl <- grep("grey",cl,value=T)
# Umgekehrte Reihenfolge
dfr1 <- dfr
for(x in 2:(L+1)) {dfr1[,x] <- dfr[,(L+3)-x];
names(dfr1)[x] <- names(dfr)[(L+3)-x]}
###
plot(x = log(a), y = rep(0,L), xlim = rev(range(log(a))), ylim=c(-20,22),
type = "l", xaxt ="n", xlab = expression(paste(lambda)),
ylab="Geschätzte Regressionskoeffizienten")
axis(1, at=log(a), labels=a,cex=1)
tmp <- nrow(dfr)
for(ll in 1:tmp){
# ll <- 1
lines(x=log(a),y=dfr[ll,2:(L+1)],type="l",pch=15-ll,col=cl[15-ll])
points(x=log(a),y=dfr[ll,2:(L+1)],type="p",pch=15-ll)
legend(x=2.8-0.7*(ll>tmp/2),y=25-2*(ifelse(ll>7,ll-7,ll)),
legend =dfr$coef[ll],pch=15-ll,bty="n",cex=0.9)
}
# Kennzeichung der gewählten lambda
v <- log(lasso.mod2$lambda.min)
lab2 <- expression(paste("ausgewähltes ",lambda," = .43"))
text(x=v+0.6,y=-8,labels=lab2)
abline(v = v,lty=2,cex=1.2)
# -------------------------------------------------------------
# Abbildung 9.5
# Auswahl Lambda anhand min(cvm)
#
xlab <- expression(paste(lambda))
plot(lasso.mod2, xlim = rev(range(log(lambda))),
ylim=c(550,1300),xlab=xlab,xaxt ="n",
ylab = "Mittleres Fehlerquadrat der Kreuzvalidierung (cvm)",
font.main=1,cex.main=1)
axis(1, at=log(a), labels=a,cex=1)
lab1 <- expression(paste(lambda," bei min(cvm)"))
text(x=log(lasso.mod2$lambda.min)+0.5,y=max(lasso.mod2$cvm)-50,
labels=lab1,cex=1)
lab2 <- expression(paste("(ausgewähltes ",lambda," = .43)"))
text(x=log(lasso.mod2$lambda.min)+0.6,y=max(lasso.mod2$cvm)-100,
labels=lab2,cex=1)
abline(v = log(lasso.mod2$lambda.min),lty=2)
text(x=log(lasso.mod2$lambda.min)-0.3,y = min(lasso.mod2$cvm)-30,
labels="min(cvm)",cex=1 )
abline(h = min(lasso.mod2$cvm),lty=2)
text <- expression(paste("Anzahl der Nicht-null-Koeffizienten (",
lambda," entsprechend)"))
mtext(text=text,side=3,line=3)
# -------------------------------------------------------------
# Abbildung 9.6
# Rohwert-Schätzwert Schule 1196 & 1217 im Vergleich
#
id <- c(1196, 1217)
par(mai=c(1.2,3,1,.5))
plot(x=rep(NA,2),y=c(1:2),xlim=c(470,610),yaxt ="n",type="l",
xlab="Erwartungswerte je nach Modell und Schulleistung",ylab="")
legend <- c("Schulleistung (TWLE)",paste0("", c("OLS1","OLS2","Lasso1",
"Lasso2","NP"),
"-Modell"))
axis(2, at=c(seq(1,1.4,0.08),seq(1.6,2,0.08)), las=1,cex=0.7,
labels=rep(legend,2))
text <- paste0("Schule ",id)
mtext(text=text,side=2,at = c(1.2,1.8),line = 10)
exp.vars <- c("TWLE",
paste0("expTWLE.", c("OLS1","OLS2","Lasso1","Lasso2","np")))
pch = c(19, 0,3,2,4,5)
ii <- 1
col = c("grey", rep("lightgrey",5))
for(vv in exp.vars){
# vv <- "TWLE"
x <- dat0[which(dat0$idschool %in% id),vv]
abline(h = c(0.92+ii*0.08,1.52+ii*0.08), lty=1+1*(ii>1),col=col[ii])
points(x=x,y=c(0.92+ii*0.08,1.52+ii*0.08),type="p",pch=pch[ii])
ii <- ii + 1
}
## End(Not run)
LSAmitR documentation built on June 1, 2022, 9:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| Kapitel 9 | R Documentation |
Kapitel 9: Fairer Vergleich in der Rueckmeldung
Description
Das ist die Nutzerseite zum Kapitel 9, Fairer Vergleich in der Rückmeldung, im Herausgeberband Large-Scale Assessment mit R: Methodische Grundlagen der österreichischen Bildungsstandardüberprüfung. Im Abschnitt Details werden die im Kapitel verwendeten R-Syntaxen zur Unterstützung für Leser/innen kommentiert und dokumentiert. Im Abschnitt Examples werden die R-Syntaxen des Kapitels vollständig wiedergegeben und gegebenenfalls erweitert.
Details
Vorbereitungen
Der zur Illustration verwendete Datensatz dat beinhaltet (imputierte)
aggregierte Leistungs- und Hintergrunddaten von 244 Schulen, bestehend aus 74
ländlichen allgemeinbildenden Pflichtschulen (APS, Stratum 1), 69 städtischen
APS (Stratum 2), 52 ländlichen allgemeinbildenden höheren Schulen (AHS, Stratum
3) und 49 städtischen AHS (Stratum 4).
Im Kapitel wird zur Bildung von Interaktionseffekten und quadratischen Termen
der Kovariaten eine neue Funktion covainteraction verwendet.
data(datenKapitel09)
dat <- datenKapitel09
covainteraction <- function(dat,covas,nchar){
for(ii in 1:(length(covas))){
vv1 <- covas[ii]
# Interaktion von vv1 mit sich selbst
subname1 <- substr(vv1,1,nchar)
intvar <- paste0(subname1, subname1)
if(vv1 == covas[1]){
dat.int <- dat[,vv1]*dat[,vv1];
newvars <- intvar } else {
dat.int <- cbind(dat.int,dat[,vv1]*dat[,vv1]);
newvars <- c(newvars,intvar)
}
# Interaktion von vv1 mit restlichen Variablen
if(ii < length(covas)){
for(jj in ((ii+1):length(covas))){
vv2 <- covas[jj]
subname2 <- substr(vv2,1,nchar)
intvar <- paste0(subname1, subname2)
newvars <- c(newvars, intvar)
dat.int <- cbind(dat.int,dat[,vv1]*dat[,vv2])
}
}
}
dat.int <- data.frame(dat.int)
names(dat.int) <- newvars
return(dat.int)
}
Abschnitt 9.2.5.1: Kovariaten und Interaktionsterme
Listing 1: Kovariatenauswahl und z-Standardisierung
Als Variablen zur Beschreibung von Kontext und Schülerzusammensetzung in den
Schulen werden in diesem Beispiel die logarithmierte Schulgröße groesse,
der Anteil an Mädchen female, der Anteil an Schülerinnen und Schülern mit
Migrationshintergrund mig und der mittlere sozioökonomische Status (SES)
sozstat eingeführt.
Die abhängige Variable ist die aggregierte
Schülerleistung der Schule TWLE. Alle Kovariaten (vars) werden
zunächst z-standardisiert (zvars).
vars <- c("groesse","female","mig","sozstat")
zvars <- paste0("z",vars)
dat[,zvars] <- scale(dat[,vars],scale = TRUE)
Listing 2: Interaktionen bilden, z-standardisieren
Zur Optimierung der Modellspezifikation werden Interaktionseffekte und
quadratische Terme der Kovariaten gebildet, dann z-standardisiert und in den
Gesamtdatensatz hinzugefügt.
Die neuen Variablennamen sind in der Liste intvars aufgelistet.
dat1 <- LSAmitR::covainteraction(dat = dat,covas = zvars,nchar = 4)
intvars <- names(dat1) # Interaktionsvariablen
dat1[,intvars] <- scale(dat1[,intvars],scale = TRUE)
dat <- cbind(dat,dat1)
Listing 3: Haupt- und Interaktionseffekte
maineff <- zvars # Haupteffekte
alleff <- c(zvars,intvars) # Haupt- und Interaktionseffekte
Abschnitt 9.2.5.2: OLS-Regression
Listing 4: OLS-Regression mit Haupteffekten
Das OLS-Regressionsmodell mit den Haupteffekten als Modellprädiktoren
(ols.mod1) für Schulen im Stratum (st) 4
fm.ols1 <- paste0("TWLE ~ ",paste(maineff,collapse=" + "))
fm.ols1 <- as.formula(fm.ols1) # Modellgleichung
st <- 4
pos <- which(dat$stratum == st) # Schulen im Stratum st
ols.mod1 <- lm(formula = fm.ols1,data = dat[pos,]) # Regression
Abschnitt 9.2.5.3: Lasso-Regression
Listing 5: Datenaufbereitung
Für die Durchführung der Lasso-Regression wird das R-Paket glmnet
(Friedman et al., 2010) eingesetzt. Die Kovariatenmatrix (Z) sowie der
Vektor der abhängigen Leistungswerte (Y) müssen definiert werden.
library(glmnet)
Z <- as.matrix(dat[pos,alleff]) # Kovariatenmatrix
Y <- dat$TWLE[pos] # Abhängige Variable
Listing 6: Bestimmung Teilmengen für Kreuzvalidierung, Lasso-Regression
Das Lasso-Verfahren wird mit der Funktion cv.glmnet() durchgeführt.
Zur Auswahl eines optimalen shrinkage λ wird das Verfahren der
K-fachen Kreuzvalidierung verwendet.
Die Schulstichprobe wird durch ID-Zuweisung (foldid) verschiedenen
Teilmengen zugewiesen.
nid <- floor(length(pos)/3) # Teilmengen definieren
foldid <- rep(c(1:nid),3,length.out=length(pos)) # Zuweisung
lasso.mod2 <- glmnet::cv.glmnet(x=Z,y=Y,alpha = 1, foldid = foldid)
Listing 7: Erwartungswerte der Schulen
Entsprechend lambda.min werden die Regressionskoeffizienten und die
entsprechenden Erwartungswerte der Schulen bestimmt.
lasso.pred2 <- predict(lasso.mod2,newx = Z,s="lambda.min")
dat$expTWLE.Lasso2[pos] <- as.vector(lasso.pred2)
Listing 8: Bestimmung R^2
Bestimmung des aufgeklärten Varianzanteils der Schulleistung R^2.
varY <- var(dat$TWLE[pos])
varY.lasso.mod2 <- var(dat$expTWLE.Lasso2[pos])
R2.lasso.mod2 <- varY.lasso.mod2/varY
Abschnitt 9.2.5.4: Nichtparametrische Regression
Listing 9: Distanzberechnung
Der erste Schritt zur Durchführung einer nichtparametrischen Regression ist die
Erstellung der Distanzmatrix zwischen Schulen. In diesem Beispiel wird die
euklidische Distanz als Distanzmaß berechnet, alle standardisierten Haupteffekte
sind eingeschlossen. Außerdem setzen wir die Gewichte von allen Kovariaten
(gi) auf 1. dfr.i in diesem Beispiel ist die Distanzmatrix für
erste Schule im Stratum.
N <- length(pos) # Anzahl Schulen im Stratum
schools <- dat$idschool[pos] # Schulen-ID
i <- 1
# Teildatensatz von Schule i
dat.i <- dat[pos[i],c("idschool","TWLE",maineff)]
names(dat.i) <- paste0(names(dat.i),".i")
# Daten der Vergleichsschulen
dat.vgl <- dat[pos[-i],c("idschool","TWLE",maineff)]
index.vgl <- match(dat.vgl$idschool,schools)
# Daten zusammenfügen
dfr.i <- data.frame("index.i"=i,dat.i,"index.vgl"=index.vgl,
dat.vgl, row.names=NULL)
# Distanz zur Schule i
dfr.i$dist <- 0
gi <- c(1,1,1,1)
for(ii in 1:length(maineff)){
vv <- maineff[ii]
pair.vv <- grep(vv, names(dfr.i), value=T)
dist.vv <- gi[ii]*((dfr.i[,pair.vv[1]]-dfr.i[,pair.vv[2]])^2)
dfr.i$dist <- dfr.i$dist + dist.vv }
Listing 10: H initiieren
p(x) = λ^x exp(-λ)/x!.
Die Gewichte w_{ik} für jedes Paar (i, k) von Schulen werden mithilfe der
Distanz, der Gauß’schen Kernfunktion (dnorm) als Transformationsfunktion
sowie einer schulspezifischen Bandweite h_i berechnet. Die Auswahl des optimalen
Werts \hat{h_i} für jede Schule i erfolgt nach Vieu (1991). Zunächst wird
ein Vektor H so gewählt, dass der optimale Wert \hat{h_i} größer
als der kleinste und kleiner als der größte Wert in H ausfällt.
Je kleiner das Intervall zwischen den Werten in H ist, desto
wahrscheinlicher ist, dass ein Listenelement den optimalen Wert erlangt.
Auf der anderen Seite korrespondiert die Rechenzeit mit der Länge von H.
Gemäß der Größe der Vergleichsgruppe wählen wir eine Länge von 30 für H,
zusätzlich wird ein sehr großer Wert (100000) für die Fälle hinzugefügt, bei
denen alle Gewichte beinahe gleich sind.
d.dist <- max(dfr.i$dist)-min(dfr.i$dist)
H <- c(seq(d.dist/100,d.dist,length=30),100000)
V1 <- length(H)
# Anzahl Vergleichsschulen
n <- nrow(dfr.i)
Listing 11: Leave-One-Out-Schätzer der jeweiligen Vergleichsschule k nach h in H
Auf Basis aller Werte in H und dem jeweils entsprechenden Gewicht w_{ik}
(wgt.h) werden die Leave-One-Out-Schätzer der jeweiligen Vergleichsschule
(pred.k) berechnet.
sumw <- 0*H # Vektor w_{ik} initiieren, h in H
av <- "TWLE"
dfr0.i <- dfr.i[,c("idschool",av)]
# Schleife über alle h-Werte
for (ll in 1:V1 ){
h <- H[ll]
# Gewicht w_{ik} bei h
dfr.i$wgt.h <- dnorm(sqrt(dfr.i$dist), mean=0, sd=sqrt(h))
# Summe von w_{ik} bei h
sumw[which(H==h)] <- sum(dfr.i$wgt.h)
# Leave-one-out-Schätzer von Y_k
for (k in 1:n){
# Regressionsformel
fm <- paste0(av,"~",paste0(maineff,collapse="+"))
fm <- as.formula(fm)
# Regressionsanalyse ohne Beitrag von Schule k
dfr.i0 <- dfr.i[-k,]
mod.k <- lm(formula=fm,data=dfr.i0,weights=dfr.i0$wgt.h)
# Erwartungswert anhand Kovariaten der Schule k berechnen
pred.k <- predict(mod.k, dfr.i)[k]
dfr0.i[k,paste0( "h_",h) ] <- pred.k
}}
# Erwartungswerte auf Basis verschiedener h-Werte
dfr1 <- data.frame("idschool.i"=dfr.i$idschool.i[1],"h"=H )
Listing 12: Kreuzvalidierungskriterium nach h in H
Zur Berechnung der Kreuzvalidierungskriterien (CVh, je kleiner, desto
reliabler sind die Schätzwerte) für alle Werte h in H verwenden wir in
diesem Beispiel die Plug-in-Bandweite nach Altman und Leger (1995) (hAL),
die mit der Funktion ALbw() des R-Pakets kerdiest aufrufbar ist.
library(kerdiest)
hAL <- kerdiest::ALbw("n",dfr.i$dist) # Plug-in Bandweite
dfr.i$cross.h <- hAL
dfr.i$crosswgt <- dnorm( sqrt(dfr.i$dist), mean=0, sd = sqrt(hAL) )
# Kreuzvalidierungskriterium CVh
vh <- grep("h_",colnames(dfr0.i),value=TRUE)
for (ll in 1:V1){
dfr1[ll,"CVh"] <- sum( (dfr0.i[,av] - dfr0.i[,vh[ll]])^2 *
dfr.i$crosswgt) / n}
Listing 13: Bestimmung des optimalen Wertes von h
Der optimale Wert von h in H (h.min) entspricht dem mit dem
kleinsten resultierenden CVh.
dfr1$min.h.index <- 0
ind <- which.min( dfr1$CVh )
dfr1$min.h.index[ind] <- 1
dfr1$h.min <- dfr1$h[ind]
Listing 14: Kleinste Quadratsumme der Schätzfehler
Kleinste Quadratsumme der Schätzfehler der nichtparametrischen Regression mit
h=h.min.
dfr1$CVhmin <- dfr1[ ind , "CVh" ]
Listing 15: Effizienzsteigerung
Die Effizienz (Steigerung der Schätzungsreliabilität) der nichtparametrischen
Regression gegenüber der linearen Regression (äquivalent zu CVh bei
h=100000).
dfr1$eff_gain <- 100 * ( dfr1[V1,"CVh"] / dfr1$CVhmin[1] - 1 )
Listing 16: Durchführung der nichtparametrischen Regression
h <- dfr1$h.min[1] # h.min
dfr.i$wgt.h <- dnorm(sqrt(dfr.i$dist),sd=sqrt(h))/
dnorm(0,sd= sqrt(h)) # w_{ik} bei h.min
dfr.i0 <- dfr.i
# Lokale Regression
mod.ii <- lm(formula=fm,data=dfr.i0,weights=dfr.i0$wgt.h)
# Kovariaten Schule i
predM <- data.frame(dfr.i[1,paste0(maineff,".i")])
names(predM) <- maineff
pred.ii <- predict(mod.ii, predM) # Schätzwert Schule i
dat$expTWLE.np[match(dfr1$idschool.i[1],dat$idschool)] <- pred.ii
Abschnitt 9.2.7, Berücksichtigung der Schätzfehler
Der Erwartungsbereich wird nach der im Buch beschriebenen Vorgehensweise bestimmt.
Listing 17: Bestimmung des Erwartungsbereichs
Bestimmung der Breite des Erwartungsbereichs aller Schulen auf Basis der Ergebnisse der OLS-Regression mit Haupteffekten.
vv <- "expTWLE.OLS1" # Variablenname
mm <- "OLS1" # Kurzname des Modells
dfr <- NULL # Ergebnistabelle
# Schleife über alle möglichen Breite von 10 bis 60
for(w in 10:60){
# Variablen für Ergebnisse pro w
var <- paste0(mm,".pos.eb",w) # Position der Schule
var.low <- paste0(mm,".eblow",w) # Untere Grenze des EBs
var.upp <- paste0(mm,".ebupp",w) # Obere Grenze des EBs
# Berechnen
dat[,var.low] <- dat[,vv]-w/2 # Untere Grenze des EBs
dat[,var.upp] <- dat[,vv]+w/2 # Obere Grenze des EBs
# Position: -1=unterhalb, 0=innerhalb, 1=oberhalb des EBs
dat[,var] <- -1*(dat$TWLE < dat[,var.low]) + 1*(dat$TWLE > dat[,var.upp])
# Verteilung der Schulpositionen
tmp <- data.frame(t(matrix(prop.table(table(dat[,var])))))
names(tmp) <- c("unterhalb","innerhalb","oberhalb")
tmp <- data.frame("ModellxBereich"=var,tmp); dfr <- rbind(dfr,tmp) }
# Abweichung zur Wunschverteilung 25-50-25
dfr1 <- dfr
dfr1[,c(2,4)] <- (dfr1[,c(2,4)] - .25)^2
dfr1[,3] <- (dfr1[,3] - .5)^2
dfr1$sumquare <- rowSums(dfr1[,-1])
# Auswahl markieren
dfr$Auswahl <- 1*(dfr1$sumquare == min(dfr1$sumquare) )
Author(s)
Giang Pham, Alexander Robitzsch, Ann Cathrice George, Roman Freunberger
References
Pham, G., Robitzsch, A., George, A. C. & Freunberger, R. (2016). Fairer Vergleich in der Rückmeldung. In S. Breit & C. Schreiner (Hrsg.), Large-Scale Assessment mit R: Methodische Grundlagen der österreichischen Bildungsstandardüberprüfung (pp. 295–332). Wien: facultas.
See Also
Zu datenKapitel09, den im Kapitel verwendeten Daten.
Zurück zu Kapitel 8, Fehlende Daten und Plausible Values.
Zu Kapitel 10, Reporting und Analysen.
Zur Übersicht.
Examples
## Not run:
library(miceadds)
library(glmnet)
library(kerdiest)
covainteraction <- function(dat,covas,nchar){
for(ii in 1:(length(covas))){
vv1 <- covas[ii]
# Interaktion von vv1 mit sich selbst
subname1 <- substr(vv1,1,nchar)
intvar <- paste0(subname1, subname1)
if(vv1 == covas[1]){
dat.int <- dat[,vv1]*dat[,vv1];
newvars <- intvar } else {
dat.int <- cbind(dat.int,dat[,vv1]*dat[,vv1]);
newvars <- c(newvars,intvar)
}
# Interaktion von vv1 mit restlichen Variablen
if(ii < length(covas)){
for(jj in ((ii+1):length(covas))){
vv2 <- covas[jj]
subname2 <- substr(vv2,1,nchar)
intvar <- paste0(subname1, subname2)
newvars <- c(newvars, intvar)
dat.int <- cbind(dat.int,dat[,vv1]*dat[,vv2])
}
}
}
dat.int <- data.frame(dat.int)
names(dat.int) <- newvars
return(dat.int)
}
data(datenKapitel09)
dat <- datenKapitel09
# Platzhalter für Leistungsschätzwerte aller Modelle
dat$expTWLE.OLS1 <- NA
dat$expTWLE.OLS2 <- NA
dat$expTWLE.Lasso1 <- NA
dat$expTWLE.Lasso2 <- NA
dat$expTWLE.np <- NA
## -------------------------------------------------------------
## Abschnitt 9.2.5, Umsetzung in R
## -------------------------------------------------------------
# -------------------------------------------------------------
# Abschnitt 9.2.5.1, Listing 1: Kovariatenauswahl und
# z-Standardisierung
#
vars <- c("groesse","female","mig","sozstat")
zvars <- paste0("z",vars)
dat[,zvars] <- scale(dat[,vars],scale = TRUE)
# -------------------------------------------------------------
# Abschnitt 9.2.5.1, Listing 2:
#
# Interaktionen bilden, z-standardisieren
dat1 <- LSAmitR::covainteraction(dat = dat,covas = zvars,nchar = 4)
intvars <- names(dat1) # Interaktionsvariablen
dat1[,intvars] <- scale(dat1[,intvars],scale = TRUE)
dat <- cbind(dat,dat1)
# -------------------------------------------------------------
# Abschnitt 9.2.5.1, Listing 3: Modellprädiktoren: Haupt- und
# Interaktionseffekte
#
maineff <- zvars # Haupteffekte
alleff <- c(zvars,intvars) # Haupt- und Interaktionseffekte
# -------------------------------------------------------------
# Abschnitt 9.2.5.2, Listing 4: OLS-Regression mit Haupteffekten
#
fm.ols1 <- paste0("TWLE ~ ",paste(maineff,collapse=" + "))
fm.ols1 <- as.formula(fm.ols1) # Modellgleichung
st <- 4
pos <- which(dat$stratum == st) # Schulen im Stratum st
ols.mod1 <- lm(formula = fm.ols1,data = dat[pos,]) # Regression
# -------------------------------------------------------------
# Abschnitt 9.2.5.3, Listing 5: Lasso-Regression
# Datenaufbereitung
#
library(glmnet)
Z <- as.matrix(dat[pos,alleff]) # Kovariatenmatrix
Y <- dat$TWLE[pos] # Abhängige Variable
# -------------------------------------------------------------
# Abschnitt 9.2.5.3, Listing 6: Lasso-Regression
# Bestimmung Teilmengen für Kreuzvalidierung, Lasso-Regression
#
nid <- floor(length(pos)/3) # Teilmengen definieren
foldid <- rep(c(1:nid),3,length.out=length(pos)) # Zuweisung
lasso.mod2 <- glmnet::cv.glmnet(x=Z,y=Y,alpha = 1, foldid = foldid)
# -------------------------------------------------------------
# Abschnitt 9.2.5.3, Listing 7: Lasso-Regression
# Erwartungswerte der Schulen
#
lasso.pred2 <- predict(lasso.mod2,newx = Z,s="lambda.min")
dat$expTWLE.Lasso2[pos] <- as.vector(lasso.pred2)
# -------------------------------------------------------------
# Abschnitt 9.2.5.3, Listing 8: Lasso-Regression
# Bestimmung R^2
#
varY <- var(dat$TWLE[pos])
varY.lasso.mod2 <- var(dat$expTWLE.Lasso2[pos])
R2.lasso.mod2 <- varY.lasso.mod2/varY
# -------------------------------------------------------------
# Abschnitt 9.2.5.4, Listing 9: Nichtparametrische Regression
# Distanzberechnung zur Schule i (Stratum st)
#
N <- length(pos) # Anzahl Schulen im Stratum
schools <- dat$idschool[pos] # Schulen-ID
i <- 1
# Teildatensatz von Schule i
dat.i <- dat[pos[i],c("idschool","TWLE",maineff)]
names(dat.i) <- paste0(names(dat.i),".i")
# Daten der Vergleichsschulen
dat.vgl <- dat[pos[-i],c("idschool","TWLE",maineff)]
index.vgl <- match(dat.vgl$idschool,schools)
# Daten zusammenfügen
dfr.i <- data.frame("index.i"=i,dat.i,"index.vgl"=index.vgl,
dat.vgl, row.names=NULL)
# Distanz zur Schule i
dfr.i$dist <- 0
gi <- c(1,1,1,1)
for(ii in 1:length(maineff)){
vv <- maineff[ii]
pair.vv <- grep(vv, names(dfr.i), value=T)
dist.vv <- gi[ii]*((dfr.i[,pair.vv[1]]-dfr.i[,pair.vv[2]])^2)
dfr.i$dist <- dfr.i$dist + dist.vv }
# -------------------------------------------------------------
# Abschnitt 9.2.5.4, Listing 10: Nichtparametrische Regression
#
# H initiieren
d.dist <- max(dfr.i$dist)-min(dfr.i$dist)
H <- c(seq(d.dist/100,d.dist,length=30),100000)
V1 <- length(H)
# Anzahl Vergleichsschulen
n <- nrow(dfr.i)
# -------------------------------------------------------------
# Abschnitt 9.2.5.4, Listing 11: Nichtparametrische Regression
# Berechnung der Leave-One-Out-Schätzer der jeweiligen
# Vergleichsschule k nach h in H
#
sumw <- 0*H # Vektor w_{ik} initiieren, h in H
av <- "TWLE"
dfr0.i <- dfr.i[,c("idschool",av)]
# Schleife über alle h-Werte
for (ll in 1:V1 ){
h <- H[ll]
# Gewicht w_{ik} bei h
dfr.i$wgt.h <- dnorm(sqrt(dfr.i$dist), mean=0, sd=sqrt(h))
# Summe von w_{ik} bei h
sumw[which(H==h)] <- sum(dfr.i$wgt.h)
# Leave-one-out-Schätzer von Y_k
for (k in 1:n){
# Regressionsformel
fm <- paste0(av,"~",paste0(maineff,collapse="+"))
fm <- as.formula(fm)
# Regressionsanalyse ohne Beitrag von Schule k
dfr.i0 <- dfr.i[-k,]
mod.k <- lm(formula=fm,data=dfr.i0,weights=dfr.i0$wgt.h)
# Erwartungswert anhand Kovariaten der Schule k berechnen
pred.k <- predict(mod.k, dfr.i)[k]
dfr0.i[k,paste0( "h_",h) ] <- pred.k
}}
# Erwartungswerte auf Basis verschiedener h-Werte
dfr1 <- data.frame("idschool.i"=dfr.i$idschool.i[1],"h"=H )
# -------------------------------------------------------------
# Abschnitt 9.2.5.4, Listing 12: Nichtparametrische Regression
# Berechnung des Kreuzvalidierungskriteriums nach h in H
#
library(kerdiest)
hAL <- kerdiest::ALbw("n",dfr.i$dist) # Plug-in Bandweite
dfr.i$cross.h <- hAL
dfr.i$crosswgt <- dnorm( sqrt(dfr.i$dist), mean=0, sd = sqrt(hAL) )
# Kreuzvalidierungskriterium CVh
vh <- grep("h_",colnames(dfr0.i),value=TRUE)
for (ll in 1:V1){
dfr1[ll,"CVh"] <- sum( (dfr0.i[,av] - dfr0.i[,vh[ll]])^2 *
dfr.i$crosswgt) / n}
# -------------------------------------------------------------
# Abschnitt 9.2.5.4, Listing 13: Nichtparametrische Regression
# Bestimmung optimales Wertes von h (h.min)
#
dfr1$min.h.index <- 0
ind <- which.min( dfr1$CVh )
dfr1$min.h.index[ind ] <- 1
dfr1$h.min <- dfr1$h[ind]
# -------------------------------------------------------------
# Abschnitt 9.2.5.4, Listing 14: Nichtparametrische Regression
# Kleinste Quadratsumme der Schätzfehler
#
dfr1$CVhmin <- dfr1[ ind , "CVh" ]
# -------------------------------------------------------------
# Abschnitt 9.2.5.4, Listing 15: Nichtparametrische Regression
# Effizienzsteigerung berechnen
#
dfr1$eff_gain <- 100 * ( dfr1[V1,"CVh"] / dfr1$CVhmin[1] - 1 )
# -------------------------------------------------------------
# Abschnitt 9.2.5.4, Listing 16: Nichtparametrische Regression
# Durchführung der nichtparametrischen Regression bei h=h.min
#
h <- dfr1$h.min[1] # h.min
dfr.i$wgt.h <- dnorm(sqrt(dfr.i$dist),sd=sqrt(h))/
dnorm(0,sd= sqrt(h)) # w_{ik} bei h.min
dfr.i0 <- dfr.i
# Lokale Regression
mod.ii <- lm(formula=fm,data=dfr.i0,weights=dfr.i0$wgt.h)
# Kovariaten Schule i
predM <- data.frame(dfr.i[1,paste0(maineff,".i")])
names(predM) <- maineff
pred.ii <- predict(mod.ii, predM) # Schätzwert Schule i
dat[match(dfr1$idschool.i[1],dat$idschool), "expTWLE.np"] <- pred.ii
## -------------------------------------------------------------
## Abschnitt 9.2.5, Umsetzung in R, Ergänzung zum Buch
## -------------------------------------------------------------
# Korrelationen zwischen Haupteffekten
cor(dat[,maineff]) # gesamt
# Pro Stratum
for(s in 1:4) print(cor(dat[which(dat$stratum == s),maineff]))
# -------------------------------------------------------------
# Abschnitt 9.2.5.2, Ergänzung zum Buch
# OLS-Regression
#
# Modellgleichung nur mit Haupteffekten
fm.ols1 <- paste0("TWLE ~ ",paste(maineff,collapse=" + "))
fm.ols1 <- as.formula(fm.ols1)
# Modellgleichung mit Haupteffekten ohne zgroesse
fm.ols1a <- paste0("TWLE ~ ",paste(setdiff(maineff,c("zgroesse")),
collapse=" + "))
fm.ols1a <- as.formula(fm.ols1a)
# Modellgleichung mit Haupt- und Interaktionseffekten
fm.ols2 <- paste0("TWLE ~ ",paste(alleff,collapse=" + "))
fm.ols2 <- as.formula(fm.ols2)
# Ergebnistabelle über 4 Strata hinweg vorbereiten
tab1 <- data.frame("Variable"=c("(Intercept)",maineff))
tab2 <- data.frame("Variable"=c("(Intercept)",alleff))
# Durchführung: Schleife über vier Strata
for(st in 1:4){
# st <- 4
# Position Schulen des Stratums st im Datensatz
pos <- which(dat$stratum == st)
#---------------------------------
# OLS-Modell 1
# Durchführung
ols.mod1 <- lm(formula = fm.ols1,data = dat[pos,])
ols.mod1a <- lm(formula = fm.ols1a,data = dat[pos,])
# Modellergebnisse anzeigen
summary(ols.mod1)
summary(ols.mod1a)
# Erwartungswerte der Schulen
dat$expTWLE.OLS1[pos] <- fitted(ols.mod1)
# Ergebnisse in Tabelle speichern
par <- summary(ols.mod1)
tab.s <- data.frame(par$coef,R2=par$r.squared,R2.adj=par$adj.r.squared)
names(tab.s) <- paste0("stratum",st,
c("_coef","_SE","_t","_p","_R2","_R2.adj"))
tab1 <- cbind(tab1, tab.s)
# Durchführung OLS-Modell 2
ols.mod2 <- lm(formula = fm.ols2,data = dat[pos,])
# Modellergebnisse anzeigen
summary(ols.mod2)
# Erwartungswerte der Schulen
dat$expTWLE.OLS2[pos] <- fitted(ols.mod2)
# Ergebnisse in Tabelle speichern
par <- summary(ols.mod2)
tab.s <- data.frame(par$coef,R2=par$r.squared,R2.adj=par$adj.r.squared)
names(tab.s) <- paste0("stratum",st,
c("_coef","_SE","_t","_p","_R2","_R2.adj"))
tab2 <- cbind(tab2, tab.s)
}
# Daten Schule 1196 ansehen
dat[which(dat$idschool == 1196),]
# Schätzwerte nach ols.mod1 und ols.mod2 vergleichen
summary(abs(dat$expTWLE.OLS1 - dat$expTWLE.OLS2))
cor.test(dat$expTWLE.OLS1,dat$expTWLE.OLS2)
# Grafische Darstellung des Vergleich (Schule 1196 rot markiert)
plot(dat$expTWLE.OLS1,dat$expTWLE.OLS2,xlim=c(380,650),ylim=c(380,650),
col=1*(dat$idschool == 1196)+1,pch=15*(dat$idschool == 1196)+1)
abline(a=0,b=1)
# -------------------------------------------------------------
# Abschnitt 9.2.5.3, Ergänzung zum Buch
# Lasso-Regression
#
library(glmnet)
# Variablen für Erwartungswerte
dat$expTWLE.Lasso2 <- dat$expTWLE.Lasso1 <- NA
# Tabelle für Modellergebnisse
tab3 <- data.frame("Variable"=c("(Intercept)",maineff))
tab4 <- data.frame("Variable"=c("(Intercept)",alleff))
for(st in 1:4){
# st <- 4
# Position Schulen des Stratums st im Datensatz
pos <- which(dat$stratum == st)
#------------------------------------------------------------#
# Lasso-Regression mit den Haupteffekten
# Kovariatenmatrix
Z <- as.matrix(dat[pos,maineff])
# Abhängige Variable
Y <- dat$TWLE[pos]
# Kreuzvalidierung: Teilmengen definieren
nid <- floor(length(pos)/3)
# Schulen zu Teilmengen zuordnen
foldid <- rep(c(1:nid),3,length.out=length(pos))
# Regression
lasso.mod1 <- cv.glmnet(x=Z,y=Y,alpha = 1, foldid = foldid)
# Ergebnisse ansehen
print(lasso.mod1)
# Lasso-Koeffizienten bei lambda.min
print(lasso.beta <- coef(lasso.mod1,s="lambda.min"))
# Erwartungswerte der Schulen
lasso.pred1 <- predict(lasso.mod1,newx = Z,s="lambda.min")
dat$expTWLE.Lasso1[pos] <- as.vector(lasso.pred1)
# R2 bestimmen
varY <- var(dat$TWLE[pos])
varY.lasso.mod1 <- var(dat$expTWLE.Lasso1[pos])
print(R2.lasso.mod1 <- varY.lasso.mod1/varY)
# Ergebnistabelle
vv <- paste0("coef.stratum",st); tab3[,vv] <- NA
tab3[lasso.beta@i+1,vv] <- lasso.beta@x
vv <- paste0("lambda.stratum",st); tab3[,vv] <- lasso.mod1$lambda.min
vv <- paste0("R2.stratum",st); tab3[,vv] <- R2.lasso.mod1
#------------------------------------------------------------#
# Lasso-Regression mit Haupt- und Interaktionseffekten
# Kovariatenmatrix
Z <- as.matrix(dat[pos,alleff])
# Regression
lasso.mod2 <- cv.glmnet(x=Z,y=Y,alpha = 1, foldid = foldid)
# Ergebnisausdruck
print(lasso.mod2)
# Lasso-Koeffizienten bei lambda.min
print(lasso.beta <- coef(lasso.mod2,s="lambda.min"))
# Erwartungswerte der Schulen
lasso.pred2 <- predict(lasso.mod2,newx = Z,s="lambda.min")
dat$expTWLE.Lasso2[pos] <- as.vector(lasso.pred2)
# R2 bestimmen
varY.lasso.mod2 <- var(dat$expTWLE.Lasso2[pos])
R2.lasso.mod2 <- varY.lasso.mod2/varY
R2.lasso.mod2
# Ergebnistabelle
vv <- paste0("coef.stratum",st); tab4[,vv] <- NA
tab4[lasso.beta@i+1,vv] <- lasso.beta@x
vv <- paste0("lambda.stratum",st); tab4[,vv] <- lasso.mod2$lambda.min
vv <- paste0("R2.stratum",st); tab4[,vv] <- R2.lasso.mod2
}
# Regressionresiduen = Schätzung von SChul- und Unterrichtseffekt
dat$resTWLE.Lasso1 <- dat$TWLE - dat$expTWLE.Lasso1
dat$resTWLE.Lasso2 <- dat$TWLE - dat$expTWLE.Lasso2
# -------------------------------------------------------------
# Abschnitt 9.2.5.4, Ergänzung zum Buch
# Nichtparametrische Regression
#
#
# Achtung: Der nachfolgende Algorithmus benötigt viel Zeit!
#
av <- "TWLE" # Abhängige Variable
dfr3 <- NULL # Ergebnistabelle
# Variable für Leistungsschätzwerte
# Schleife über 4 Strata
for(st in 1:4){
# st <- 1
pos <- which(dat$stratum == st)
N <- length(pos)
schools <- dat$idschool[pos]
###
# Distanzmatrix dfr für alle Schulen im Stratum erstellen
dfr <- NULL
for (i in 1:N){
# i <- 1
# Teildatensatz von Schule i
dat.i <- dat[pos[i],c("idschool","TWLE",maineff)]
# Daten der Vergleichsgruppe
dat.vgl <- dat[pos[-i],c("idschool","TWLE",maineff)]
# Variablennamen von dat.vgl umbenennen
# names(dat.vgl) <- paste0("vgl.",names(dat.vgl))
# Variablennamen von dat.i umbenennen
names(dat.i) <- paste0(names(dat.i),".i")
# Daten zusammenfügen
index.vgl <- match(dat.vgl$idschool,schools)
dfr.i <- data.frame("index.i"=i,dat.i,
"index.vgl"=index.vgl,dat.vgl,
row.names=NULL)
# Distanz zur i
dfr.i$dist <- 0
gi <- c(1,1,1,1)
for(ii in 1:length(maineff)){
vv <- maineff[ii]
pair.vv <- grep(vv, names(dfr.i), value=T)
dist.vv <- gi[ii]*((dfr.i[,pair.vv[1]]-dfr.i[,pair.vv[2]])^2)
dfr.i$dist <- dfr.i$dist + dist.vv
}
print(i) ; flush.console()
dfr <- rbind( dfr , dfr.i )
}
dfr1 <- index.dataframe( dfr , systime=TRUE )
###
# h-Auswahl und Nichtparametrische Regression pro Schule i
dfr1.list <- list()
for (i in 1:N){
# i <- 1
dfr.i <- dfr[ dfr$index.i == i , ]
n <- nrow(dfr.i)
# Startwertliste für h initiieren
d.dist <- max(dfr.i$dist)-min(dfr.i$dist)
H <- c(seq(d.dist/100,d.dist,length=30),100000)
V1 <- length(H) # Anzahl der Startwerte in H
# Startwerte: Summe von w_ik
sumw <- 0*H
dfr0.i <- dfr.i[,c("idschool",av)]
# Schleife über alle h-Werte
for (ll in 1:V1 ){
h <- H[ll]
# Gewicht w_ik bei h
dfr.i$wgt.h <- dnorm(sqrt(dfr.i$dist), mean=0, sd=sqrt(h))
# Summe von w_ik bei h
sumw[which(H==h)] <- sum(dfr.i$wgt.h)
# Leave-one-out-Schätzer von Y_k
for (k in 1:n){
# Regressionsformel
fm <- paste0(av,"~",paste0(maineff,collapse="+"))
fm <- as.formula(fm)
# Regressionsanalyse ohne Beitrag von Schule k
dfr.i0 <- dfr.i[-k,]
mod.k <- lm(formula=fm,data=dfr.i0,weights=dfr.i0$wgt.h)
# Erwartungswert anhand Kovariaten der Schule k berechnen
pred.k <- predict(mod.k, dfr.i)[k]
dfr0.i[k,paste0( "h_",h) ] <- pred.k
}
print(paste0("i=",i,", h_",ll))
}
# Erwartungswerte auf Basis verschiedener h-Werte
dfr1 <- data.frame("idschool.i"=dfr.i$idschool.i[1],"h"=H )
# Berechnung des Kreuzvalidierungskriteriums
library(kerdiest)
hAL <- kerdiest::ALbw("n",dfr.i$dist) # Plug-in Bandbreite nach Altman und
# Leger
name <- paste0( "bandwidth_choice_school" , dfr.i$idschool.i[1] ,
"_cross.h_" , round2(hAL,1) )
# Regressionsgewichte auf Basis cross.h
dfr.i$cross.h <- hAL
dfr.i$crosswgt <- dnorm( sqrt(dfr.i$dist), mean=0, sd = sqrt(hAL) )
dfr.i <- index.dataframe( dfr.i , systime=TRUE )
# Kreuzvalidierungskriterium CVh
vh <- grep("h_",colnames(dfr0.i),value=TRUE)
for (ll in 1:V1){
# ll <- 5
dfr1[ll,"CVh"] <- sum( (dfr0.i[,av] - dfr0.i[,vh[ll]])^2 *
dfr.i$crosswgt) / n
print(ll)
}
# Bestimmung h.min
dfr1$min.h.index <- 0
ind <- which.min( dfr1$CVh )
dfr1$min.h.index[ind ] <- 1
dfr1$h.min <- dfr1$h[ind]
# Kleinste Quadratsumme der Schätzfehler
dfr1$CVhmin <- dfr1[ ind , "CVh" ]
# Effizienzsteigerung berechnen
dfr1$eff_gain <- 100 * ( dfr1[V1,"CVh"] / dfr1$CVhmin[1] - 1 )
# h auswählen
h <- dfr1$h.min[1]
# Gewichte anhand h berechnen
dfr.i$wgt.h <- dnorm( sqrt( dfr.i$dist ) , sd = sqrt( h) ) /
dnorm( 0 , sd = sqrt( h) )
dfr.i0 <- dfr.i
mod.ii <- lm(formula = fm,data = dfr.i0,weights = dfr.i0$wgt.h)
# Leistungsschätzwerte berechnen
predM <- data.frame(dfr.i[1,paste0(maineff,".i")])
names(predM) <- maineff
pred.ii <- predict( mod.ii , predM )
dfr1$fitted_np <- pred.ii
dfr1$h.min_sumwgt <- sum( dfr.i0$wgt.h )
dfr1$h_sumwgt <- sumw
# Leistungsschätzwerte zum Datensatz hinzufügen
dat$expTWLE.np[match(dfr1$idschool.i[1],dat$idschool)] <- pred.ii
dfr1.list[[i]] <- dfr1
}
###
# Ergebnisse im Stratum st zusammenfassen
dfr2 <- NULL
for(i in 1:length(dfr1.list)){
dat.ff <- dfr1.list[[i]]
dfr2.ff <- dat.ff[1,c("idschool.i","h.min","fitted_np","h.min_sumwgt",
"CVhmin","eff_gain")]
dfr2.ff$CVlinreg <- dat.ff[V1,"CVh"]
names(dfr2.ff) <- c("idschool","h.min","fitted_np","h.min_sumwgt",
"CVhmin","eff_gain","CVlinreg")
dfr2 <- rbind(dfr2, dfr2.ff)
print(i)
}
#---------------------------------------------------##
# R2 berechnen
varY <- var(dat$TWLE[pos])
varY.np <- var(dat$expTWLE.np[pos])
dfr2$R2.np <- varY.np/varY
#---------------------------------------------------##
# Zur Gesamtergebnistabelle
dfr3 <- rbind(dfr3,cbind("Stratum"=st,dfr2))
}
# Effizienz der NP-Regression gegenüber OLS-Regression
summary(dfr3$eff_gain)
table(dfr3$eff_gain > 5)
table(dfr3$eff_gain > 10)
table(dfr3$eff_gain > 20)
# Regressionsresiduen
dat$resTWLE.np <- dat$TWLE - dat$expTWLE.np
## -------------------------------------------------------------
## Abschnitt 9.2.6, Ergänzung zum Buch
## Ergebnisse im Vergleich
## -------------------------------------------------------------
# Output-Variablen
out <- grep("expTWLE",names(dat),value=T)
lt <- length(out)
# Korrelationsmatrix
tab <- tab1 <- as.matrix(round2(cor(dat[,out]),3))
# Varianzmatrix
tab2 <- as.matrix(round2(sqrt(var(dat[,out])),1))
tab3 <- matrix(NA,lt,lt)
# Differenzmatrix
for(ii in 1:(lt-1))
for(jj in (ii+1):lt) tab3[ii,jj] <- round2(mean(abs(dat[,out[jj]] -
dat[,out[ii]])),1)
tab4 <- matrix(NA,lt,lt)
# Differenzmatrix
for(ii in 1:(lt-1))
for(jj in (ii+1):lt) tab4[ii,jj] <- round2(sd(abs(dat[,out[jj]] -
dat[,out[ii]])),1)
# Ergebnistabelle
diag(tab) <- diag(tab2)
tab[upper.tri(tab)] <- tab3[upper.tri(tab3)]
# R2 Gesamt
varY <- var(dat$TWLE)
varexp.OLS1 <- var(dat$expTWLE.OLS1); R2.OLS1 <- varexp.OLS1/varY
varexp.OLS2 <- var(dat$expTWLE.OLS2); R2.OLS2 <- varexp.OLS2/varY
varexp.Lasso1 <- var(dat$expTWLE.Lasso1); R2.Lasso1 <- varexp.Lasso1/varY
varexp.Lasso2 <- var(dat$expTWLE.Lasso2); R2.Lasso2 <- varexp.Lasso2/varY
varexp.np <- var(dat$expTWLE.np); R2.np <- varexp.np/varY
R2 <- c(R2.OLS1,R2.OLS2,R2.Lasso1,R2.Lasso2,R2.np)
tab <- cbind(tab,R2)
# R2 pro Stratum
dat0 <- dat
for(st in 1:4){
# st <- 1
dat <- dat0[which(dat0$stratum == st),]
varY <- var(dat$TWLE)
varexp.OLS1 <- var(dat$expTWLE.OLS1); R2.OLS1 <- varexp.OLS1/varY
varexp.OLS2 <- var(dat$expTWLE.OLS2); R2.OLS2 <- varexp.OLS2/varY
varexp.Lasso1 <- var(dat$expTWLE.Lasso1); R2.Lasso1 <- varexp.Lasso1/varY
varexp.Lasso2 <- var(dat$expTWLE.Lasso2); R2.Lasso2 <- varexp.Lasso2/varY
varexp.np <- var(dat$expTWLE.np); R2.np <- varexp.np/varY
R2 <- c(R2.OLS1,R2.OLS2,R2.Lasso1,R2.Lasso2,R2.np)
tab <- cbind(tab,R2)
}
colnames(tab)[7:10] <- paste0("R2_stratum",1:4)
## -------------------------------------------------------------
## Abschnitt 9.2.7, Berücksichtigung der Schätzfehler
## -------------------------------------------------------------
# -------------------------------------------------------------
# Abschnitt 9.2.7, Listing 17: Bestimmung des Erwartungsbereichs
#
vv <- "expTWLE.OLS1" # Variablenname
mm <- "OLS1" # Kurzname des Modells
dfr <- NULL # Ergebnistabelle
# Schleife über alle möglichen Breite von 10 bis 60
for(w in 10:60){
# Variablen für Ergebnisse pro w
var <- paste0(mm,".pos.eb",w) # Position der Schule
var.low <- paste0(mm,".eblow",w) # Untere Grenze des EBs
var.upp <- paste0(mm,".ebupp",w) # Obere Grenze des EBs
# Berechnen
dat[,var.low] <- dat[,vv]-w/2 # Untere Grenze des EBs
dat[,var.upp] <- dat[,vv]+w/2 # Obere Grenze des EBs
# Position: -1=unterhalb, 0=innerhalb, 1=oberhalb des EBs
dat[,var] <- -1*(dat$TWLE < dat[,var.low]) + 1*(dat$TWLE > dat[,var.upp])
# Verteilung der Schulpositionen
tmp <- data.frame(t(matrix(prop.table(table(dat[,var])))))
names(tmp) <- c("unterhalb","innerhalb","oberhalb")
tmp <- data.frame("ModellxBereich"=var,tmp); dfr <- rbind(dfr,tmp) }
# Abweichung zur Wunschverteilung 25-50-25
dfr1 <- dfr
dfr1[,c(2,4)] <- (dfr1[,c(2,4)] - .25)^2
dfr1[,3] <- (dfr1[,3] - .5)^2
dfr1$sumquare <- rowSums(dfr1[,-1])
# Auswahl markieren
dfr$Auswahl <- 1*(dfr1$sumquare == min(dfr1$sumquare) )
# -------------------------------------------------------------
# Abschnitt 9.2.7, Ergänzung zum Buch
# Bestimmung des Erwartungsbereichs
#
# Ergebnisse aller Schulen werden aus Ursprungsdatensatz geladen.
dat <- datenKapitel09
# Liste der Erwartungswerte-Variablen
exp.vars <- grep("expTWLE",names(dat),value=T)
# Modellnamen
m.vars <- gsub("expTWLE.","",exp.vars, fixed = TRUE)
# Liste der Ergebnistabelle
list0 <- list()
# Ergebnisse
tab.erg <- NULL
# Schleife über alle Erwartungswerte aller Modelle
for(ii in 1:length(exp.vars)){
# ii <- 1
vv <- exp.vars[ii]
mm <- m.vars[ii]
# Ergebnistabelle
dfr <- NULL
# Schleife über alle möglichen Breite von 10 bis 60
for(w in 10:60){
# eb <- 10
var <- paste0(mm,".pos.eb",w) # Position der Schule
var.low <- paste0(mm,".eblow",w) # Untere Grenze des EBs
var.upp <- paste0(mm,".ebupp",w) # Obere Grenze des EBs
# Untere Grenze des EBs = Erwartungswert - w/2
dat[,var.low] <- dat[,vv]-w/2
# Obere Grenze des EBs = Erwartungswert + w/2
dat[,var.upp] <- dat[,vv]+w/2
# Position der Schule bestimmen
# -1 = unterhalb, 0 = innterhalb, 1 = oberhalb des EBs
dat[,var] <- -1*(dat$TWLE < dat[,var.low]) + 1*(dat$TWLE > dat[,var.upp])
# Verteilung der Positionen
tmp <- data.frame(t(matrix(prop.table(table(dat[,var])))))
names(tmp) <- c("unterhalb","innerhalb","oberhalb")
tmp <- data.frame("ModellxBereich"=var,tmp)
dfr <- rbind(dfr,tmp)
}
# Vergleich mit Wunschverteilung 25-50-25
dfr1 <- dfr
dfr1[,c(2,4)] <- (dfr1[,c(2,4)] - .25)^2
dfr1[,3] <- (dfr1[,3] - .5)^2
dfr1$sumquare <- rowSums(dfr1[,-1])
# Auswahl markieren
dfr$Auswahl <- 1*(dfr1$sumquare == min(dfr1$sumquare) )
# Zum Liste hinzufügen
list0[[ii]] <- dfr
print(dfr[which(dfr$Auswahl == 1),])
tab.erg <- rbind(tab.erg, dfr[which(dfr$Auswahl == 1),])
}
# Nur gewählte Ergebnisse im Datensatz beibehalten
all.vars <- grep("eb",names(dat),value=T)
# Untere und Obere Grenze mit speichern
eb.vars <- tab.erg[,1]
low.vars <- gsub("pos.eb","eblow",eb.vars)
upp.vars <- gsub("pos.eb","ebupp",eb.vars)
del.vars <- setdiff(all.vars, c(eb.vars,low.vars,upp.vars))
dat <- dat[,-match(del.vars,names(dat))]
## -------------------------------------------------------------
## Appendix: Abbildungen
## -------------------------------------------------------------
# -------------------------------------------------------------
# Abbildung 9.4
#
# Koeffizienten bei der ersten 50 lambdas ausdrucken
# Stratum 4
lambda <- lasso.mod2$lambda[1:50]
a <- round2(lambda,2)
a1 <- a[order(a)]
L <- length(a)
dfr <- NULL
for(ll in 1:L){
dfr.ll <- as.matrix(coef(lasso.mod2,newx = Z,s=a[ll] ))
colnames(dfr.ll) <- paste0("a_",ll)
dfr.ll <- data.frame("coef"=rownames(dfr.ll),dfr.ll)
rownames(dfr.ll) <- NULL
if(ll == 1) dfr <- dfr.ll else dfr <- merge(dfr, dfr.ll)
}
# Ohne Intercept
dfr <- dfr[-1,]
rownames(dfr) <- 1:nrow(dfr)
cl <- colors()
cl <- grep("grey",cl,value=T)
# Umgekehrte Reihenfolge
dfr1 <- dfr
for(x in 2:(L+1)) {dfr1[,x] <- dfr[,(L+3)-x];
names(dfr1)[x] <- names(dfr)[(L+3)-x]}
###
plot(x = log(a), y = rep(0,L), xlim = rev(range(log(a))), ylim=c(-20,22),
type = "l", xaxt ="n", xlab = expression(paste(lambda)),
ylab="Geschätzte Regressionskoeffizienten")
axis(1, at=log(a), labels=a,cex=1)
tmp <- nrow(dfr)
for(ll in 1:tmp){
# ll <- 1
lines(x=log(a),y=dfr[ll,2:(L+1)],type="l",pch=15-ll,col=cl[15-ll])
points(x=log(a),y=dfr[ll,2:(L+1)],type="p",pch=15-ll)
legend(x=2.8-0.7*(ll>tmp/2),y=25-2*(ifelse(ll>7,ll-7,ll)),
legend =dfr$coef[ll],pch=15-ll,bty="n",cex=0.9)
}
# Kennzeichung der gewählten lambda
v <- log(lasso.mod2$lambda.min)
lab2 <- expression(paste("ausgewähltes ",lambda," = .43"))
text(x=v+0.6,y=-8,labels=lab2)
abline(v = v,lty=2,cex=1.2)
# -------------------------------------------------------------
# Abbildung 9.5
# Auswahl Lambda anhand min(cvm)
#
xlab <- expression(paste(lambda))
plot(lasso.mod2, xlim = rev(range(log(lambda))),
ylim=c(550,1300),xlab=xlab,xaxt ="n",
ylab = "Mittleres Fehlerquadrat der Kreuzvalidierung (cvm)",
font.main=1,cex.main=1)
axis(1, at=log(a), labels=a,cex=1)
lab1 <- expression(paste(lambda," bei min(cvm)"))
text(x=log(lasso.mod2$lambda.min)+0.5,y=max(lasso.mod2$cvm)-50,
labels=lab1,cex=1)
lab2 <- expression(paste("(ausgewähltes ",lambda," = .43)"))
text(x=log(lasso.mod2$lambda.min)+0.6,y=max(lasso.mod2$cvm)-100,
labels=lab2,cex=1)
abline(v = log(lasso.mod2$lambda.min),lty=2)
text(x=log(lasso.mod2$lambda.min)-0.3,y = min(lasso.mod2$cvm)-30,
labels="min(cvm)",cex=1 )
abline(h = min(lasso.mod2$cvm),lty=2)
text <- expression(paste("Anzahl der Nicht-null-Koeffizienten (",
lambda," entsprechend)"))
mtext(text=text,side=3,line=3)
# -------------------------------------------------------------
# Abbildung 9.6
# Rohwert-Schätzwert Schule 1196 & 1217 im Vergleich
#
id <- c(1196, 1217)
par(mai=c(1.2,3,1,.5))
plot(x=rep(NA,2),y=c(1:2),xlim=c(470,610),yaxt ="n",type="l",
xlab="Erwartungswerte je nach Modell und Schulleistung",ylab="")
legend <- c("Schulleistung (TWLE)",paste0("", c("OLS1","OLS2","Lasso1",
"Lasso2","NP"),
"-Modell"))
axis(2, at=c(seq(1,1.4,0.08),seq(1.6,2,0.08)), las=1,cex=0.7,
labels=rep(legend,2))
text <- paste0("Schule ",id)
mtext(text=text,side=2,at = c(1.2,1.8),line = 10)
exp.vars <- c("TWLE",
paste0("expTWLE.", c("OLS1","OLS2","Lasso1","Lasso2","np")))
pch = c(19, 0,3,2,4,5)
ii <- 1
col = c("grey", rep("lightgrey",5))
for(vv in exp.vars){
# vv <- "TWLE"
x <- dat0[which(dat0$idschool %in% id),vv]
abline(h = c(0.92+ii*0.08,1.52+ii*0.08), lty=1+1*(ii>1),col=col[ii])
points(x=x,y=c(0.92+ii*0.08,1.52+ii*0.08),type="p",pch=pch[ii])
ii <- ii + 1
}
## End(Not run)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.