## R Code to generate graphics for ''To Work or Not to Work?''
## Copyright (c) 2008 Aleksandr Andreev.
##
## This code is FREE SOFTWARE, you may use, modify and / or distribute it
## provided that you give credit to the original author.
## While this code is provided in the hope that it may be useful,
## it comes with ABSOLUTELY NO WARRANTY.
##
library(TeachingDemos)
library(HH)

data 	<- read.csv("postestimation.csv")
data$a001ter    <- as.factor(data$a001ter)
data$male       <- as.factor(data$male)
data$residency  <- as.factor(data$residency)
data$status     <- as.factor(data$status)

attach(data)
levels(male)      <- c("Women", "Men")
levels(residency) <- c("Village", "PGT", "City", "Large City")

# Generate a plot of probability by age and sex
par(mfrow=c(2, 2))
plot(supsmu(age, probworkinghealthymale, bass=8), xlab='Age', ylab='Probability', main='Healthy Population', ylim=c(0,1), type='l')
lines(supsmu(age, probworkinghealthyfemale, bass=8), lty=2)
plot(supsmu(age, probworkinggroup1male, bass=10), xlab='Age', ylab='Probability', main='Group I Disabled Population', ylim=c(0,1), type='l')
lines(supsmu(age, probworkinggroup1female, bass=10), lty=2)
plot(supsmu(age, probworkinggroup2male, bass=10), xlab='Age', ylab='Probability', main='Group II Disabled Population', ylim=c(0,1), type='l')
lines(supsmu(age, probworkinggroup2female, bass=10), lty=2)
plot(supsmu(age, probworkinggroup3male, bass=10), xlab='Age', ylab='Probability', main='Group III Disabled Population', ylim=c(0,1), type='l')
lines(supsmu(age, probworkinggroup3female, bass=10), lty=2)

par(xpd=NA)
tmp <- cnvrt.coords(0.47,0.08, 'tdev')$usr 
legend(tmp, legend=c("Men", "Women"), lty=c(1, 2), bty='n')

# generate the residency plot
par(mfrow=c(2, 2))
at <- c(1.5, 2, 3.25, 3.75, 5, 5.5, 6.75, 7.25)
boxplot(probworkinghealthy ~ interaction(male, residency), boxwex=0.2,  at = at, axes=F, ylab="Probability", main="Healthy Population", col=c("grey", "white"))
axis(2)
axis(1, labels=rep(c("Women", "Men"), 4), at=at, cex.axis=.7, las=3)
 mtext(side = 1, line=3.2, c("Village", "PGT", "City", "Large City"), at=c(1.75, 3.5, 5.25, 7), cex=.8)
box()

boxplot(probworkinggroup1 ~ interaction(male, residency), boxwex=0.2,  at = at, axes=F, ylab="Probability", main="Group I Disabled Population", col=c("grey", "white"))
axis(2)
axis(1, labels=rep(c("Women", "Men"), 4), at=at, cex.axis=.7, las=3)
 mtext(side = 1, line=3.2, c("Village", "PGT", "City", "Large City"), at=c(1.75, 3.5, 5.25, 7), cex=.8)
box()

boxplot(probworkinggroup2 ~ interaction(male, residency), boxwex=0.2,  at = at, axes=F, ylab="Probability", main="Group II Disabled Population", col=c("grey", "white"), ylim=c(0,1))
axis(2)
axis(1, labels=rep(c("Women", "Men"), 4), at=at, cex.axis=.7, las=3)
 mtext(side = 1, line=3.2, c("Village", "PGT", "City", "Large City"), at=c(1.75, 3.5, 5.25, 7), cex=.8)
box()

boxplot(probworkinggroup3 ~ interaction(male, residency), boxwex=0.2,  at = at, axes=F, ylab="Probability", main="Group III Disabled Population", col=c("grey", "white"))
axis(2)
axis(1, labels=rep(c("Women", "Men"), 4), at=at, cex.axis=.7, las=3)
 mtext(side = 1, line=3.2, c("Village", "PGT", "City", "Large City"), at=c(1.75, 3.5, 5.25, 7), cex=.8)
box()

# generate the family status plots
par(mfrow=c(2, 2))
at <- c(1.5, 2, 3.25, 3.75, 5, 5.5, 6.75, 7.25)
boxplot(probworkinghealthy ~ interaction(male, status), boxwex=0.2,  at = at, axes=F, ylab="Probability", main="Healthy Population", col=c("grey", "white"))
axis(2)
axis(1, labels=rep(c("Women", "Men"), 4), at=at, cex.axis=.7, las=3)
 mtext(side = 1, line=3.2, c("Single", "Married", "Divorced", "Widowed"), at=c(1.5, 3.25, 5.25, 7.25), cex=.75)
