Die Rekonstruktion der relativen Einkommensentwicklung sozialer Gruppen auf der Grundlage von gruppierten Einkommensangaben - Replikations-Code
Inhaltsverzeichnis
Datenaufbereitung
(Der im folgenden beschriebene Code ist auch in der Datei setup-allbus-cumul.R enthalten.)
Zunächst laden wir die vom GESIS-Datenarchiv herunter geladene SAV-Datei in
unsere R-Sitzung. Dazu aktivieren wir zunächst das Zusatz-Paket "memisc", das
unter https://cran.r-project.org/package=memisc herunter geladen werden kann.
library(memisc)
Anschließend entpacken wir die SAV-Datei aus dem ZIP-Archiv und machen die
entpackte mit R bekannt.
ZA4582.sav <- unzip("ZA4582_v1-0-0.sav.zip", exdir="/tmp") ZA4582.sav <- spss.system.file(ZA4582.sav)
Wir lassen uns dann umgehend eine Zusammenfassung der Daten in der SAV-Datei ausgeben.
ZA4582.sav
SPSS system file '/tmp/ZA4582_v1-0-0.sav' with 1936 variables and 61194 observations
Da wir nicht alle Variablen in dem Datensatz benötigen, bilden wir ein "Subset"
der Daten mit Hilfe der gleichnamigen Funktion subset(). Da die Namen der
Variablen in dem gelieferten Datensatz alles andere als mnemonisch sind, lassen
wir subset() die Variablen sogleich umbenennen.
allbus.work <- subset(ZA4582.sav, select = c( year = year, respid = respid, ost_west = ost_west, german = german, samptype = samptype, partypref = v24, voteint = v25, voteint_berlin = v26, birthyear = v727, age = v729, age_categ = v730, sex = v731, terwey_occmajor = v863, terwey_goldthorpe = v874, pers_income_open = v913, pers_income_list = v914, pers_income = v915, pers_income_categ = v916, hh_income_open = v918, hh_income_list = v919, hh_income = v920, hh_income_categ = v921, oecd_income = v924, oecd_income_categ = v925, educ.certif = v746, wghtpt = wghtpt, wghtpow = wghtpow, wghtptow = wghtptow, wghtht = wghtht, wghthow = wghthow, wghthtow = wghthtow ) )
Einige der Variablen müssen zunächst als metrisch (intervall- oder ratio-skaliert) deklariert werden, da einige ihrer Werte ettiketiert sind und das normalerweise dazu führt, dass die Variablen als kategorial (nominal oder ordinal) interpretiert werden.
allbus.work <- within(allbus.work,{ foreach(v = c(birthyear,age, pers_income_open, hh_income_open, pers_income, hh_income, oecd_income, wghtpt, wghtpow, wghtptow, wghtht, wghthow, wghthtow ), measurement(v) <- "interval") })
Nachdem wir unseren Arbeitsdatensatz erstellt haben, können wir ihn für die
weiteren Analysen aufbereiten. Zunächst rekodieren wir die Variablen
pers_income_categ, hh_income_categ, und oecd_income_categ in jeweils zwei
neue numerische Variablen so um, dass jeweils eine den Wert der unteren Grenze
des durch den jeweiligen Code repräsentierten Einkommensintervalls enthält und
eine andere den Wert der oberen Grenze.
allbus.work <- within(allbus.work,{ foreach(v = c(pers_income_categ, hh_income_categ, oecd_income_categ), u = c(pers_income_upr, hh_income_upr, oecd_income_upr), l = c(pers_income_lwr, hh_income_lwr, oecd_income_lwr),{ u <- recode(v, 0 -> 0, 1 -> 200, 2 -> 300, 3 -> 400, 4 -> 500, 5 -> 625, 6 -> 750, 7 -> 875, 8 -> 1000, 9 -> 1125, 10 -> 1250, 11 -> 1375, 12 -> 1500, 13 -> 1750, 14 -> 2000, 15 -> 2250, 16 -> 2500, 17 -> 2750, 18 -> 3000, 19 -> 4000, 20 -> 5000, 21 -> 7500, 22 -> Inf) missing.values(u) <- NULL labels(u) <- NULL description(u) <- "" l <- recode(v, 0:1 -> 0, 2 -> 200, 3 -> 300, 4 -> 400, 5 -> 500, 6 -> 625, 7 -> 750, 8 -> 875, 9 -> 1000, 10 -> 1125, 11 -> 1250, 12 -> 1375, 13 -> 1500, 14 -> 1750, 15 -> 2000, 16 -> 2250, 17 -> 2500, 18 -> 2750, 19 -> 3000, 20 -> 4000, 21 -> 5000, 22 -> 7500) missing.values(l) <- NULL labels(l) <- NULL description(l) <- "" }) })
Mit description() erhalten wir eine Auflistung der Variablen mit ihren
Ettiketten. Wir schreiben das Ergebnis des Aufrufs dieser Funktion umgehend in
eine Datei namens allbus-work-desc.txt.
Write(description(allbus.work), file="allbus-work-desc.txt")
Ein Codebuch erhalten wir mit codebook(). Das Ergebnis des entsprechenden
Funktionsaufrufes schreiben wir auch in eine Datei names allbus-work-cdbk.txt
Write(codebook(allbus.work), file="allbus-work-cdbk.txt")
Jetzt ist der Datensatz soweit vorbereitet, dass wir ihn für die weitere Verwendung abspeichern können.
save(allbus.work, file="allbus-work.RData")
Untersuchung der Verteilung des kategorisierten Einkommens
Im folgenden wird der Code vorgestellt und demonstiert, der die Graphiken in Abbildung 1, sowei einige ergänzende Graphiken erstellt. Er ist in der Datei income-dist.R enthalten.
Vorbereitungen
Wenn wir mit den folgenden Code in einer neuen R-Sitzung ausführen wollen, müssen wir das benötigte R-Paket "memisc" erneut aktivieren, und den Arbeitsdatensatz neu laden.
library(memisc) load("allbus-work.RData")
Aus der R-Script-Datei plot-stepfun.R holen wir uns nun eine Funktion, die uns erlaubt, elegant und schnell Stufenfunktionen zu erstellen.
source("plot-stepfun.R")
Wir definieren noch eine Hilfsfunktion, die aus einem Datenvektor alle endlichen Werte herausfiltert. Dies ist
keep.finite <- function(x){ keep <- is.finite(x) x[keep] }
Die im Folgenden definierte Funktion income_plot ist das "Arbeitspferd" für
die Herstellung der folgenden Graphiken. Sie erzeugt einen Diagramm der
empirischen kumulativen Verteilungsfunktion des Einkommens, entweder in auf der
Originalskala der Daten (für log=FALSE) oder auf einer logarithmierten Skala.
Die empirische kumulative Verteilungsfunktion wird dann jeweils mit der
kumulativen Verteilungsfunktion der standardisierten (Log-)Normalverteilung
verglichen.
income_plot <- function(x,log=FALSE,lty=3,...){ x <- keep.finite(x) if(log){ x <- keep.finite(log(x)) m <- mean(x) s <- sd(x) plot.ecdf(ecdf(x),xaxt="n",...,xlim=c(log(70),log(70000))) curve(pnorm(x,mean=m,sd=s),add=TRUE, lty=lty, from=0.001) plot.ecdf(ecdf(x),...,add=TRUE) axis(1,at=log(c(50,100,200,500,1000,2000,5000,10000,20000,50000)), labels=c(50,100,200,500,1000,2000,5000,10000,20000,50000) ) } else{ plot.ecdf(ecdf(x),...) log.x <- keep.finite(log(x)) m <- mean(log.x) s <- sd(log.x) curve(pnorm(log(x),mean=m,sd=s),add=TRUE, lty=lty, from=0.001) plot.ecdf(ecdf(x),...,add=TRUE) } }
Verteilung des Einkommens in 1986
Hier erzeugen wir uns einen Teildatensatz, der Daten nur für das Jahr 1986 enthält.
allbus.work.1986 <- subset(allbus.work, year == 1986) allbus.work.1986 <- as.data.frame(allbus.work.1986)
Der folgende Code erzeugt ein Diagramm der Verteilung der oberen Grenze der Kategorien des Haushalts-Nettoeinkommens in 1986.
with(allbus.work.1986, income_plot(hh_income_upr, ylab="Kumulative Verteilung",pch=21,lty=2, xlab="Haushaltseinkommen",bg="gray80", main=NULL))
Der folgende Code erzeugt ein Diagramm der Verteilung des offen abgefragten Haushalts-Nettoeinkommens in 1986.
with(allbus.work.1986, income_plot(hh_income_open, ylab="Kumulative Verteilung",pch=21,lty=2, xlab="Haushaltseinkommen",bg="gray80", main=NULL))
Der folgende Code erzeugt ein Diagramm der Verteilung der oberen Grenze der Kategorien des Haushalts-Nettoeinkommens in 1986 (auf einer logarithmischen Skala).
with(allbus.work.1986, income_plot(hh_income_upr, log=TRUE, ylab="Kumulative Verteilung",pch=21,lty=2, xlab="Haushaltseinkommen (Log-Skala)",bg="gray80", main=NULL))
Der folgende Code erzeugt ein Diagramm der Verteilung des offen abgefragten Haushalts-Nettoeinkommens in 1986 (auf einer logarithmischen Skala).
with(allbus.work.1986, income_plot(hh_income_open, log=TRUE, ylab="Kumulative Verteilung",pch=21,lty=2, xlab="Haushaltseinkommen (Log-Skala)",bg="gray80", main=NULL))
Verteilung des Einkommens in 2014
Hier erzeugen wir uns einen Teildatensatz, der Daten nur für das Jahr 2014 enthält.
allbus.work.2014 <- subset(allbus.work, year == 2014) allbus.work.2014 <- as.data.frame(allbus.work.2014)
Der folgende Code erzeugt ein Diagramm der Verteilung der oberen Grenze der Kategorien des Haushalts-Nettoeinkommens in 2014.
with(allbus.work.2014, income_plot(hh_income_upr, ylab="Kumulative Verteilung",pch=21,lty=2, xlab="Haushaltseinkommen",bg="gray80", main=NULL))
Der folgende Code erzeugt ein Diagramm der Verteilung des offen abgefragten Haushalts-Nettoeinkommens in 2014.
with(allbus.work.2014, income_plot(hh_income_open, ylab="Kumulative Verteilung",pch=21,lty=2, xlab="Haushaltseinkommen",bg="gray80", main=NULL))
Der folgende Code erzeugt ein Diagramm der Verteilung der oberen Grenze der Kategorien des Haushalts-Nettoeinkommens in 2014 (auf einer logarithmischen Skala).
with(allbus.work.2014, income_plot(hh_income_upr, log=TRUE, ylab="Kumulative Verteilung",pch=21,lty=2, xlab="Haushaltseinkommen (Log-Skala)",bg="gray80", main=NULL))
Der folgende Code erzeugt ein Diagramm der Verteilung des offen abgefragten Haushalts-Nettoeinkommens in 2014 (auf einer logarithmischen Skala).
with(allbus.work.2014, income_plot(hh_income_open, log=TRUE, ylab="Kumulative Verteilung",pch=21,lty=2, xlab="Haushaltseinkommen (Log-Skala)",bg="gray80", main=NULL))
Die Simulation von Zufallszahlen aus einer trunkierten Normalverteilung
Im folgenden wird der Code vorgestellt und demonstiert, der die Graphik in Abbildung 2 erstellt. Er ist in der Datei demo-normscores.R enthalten.
Die Datei normscores.R enthält Code unter anderem zur Berechnung des
begingten Erwartungswertes einer normalverteilten Zufallvariablen, gegeben dass
ihr Wert in einem bestimmten Intervall liegt. Allerdings enthält die Datei auch
Code zur Erzeugung von Zufallszahlen, die einer auf ein Intervall restringierten
Normalverteilung folgen. Diese Zufallszahlen, die sich zum Zweck der Imputation
eignen, werden im Folgenden demonstriert.
source("normscores.R")
Mit der in der Datei normscores.R definierten Funktion
rnorm_intvl()
erzeugen wir jetzt 50.000 Zufallszahlen, die im Intervall [-1;+2] liegen:
y.star <- rnorm_intvl(50000,-1,2)
Alsdann erzeugen wir ein Histogramm dieser Zufallszahlen.
hist(y.star,breaks=25, col="gray70", xlab=expression(y^"*"), ylab="Häufigkeit",main=NULL)
Anwendung: Relative Verteilung des Einkommens von Angehörigen unterschiedlicher Erwerbsklassen
(Der im folgenden beschriebene Code ist auch in der Datei income-gclass.R enthalten.)
Vorbereitungen
Um den folgenden Code unabhängig vom vorangegangenen ausführen zu können,
aktivieren wir das R-Paket "memisc". Weiterhin laden wir aus der R-Scriptdatei
normscores.R die function normscores(), mit deren Hilfe wir Scores für
latente standard-normalverteilte Einkommensvariable ermitteln. Um die
Standardisierung separat jeweils für jedes Erhebungsjahr und für Ost und West
vornehmen zu können, benötigen wir ferner die Funktion groupwise(), die wir
aus der R-Scriptdatei groupwise.R laden. Abschließend laden wir den aus dem
kumulierten ALLBUS extrahierten Arbeitsdatensatz.
library(memisc) source("normscores.R") source("groupwise.R") load("allbus-work.RData")
Nachdem wir die benötigten Funkionen und Daten geladen haben, können wir jetzt
zur Aufbereitung der Daten schreiten. Das Objekt allbus.work ist ein
"data.set"-Objekt, d.h. es gehört einer im Paket "memisc" definierten Klasse
an. Objekte der Klasse "data.set" erlauben es, den darin enthaltenen Variablen
Werte-Ettiketten zuzuordnen oder für sie ein bestimmtes Skalenniveau
zuzuordnen. Diese den Variablen zugeordneten Attribute entsprechend dem, was man
bereits aus statistischen Software-Paketen wie SPSS oder Stata kennt, aber
nicht bzw. nicht ohne das "memisc"-Paket von R unterstützt wird. Diese
Attribute bestimmen unter anderem, wie diese Variablen behandelt werden, wenn
ein "data.set"-Objekt in ein für die Datenanalyse mit R vorgesehenen "data
frame" (d.h. ein Objekt aus der Klasse "data.frame") umgewandelt wird.
Konkret entfernt der folgende Code die Werte-Ettiketten von der Variable year
und ordnet ihr ein intervallskaliertes Messniveau zu. Der Variablen Variablen
ost_west werden passende Werte-Ettiketten zugeordnet. Es wird eine Faktorvariable
year.ow konstruiert, die später verwendet wird, um für beide Teile
Deutschlands und jedes Jahr standardisiserte Einkommens-Scores zu erhalten.
Eine weitere Vorbereitung geschieht in der Umkodierung der Variablen
terwey_goldthorpe, deren Codes auf 10 Zusammengefasst werden und der Variablen
gclass zugeordnet werden. In dem dazu verwendeten Aufruf der Funktion
recode() werden den neuen Codes auch gleich Werte-Ettiketten zugeordnet.
Schließlich wird eine kategorisierte Altersvariable mit dem Namen age_categ
konstruiert, die es erlaubt, die nach Alter unterschiedliche Entwicklung des
Einkommens nachzuzeichnen.
Schließlich werden vier Variablen hh_income_sc, hh_income_lsc,
pers_income_sc, und pers_income_sc
konstruiert, die aus dem kategorisierten
Haushalts- bzw. persönlichen Einkommen standardisierte Scores für die dem
Einkommen zu Grunde liegende latente Variable enthalten, und zwar jeweils in der
normalverteilten und der log-normalverteilten Variante.
allbus.work <- within(allbus.work,{ labels(year) <- NULL measurement(year) <- "interval" labels(ost_west) <- c("Alte Bundesländer"=1, "Neue Bundesländer"=2) year.ow <- interaction(as.factor(ost_west),as.factor(year)) gclass <- recode(terwey_goldthorpe, "Un/angelernte Arbeiter" = 1 <- 9, "Facharbeiter" = 2 <- 8, "Vorarbeiter, Techniker" = 3 <- 7, "Nichtmanuelle Routine" = 4 <- 11, "Einfache Bürotätigkeit" = 5 <- 3, "Untere Dienstklasse" = 6 <- 2, "Obere Diensklasse" = 7 <- 1, "Selbständige" = 8 <- 5, "Unternehmenseigner" = 9 <- 4, "Landwirtschaftliche Berufe" = 10 <- c(6,10,12) ) age_categ <- cases("18-25" = (18 <= age & age <= 25), "26-35" = (26 <= age & age <= 35), "36-45" = (36 <= age & age <= 45), "46-55" = (46 <= age & age <= 55), "56-65" = (56 <= age & age <= 65), "66+" = (66 <= age) ) hh_income_sc <- groupwise(normscores,hh_income_categ,by=year.ow) hh_income_lsc <- exp(hh_income_sc) pers_income_sc <- groupwise(normscores,pers_income_categ,by=year.ow) pers_income_lsc <- exp(pers_income_sc) })
Für die spätere Analyse der Einkommensentwicklung in der Klasse der un- und angelernten Arbeiter bilden wir einen entsprechenden Teildatensatz.
allbus.semi.unskilled <- subset(allbus.work, gclass==1)
Nach dem wir die Daten aufbereitet haben, wandeln wir sie in
"data.frame"-Objekte um.
allbus.work <- as.data.frame(allbus.work) allbus.semi.unskilled <- as.data.frame(allbus.semi.unskilled)
Haushaltseinkommen nach Erwerbsklasse
Zur Analyse der Entwicklung des Einkommens nach Erwerbsklasse fassen wir die verschiedenen Einkommensvariablen nach Jahr, Landesteil und Erwerbsklasse zusammen: Für jede Kombination von Werten dieser drei Variablen bilden wir den Mittelwert der Einkommensvariablen.
income.gclass <- Aggregate( c( hh_income=mean(hh_income,na.rm=TRUE), hh_income_lsc=mean(hh_income_lsc,na.rm=TRUE), hh_income_sc=mean(hh_income_sc,na.rm=TRUE), )~year+ost_west+gclass, data=allbus.work )
Bevor wir zur graphischen Darstellung der Einkommensentwicklung schreiten, aktivieren wir das dazu nötige R-Zusatzpaket "lattice".
library(lattice)
Der folgende Aufruf der Funktion xyplot() erzeugt Liniendiagramme mit
Kurvenglättung (type="smooth") des Haushaltseinkommens in offener Abfrage
(hh_income) nach Erhebungsjahr (year), aufgeschlüsselt nach Erwerbsklasse
(groups=gclass).
(hh_income.xyplot <- xyplot(hh_income~year, data=income.gclass, subset=(ost_west=="Alte Bundesländer"), type="smooth",span=1/2, groups=gclass, lwd=2, ylab="Haushaltseinkommen", xlab="Jahr", auto.key=list(points=FALSE,lines=TRUE,space="top", columns=2), par.settings=list( par.main.text = list( font=3,cex=1 ), superpose.line = list( lty=1:5, lwd=2 ) ) ))
Der nächste Aufruf der Funktion xyplot() zeigt ein dem vorangegangenen
analoges Diagramm mit einer anderen abhängigen Variable, den standardisierten
normalverteilten Einkommensscores. Die allgemeine Einkommensentwicklung ist hier
gewissermaßen herausgerechnet, dadurch, dass das die Score-Mittelwerte separat
für jedes Jahr ermittelt werden.
(hh_income_sc.xyplot <- xyplot(hh_income_sc~year, data=income.gclass, subset=(ost_west=="Alte Bundesländer"), type="smooth",span=1/2, #type="l", groups=gclass, lwd=2, ylab="Einkommensscores (Normalverteilung)", xlab="Jahr", auto.key=list(points=FALSE,lines=TRUE,space="top", columns=2), par.settings=list( par.main.text = list( font=3,cex=1 ), superpose.line = list( lty=1:5, lwd=2 ) ) ))
Der dritte Aufruf der Funktion xyplot() stellt die exponentiierten
standardisierten Einkommensscores dar. Deren Verteilung zu jeden Zeitpunkt ist
der Verteilung des "rohen" Einkommens ähnlich. Aber immer noch ist die
Gesamtentwicklung herausgerechnet, so dass die relative Entwicklung des nach
Erwerbsklassen aufgeschlüsselten Einkommens deutlicher hervortritt.
(hh_income_lsc.xyplot <- xyplot(hh_income_lsc~year, data=income.gclass, subset=(ost_west=="Alte Bundesländer"), type="smooth",span=1/2, #type="l", groups=gclass, lwd=2, ylab="Einkommensscores (Log-Normalverteilung)", xlab="Jahr", auto.key=list(points=FALSE,lines=TRUE,space="top", columns=2), par.settings=list( par.main.text = list( font=3,cex=1 ), superpose.line = list( lty=1:5, lwd=2 ) ) ))
Persönliches Einkommen von Angehörigen der Klasse der un- und angelernten Arbeiter nach Alter
In gleicher Weise wie wir zuvor Scores für das Haushaltseinkommen konstruiert haben, die für jedes Jahr, sowie für Ost- und West-Deutschland standardisiert haben, konstruieren wir jetzt Scores für das persönliche Einkommen, für jedes Jahr und jeden Landesteil, die zusätzlich noch mit Bezug auf die Klasse der un- und angelernten Arbeiter standardisiert sind. Die Scores erlauben es, die relative Einkommensentwicklung der verschiedenen Altersgruppen im Vergleich innerhalb dieser Klasse nachzuzeichnen.
income.semi.unskilled.age <- Aggregate( c( p_income=mean(pers_income,na.rm=TRUE), p_income_lsc=mean(pers_income_lsc,na.rm=TRUE), p_income_sc=mean(pers_income_sc,na.rm=TRUE) )~year+ost_west+age_categ, data=allbus.semi.unskilled )
Der folgende Aufruf der Funktion xyplot() erzeugt Liniendiagramme mit
Kurvenglättung (type="smooth") des persönlichen Einkommens in offener Abfrage
(pers_income) nach Erhebungsjahr (year), aufgeschlüsselt nach Altersgruppe
(groups=age_categ).
(p_income.susk.age.xyplot <- xyplot(p_income~year, data=income.semi.unskilled.age, subset=(ost_west=="Alte Bundesländer"), type="smooth",span=1/2, groups=age_categ, lwd=2, ylab="Persönl. Einkommen", xlab="Jahr", auto.key=list(points=FALSE,lines=TRUE,space="top", columns=2), par.settings=list( par.main.text = list( font=3,cex=1 ), superpose.line = list( lty=1:5, lwd=2 ) ) ))
Der folgende Aufruf von xyplot() erzeugt ein Diagramm der Entwicklung des
persönlichen Einkommens innerhalb der Klasse der un-/angelernten Arbeiter
gemessen durch die Mittelwerte der standardisierten normalverteilten Scores.
(p_income_sc.susk.age.xyplot <- xyplot(p_income_sc~year, data=income.semi.unskilled.age, subset=(ost_west=="Alte Bundesländer"), type="smooth",span=1/2, groups=age_categ, lwd=2, ylab="Einkommensscores (Normalverteilung)", xlab="Jahr", auto.key=list(points=FALSE,lines=TRUE,space="top", columns=2), par.settings=list( par.main.text = list( font=3,cex=1 ), superpose.line = list( lty=1:5, lwd=2 ) ) ))
In der letzten Graphik wird die Entwicklung des persönlichen Einkommens innerhalb der Erwerbsklasse der un- und angelernten Arbeiter wieder auf der log-normalverteilten Skala dargestellt die der eigentlichen Verteilung des Einkommens entspricht – abgesehen davon, dass das Einkommen mit Bezug auf Jahr, Landesteil und Erwerbsklasse standardisiert ist.
(p_income_lsc.susk.age.xyplot <- xyplot(p_income_lsc~year, data=income.semi.unskilled.age, subset=(ost_west=="Alte Bundesländer"), type="smooth",span=1/2, groups=age_categ, lwd=2, ylab="Einkommensscores (Log-Normalverteilung)", xlab="Jahr", auto.key=list(points=FALSE,lines=TRUE,space="top", columns=2), par.settings=list( par.main.text = list( font=3,cex=1 ), superpose.line = list( lty=1:5, lwd=2 ) ) ))