box()

boxplot(probworkinggroup1 ~ interaction(male, status), boxwex=0.2,  at = at, axes=F, ylab="Probability", main="Group I Disabled Population", col=c("grey", "white"))
axis(2)
axis(1, labels=rep(c("Women", "Men"), 4), at=at, cex.axis=.7, las=3)
 mtext(side = 1, line=3.2, c("Single", "Married", "Divorced", "Widowed"), at=c(1.5, 3.25, 5.25, 7.25), cex=.75)
box()

boxplot(probworkinggroup2 ~ interaction(male, status), boxwex=0.2,  at = at, axes=F, ylab="Probability", main="Group II Disabled Population", col=c("grey", "white"), ylim=c(0,1))
axis(2)
axis(1, labels=rep(c("Women", "Men"), 4), at=at, cex.axis=.7, las=3)
 mtext(side = 1, line=3.2, c("Single", "Married", "Divorced", "Widowed"), at=c(1.5, 3.25, 5.25, 7.25), cex=.75)
box()

boxplot(probworkinggroup3 ~ interaction(male, status), boxwex=0.2,  at = at, axes=F, ylab="Probability", main="Group III Disabled Population", col=c("grey", "white"))
axis(2)
axis(1, labels=rep(c("Women", "Men"), 4), at=at, cex.axis=.7, las=3)
 mtext(side = 1, line=3.2, c("Single", "Married", "Divorced", "Widowed"), at=c(1.5, 3.25, 5.25, 7.25), cex=.75)
box()

# MMC plot comparing Tyva with Moscow
data3 <- read.csv("tempost.csv")
data3$region <- as.factor(data3$region)
attach(data3)
levels(region) <- c("Moscow Men", "Moscow Women", "Tyva Men", "Tyva Women")

# ANOVA model
regions.aov <- aov(probability ~ region)

par (cex = .7)
plot(glht.mmc(regions.aov, linfct = mcp(region = "Tukey")), main="Mean-Mean Comparison for Healthy Individuals")

## Generate plots of changes in probability
## following decrease of disability pensions
data2 <- read.csv("differences.csv")
data2$a001ter   <- as.factor(data2$a001ter)
data2$male      <- as.factor(data2$male)
data2$residency <- as.factor(data2$residency)
data2$status    <- as.factor(data2$status)

attach(data2)

par(mfrow=c(3, 3))
plot(supsmu(age, diffgroup1male, bass=10), xlab='Age', ylab='Change in Probability', main='Group I Disabled Population', ylim=c(-0.02, 0.01), type='l')
lines(supsmu(age, diffgroup1female, bass=10), lty=2)
abline(0, 0, lty=3)
plot(supsmu(age, diffgroup2male, bass=10), xlab='Age', ylab='Change in Probability', main='Group II Disabled Population', ylim=c(-0.02, 0.01), type='l')
lines(supsmu(age, diffgroup2female, bass=10), lty=2)
abline(0, 0, lty=3)
plot(supsmu(age, diffgroup3male, bass=10), xlab='Age', ylab='Change in Probability', main='Group III Disabled Population', ylim=c(-0.02, 0.01), type='l')
lines(supsmu(age, diffgroup3female, bass=10), lty=2)
abline(0, 0, lty=3)
tmp <- cnvrt.coords(0.7,0.3, 'tdev')$usr 
legend(tmp, legend=c("Men", "Women"), lty=c(1, 2), bty='y')

at <- c(1.5, 2, 3.25, 3.75, 5, 5.5, 6.75, 7.25)
boxplot(diffgroup1 ~ interaction(male, residency), boxwex=0.2,  at = at, axes=F, ylab="Change in Probability", main="Group I Disabled Population", col=c("grey", "white"), ylim=c(-0.02, 0.01))
abline(0, 0, lty=3)
axis(2)
axis(1, labels=rep(c("Women", "Men"), 4), at=at, cex.axis=.7, las=3)
 mtext(side = 1, line=3.2, c("Village", "PGT", "City", "Large City"), at=c(1.75, 3.5, 5.25, 7), cex=.8)
box()

boxplot(diffgroup2 ~ interaction(male, residency), boxwex=0.2,  at = at, axes=F, ylab="Change in Probability", main="Group II Disabled Population", col=c("grey", "white"), ylim=c(-0.02, 0.01))
abline(0, 0, lty=3)
axis(2)
axis(1, labels=rep(c("Women", "Men"), 4), at=at, cex.axis=.7, las=3)
 mtext(side = 1, line=3.2, c("Village", "PGT", "City", "Large City"), at=c(1.75, 3.5, 5.25, 7), cex=.8)
box()

boxplot(diffgroup3 ~ interaction(male, residency), boxwex=0.2,  at = at, axes=F, ylab="Change in Probability", main="Group III Disabled Population", col=c("grey", "white"), ylim=c(-0.02, 0.01))
abline(0, 0, lty=3)
axis(2)
axis(1, labels=rep(c("Women", "Men"), 4), at=at, cex.axis=.7, las=3)
 mtext(side = 1, line=3.2, c("Village", "PGT", "City", "Large City"), at=c(1.75, 3.5, 5.25, 7), cex=.8)
box()

at <- c(1.5, 2, 3.25, 3.75, 5, 5.5, 6.75, 7.25)
boxplot(diffgroup1 ~ interaction(male, status), boxwex=0.2,  at = at, axes=F, ylab="Change in Probability", main="Group I Disabled Population", col=c("grey", "white"), ylim=c(-0.02, 0.01))
abline(0, 0, lty=3)
axis(2)
axis(1, labels=rep(c("Women", "Men"), 4), at=at, cex.axis=.7, las=3)
 mtext(side = 1, line=3.2, c("Single", "Married", "Divorced", "Widowed"), at=c(1.5, 3.25, 5.25, 7.25), cex=.75)
box()

boxplot(diffgroup2 ~ interaction(male, status), boxwex=0.2,  at = at, axes=F, ylab="Change in Probability", main="Group II Disabled Population", col=c("grey", "white"), ylim=c(-0.02, 0.01))
abline(0, 0, lty=3)
axis(2)
axis(1, labels=rep(c("Women", "Men"), 4), at=at, cex.axis=.7, las=3)
 mtext(side = 1, line=3.2, c("Single", "Married", "Divorced", "Widowed"), at=c(1.5, 3.25, 5.25, 7.25), cex=.75)
box()

boxplot(diffgroup3 ~ interaction(male, status), boxwex=0.2,  at = at, axes=F, ylab="Change in Probability", main="Group III Disabled Population", col=c("grey", "white"), ylim=c(-0.02, 0.01))
abline(0, 0, lty=3)
axis(2)
axis(1, labels=rep(c("Women", "Men"), 4), at=at, cex.axis=.7, las=3)
 mtext(side = 1, line=3.2, c("Single", "Married", "Divorced", "Widowed"), at=c(1.5, 3.25, 5.25, 7.25), cex=.75)
box()

## Generate Map of the Russian Federation
## Hat-tip: Thanks to Roger Bivand on the r-help list for much of this code
##
library(maptools)
library(sp)
library(lattice)

## read in the (edited) map shape file
FS <- readShapePoly("/home/sasha/maps/rusobls.shp", IDvar="SUBJECT")

# Note, due to the length limits for strings in DBF files, the following
## subjects have different names in NOBUS and the shape file
## you MUST!!! rename all of them in NOBUS before proceeding
## NOBUS NAME				SHAPE FILE NAME
## kaliningradskaya oblast		kalinigradskaya obl
## republic bashkortostan		bashkortostan
## republic severnaya osetiya		severnaya osetia
## karachaevo-cherkesskaya republic	karachaevo republic
## kabardino-balkarskaya republic	kabardino republic
## republic saha (yakutiya)		republic saha

################ HEALTHY POPULATION #############################
## generate the mean of the variable in question

agg1 <- tapply(probworkinghealthymale, a001ter, mean, na.rm=TRUE)
agg2 <- data.frame(mean_probworkinghealthymale=agg1, row.names=names(agg1))

## extract FS data slot
FSd <- as(FS, "data.frame")

## merge the two datasets
FS1d <- merge(FSd, agg2, by="row.names", all=TRUE)

## do the same for the disabled (group 3) population
agg3 <- tapply(probworkinghealthyfemale, a001ter, mean, na.rm=TRUE)
agg4 <- data.frame(mean_probworkinghealthyfemale=agg3, row.names=names(agg3))
FS2d <- merge(FSd, agg4, by="row.names", all=TRUE)

## rename the rows
## STOP! IF THIS ERRORS OUT, YOU DID SOMETHING WRONG AND CANNOT PROCEED
## CHECK THAT ALL FS1d$SUBJECT ARE UNIQUE
row.names(FS1d) <- FS1d$SUBJECT
row.names(FS2d) <- FS2d$SUBJECT

## manually replace the following observations
FS1d["yamalo-nenetskiy ao","mean_probworkinghealthymale"]<-FS1d["tiumenskaya oblast", "mean_probworkinghealthymale"]
FS1d["khanti-mansi ao","mean_probworkinghealthymale"]<-FS1d["tiumenskaya oblast", "mean_probworkinghealthymale"]
FS1d["komi-permyak ao","mean_probworkinghealthymale"]<-FS1d["permskaya oblast", "mean_probworkinghealthymale"]
FS1d["nenetskiy ao","mean_probworkinghealthymale"]<-FS1d["arkhangelskaya oblast", "mean_probworkinghealthymale"]
FS1d["aginsky buryatsky ao","mean_probworkinghealthymale"]<-FS1d["republic buryatiya", "mean_probworkinghealthymale"]
FS1d["ust-ordynsky ao","mean_probworkinghealthymale"]<-FS1d["irkutskaya oblast", "mean_probworkinghealthymale"]
FS1d["koryakskiy ao","mean_probworkinghealthymale"]<-FS1d["kamchatskaya oblast", "mean_probworkinghealthymale"]
FS1d["evenkiyskiy ao","mean_probworkinghealthymale"]<-FS1d["krasnoyarskiy krai", "mean_probworkinghealthymale"]
FS1d["taimyrskiy ao","mean_probworkinghealthymale"]<-FS1d["krasnoyarskiy krai", "mean_probworkinghealthymale"]

## finally, there are some geographically disjoint areas
FS1d["nijegorodskaya obl","mean_probworkinghealthymale"]<-FS1d["nijegorodskaya oblast", "mean_probworkinghealthymale"]
FS1d["saratovskaya (east)","mean_probworkinghealthymale"]<-FS1d["saratovskaya oblast", "mean_probworkinghealthymale"]

## SANITY CHECK. THESE MUST MACH! IF THERE ARE NA'S, STOP, DO NOT PROCEED!!
match(row.names(FS1d), sapply(slot(FS, "polygons"), function(x) slot(x, "ID")))

## again, same thing for disabled
FS2d["yamalo-nenetskiy ao","mean_probworkinghealthyfemale"]<-FS2d["tiumenskaya oblast", "mean_probworkinghealthyfemale"]
FS2d["khanti-mansi ao","mean_probworkinghealthyfemale"]<-FS2d["tiumenskaya oblast", "mean_probworkinghealthyfemale"]
FS2d["komi-permyak ao","mean_probworkinghealthyfemale"]<-FS2d["permskaya oblast", "mean_probworkinghealthyfemale"]
FS2d["nenetskiy ao","mean_probworkinghealthyfemale"]<-FS2d["arkhangelskaya oblast", "mean_probworkinghealthyfemale"]
FS2d["aginsky buryatsky ao","mean_probworkinghealthyfemale"]<-FS2d["republic buryatiya", "mean_probworkinghealthyfemale"]
FS2d["ust-ordynsky ao","mean_probworkinghealthyfemale"]<-FS2d["irkutskaya oblast", "mean_probworkinghealthyfemale"]
FS2d["koryakskiy ao","mean_probworkinghealthyfemale"]<-FS2d["kamchatskaya oblast", "mean_probworkinghealthyfemale"]
FS2d["evenkiyskiy ao","mean_probworkinghealthyfemale"]<-FS2d["krasnoyarskiy krai", "mean_probworkinghealthyfemale"]
FS2d["taimyrskiy ao","mean_probworkinghealthyfemale"]<-FS2d["krasnoyarskiy krai", "mean_probworkinghealthyfemale"]

## finally, there are some geographically disjoint areas
FS2d["nijegorodskaya obl","mean_probworkinghealthyfemale"]<-FS2d["nijegorodskaya oblast", "mean_probworkinghealthyfemale"]
FS2d["saratovskaya (east)","mean_probworkinghealthyfemale"]<-FS2d["saratovskaya oblast", "mean_probworkinghealthyfemale"]

## Merge men and women
FS3d <- merge(FS1d, FS2d, by="row.names", all=TRUE)
row.names(FS3d) <- FS3d$SUBJECT.y

## MERGE INTO SPATIAL
FS1 <- SpatialPolygonsDataFrame(as(FS, "SpatialPolygons"), data=FS3d)

## plot
trellis.par.set(sp.theme()) # sets color ramp to bpy.colors()

## scale
scale = list("SpatialPolygonsRescale", layout.scale.bar(), which=2, offset=c(-3000000,3000000), fill=c("transparent", "black"), scale=1000000)
text1 = list("sp.text", c(-3000000, 3150000), "0", cex=.5, which=2)
text2 = list("sp.text", c(-2000000, 3150000), "1000 km", cex=.5, which=2)
text3 = list("sp.text", c(-2000000, 2850000), "Albers Equal-Area Projection", cex=.5, which=2)
text4 = list("sp.text", c(2200000, -1200000), "GIS data Copyright 1998 Alexander Perepechko and Dmitry Sharkov.", cex=.33, which=2)

spplot(FS1, c("mean_probworkinghealthymale", "mean_probworkinghealthyfemale"), names.attr=c("Men", "Women"), col.regions=rev(heat.colors(256)), colorkey=list(space="bottom"), main="Probability of Working by Subject of the Federation", sub="Healthy Population", as.table=TRUE, sp.layout = list(scale, text1, text2, text3, text4))

################ GROUP 3 DISABLED POPULATION #############################
## generate the mean of the variable in question

agg1 <- tapply(probworkinggroup3male, a001ter, mean, na.rm=TRUE)
agg2 <- data.frame(mean_probworkinghealthymale=agg1, row.names=names(agg1))

## extract FS data slot
FSd <- as(FS, "data.frame")

## merge the two datasets
FS1d <- merge(FSd, agg2, by="row.names", all=TRUE)

## do the same for the disabled (group 3) population
agg3 <- tapply(probworkinggroup3female, a001ter, mean, na.rm=TRUE)
agg4 <- data.frame(mean_probworkinghealthyfemale=agg3, row.names=names(agg3))
FS2d <- merge(FSd, agg4, by="row.names", all=TRUE)

## rename the rows
## STOP! IF THIS ERRORS OUT, YOU DID SOMETHING WRONG AND CANNOT PROCEED
## CHECK THAT ALL FS1d$SUBJECT ARE UNIQUE
row.names(FS1d) <- FS1d$SUBJECT
row.names(FS2d) <- FS2d$SUBJECT

## manually replace the following observations
FS1d["yamalo-nenetskiy ao","mean_probworkinghealthymale"]<-FS1d["tiumenskaya oblast", "mean_probworkinghealthymale"]
FS1d["khanti-mansi ao","mean_probworkinghealthymale"]<-FS1d["tiumenskaya oblast", "mean_probworkinghealthymale"]
FS1d["komi-permyak ao","mean_probworkinghealthymale"]<-FS1d["permskaya oblast", "mean_probworkinghealthymale"]
FS1d["nenetskiy ao","mean_probworkinghealthymale"]<-FS1d["arkhangelskaya oblast", "mean_probworkinghealthymale"]
FS1d["aginsky buryatsky ao","mean_probworkinghealthymale"]<-FS1d["republic buryatiya", "mean_probworkinghealthymale"]
FS1d["ust-ordynsky ao","mean_probworkinghealthymale"]<-FS1d["irkutskaya oblast", "mean_probworkinghealthymale"]
FS1d["koryakskiy ao","mean_probworkinghealthymale"]<-FS1d["kamchatskaya oblast", "mean_probworkinghealthymale"]
FS1d["evenkiyskiy ao","mean_probworkinghealthymale"]<-FS1d["krasnoyarskiy krai", "mean_probworkinghealthymale"]
FS1d["taimyrskiy ao","mean_probworkinghealthymale"]<-FS1d["krasnoyarskiy krai", "mean_probworkinghealthymale"]

## finally, there are some geographically disjoint areas
FS1d["nijegorodskaya obl","mean_probworkinghealthymale"]<-FS1d["nijegorodskaya oblast", "mean_probworkinghealthymale"]
FS1d["saratovskaya (east)","mean_probworkinghealthymale"]<-FS1d["saratovskaya oblast", "mean_probworkinghealthymale"]

## SANITY CHECK. THESE MUST MACH! IF THERE ARE NA'S, STOP, DO NOT PROCEED!!
match(row.names(FS1d), sapply(slot(FS, "polygons"), function(x) slot(x, "ID")))

## again, same thing for disabled
FS2d["irkutskaya oblast","mean_probworkinghealthyfemale"]<-0
FS2d["yamalo-nenetskiy ao","mean_probworkinghealthyfemale"]<-FS2d["tiumenskaya oblast", "mean_probworkinghealthyfemale"]
FS2d["khanti-mansi ao","mean_probworkinghealthyfemale"]<-FS2d["tiumenskaya oblast", "mean_probworkinghealthyfemale"]
FS2d["komi-permyak ao","mean_probworkinghealthyfemale"]<-FS2d["permskaya oblast", "mean_probworkinghealthyfemale"]
FS2d["nenetskiy ao","mean_probworkinghealthyfemale"]<-FS2d["arkhangelskaya oblast", "mean_probworkinghealthyfemale"]
FS2d["aginsky buryatsky ao","mean_probworkinghealthyfemale"]<-FS2d["republic buryatiya", "mean_probworkinghealthyfemale"]
FS2d["ust-ordynsky ao","mean_probworkinghealthyfemale"]<-FS2d["irkutskaya oblast", "mean_probworkinghealthyfemale"]
FS2d["koryakskiy ao","mean_probworkinghealthyfemale"]<-FS2d["kamchatskaya oblast", "mean_probworkinghealthyfemale"]
FS2d["evenkiyskiy ao","mean_probworkinghealthyfemale"]<-FS2d["krasnoyarskiy krai", "mean_probworkinghealthyfemale"]
FS2d["taimyrskiy ao","mean_probworkinghealthyfemale"]<-FS2d["krasnoyarskiy krai", "mean_probworkinghealthyfemale"]

## finally, there are some geographically disjoint areas
FS2d["nijegorodskaya obl","mean_probworkinghealthyfemale"]<-FS2d["nijegorodskaya oblast", "mean_probworkinghealthyfemale"]
FS2d["saratovskaya (east)","mean_probworkinghealthyfemale"]<-FS2d["saratovskaya oblast", "mean_probworkinghealthyfemale"]

## Merge men and women
FS3d <- merge(FS1d, FS2d, by="row.names", all=TRUE)
row.names(FS3d) <- FS3d$SUBJECT.y

## MERGE INTO SPATIAL
FS1 <- SpatialPolygonsDataFrame(as(FS, "SpatialPolygons"), data=FS3d)

## plot
trellis.par.set(sp.theme()) # sets color ramp to bpy.colors()

## scale
scale = list("SpatialPolygonsRescale", layout.scale.bar(), which=2, offset=c(-3000000,3000000), fill=c("transparent", "black"), scale=1000000)
text1 = list("sp.text", c(-3000000, 3150000), "0", cex=.5, which=2)
text2 = list("sp.text", c(-2000000, 3150000), "1000 km", cex=.5, which=2)
text3 = list("sp.text", c(-2000000, 2850000), "Albers Equal-Area Projection", cex=.5, which=2)
text4 = list("sp.text", c(2200000, -1200000), "GIS data Copyright 1998 Alexander Perepechko and Dmitry Sharkov.", cex=.33, which=2)

spplot(FS1, c("mean_probworkinghealthymale", "mean_probworkinghealthyfemale"), names.attr=c("Men", "Women"), col.regions=rev(heat.colors(256)), colorkey=list(space="bottom"), main="Probability of Working by Subject of the Federation", sub="Group III Disabled Population", as.table=TRUE, sp.layout = list(scale, text1, text2, text3, text4))

###  Disability incidence
 agg1 <- tapply(disabdummy, a001ter, mean, na.rm=TRUE)
# then do the above
# then
scale = list("SpatialPolygonsRescale", layout.scale.bar(), which=1, offset=c(-3000000,3000000), fill=c("transparent", "black"), scale=1000000)
text1 = list("sp.text", c(-3000000, 3150000), "0", cex=.75, which=1)
text2 = list("sp.text", c(-2000000, 3150000), "1000 km", cex=.75, which=1)
text3 = list("sp.text", c(-2000000, 2850000), "Albers Equal-Area Projection", cex=.75, which=1)
text4 = list("sp.text", c(2200000, -1200000), "GIS data Copyright 1998 Alexander Perepechko and Dmitry Sharkov.", cex=.5, which=1)

spplot(FS1, "mean_probworkinghealthymale", col.regions=rev(heat.colors(256)), main="Disability Incidence by Subject of the Federation", sub="Working Age Population", sp.layout = list(scale, text1, text2, text3, text4))