Nothing
mrp.zurich <-
function(response.model, augment.data=NA, augment.row=0, uncertainty=FALSE, Number.sim=1000, output=2){
ptm <- proc.time()
alarm <-0
aug <- 1
if(missing(augment.data)) aug <- 0
# read in packages
#library(lme4)
#library(blme)
#library(arm)
# read in RESPONSE model
model.response <- response.model
beta.fe <- fixef(model.response)
alpha.re <- ranef(model.response)
region.pointer <- c(4,2,6,6,6,6,6,5,6,2,2,3,3,5,5,5,5,5,3,5,7,1,1,2,1,2)
canton.pointer <- c(1:26)
alpha.var <- names(model.response@flist)
check.a <- c("age", "cantonnr", "education", "region","woman")==sort(alpha.var) #
S <- length(which(check.a==FALSE))
if (S!=0) stop(gettextf("The names of the random effects do not match! \n Please use 'age', 'cantonnr', 'education', 'region', and 'woman'.")) #
beta.var <- names(fixef(model.response))
# get the data for L2
dim22 <- dim(model.response@frame)[2]
dim21 <- dim22 - 4 # - 4
data.level2 <- model.response@frame[,-c(1,dim21:dim22)]
Knr <- model.response@frame$cantonnr
dimX <- 2
ifelse(dim(model.response@frame)[2]>7, dimX <- dim(data.level2)[2] + 1, dimX <- dimX) #>7
X2.26 <- matrix(NA,26,dimX)
# more than one X on L2
if (dim(model.response@frame)[2]>7){ #>7
for (ww in 1:26){
X2.26[ww,] <- c(1,colMeans(data.level2[Knr==ww,]))
}
}
# exactly one X on L2
if (dim(model.response@frame)[2]==7){ #==7
for (ww in 1:26){
X2.26[ww,] <- c(1,mean(data.level2[Knr==ww]))
}
}
# No X on L2
if (dim(model.response@frame)[2]==6){ #==6
X2.26 <- as.matrix(rep(1,26))
}
## Augment data step:
if (aug==1){
X2.26[augment.row,] <- augment.data
}
if (uncertainty==TRUE){
if (output>0){
cat("----------------------------------------------------\n You are requesting estimates AND uncertainty thereof \n")
cat("----------------------------------------------------\n")
}
if (Number.sim<25){
stop(gettextf("To estimate the uncertainty one needs a sufficient number of simulations. \n You chose a too low value for 'Number.sim'."))
}
# use SIM command to get N.sim draws of RE's & FE's
N.sim <- Number.sim
sim.block <- sim(model.response,N.sim)
educ.a.block <- ranef(sim.block)$education[,,1]
age.a.block <- ranef(sim.block)$age[,,1]
woman.a.block <- ranef(sim.block)$woman[,,1]
ctnr.a.block0 <- ranef(sim.block)$cantonnr[,,1]
ctnr.a.block <- matrix(NA,N.sim,26)
for (r in 1:26){
dd <- which(as.numeric(colnames(ctnr.a.block0))==r)
ctnr.a.block[,r] <- ctnr.a.block0[,dd]
}
region.a.block <- ranef(sim.block)$region[,,1]
beta.a.block <- fixef(sim.block)
# Help objects for the census information (educ6, age4, woman2)
education.pointer <- c(rep(1,8),rep(2,8),rep(3,8),rep(4,8),rep(5,8),rep(6,8)) # for each of the 48 SV
age.pointer <- rep(c(1,1,2,2,3,3,4,4),6) # ditto
woman.pointer <- rep(c(1,2),24)
check.vector <- matrix(NA,48,N.sim) # remember structure of census: educ (6), age (4), woman (2)
for (j in 1:N.sim){
for (i in 1:48){
check.vector[i,j] <- educ.a.block[j,education.pointer[i]] + age.a.block[j,age.pointer[i]] + woman.a.block[j,woman.pointer[i]]
}
}
# create individual level contribution on latent, 48*24
ylat.lev1.48.26 <- rbind(check.vector, check.vector, check.vector, check.vector, check.vector, check.vector, check.vector, check.vector, check.vector, check.vector, check.vector, check.vector, check.vector, check.vector, check.vector, check.vector, check.vector, check.vector, check.vector, check.vector, check.vector, check.vector, check.vector, check.vector, check.vector, check.vector)
# create ceontext level contribution on latent, 48*24
ylat.lev2 <- X2.26 %*% t(beta.a.block) + t(ctnr.a.block) + t(region.a.block)[region.pointer,]
ylat.lev2.48.26 <- matrix(NA, 1248, N.sim)
for (zz in 1:26){
aa <- (zz-1)*48 + 1
bb <- zz*48
for (jj in aa:bb){
ylat.lev2.48.26[jj,] <- ylat.lev2[zz,]
}
}
y.lat <- ylat.lev2.48.26 + ylat.lev1.48.26
# read-in census data
#data("Census48.rda")
census48 <- census$census48.MAZH2013
colStandard <- function(x){
x <- x/sum(x)
}
weight.48.26 <- apply(census48,2,colStandard)
# generate predictions
phat.weighted <- matrix(NA, 1248, N.sim)
for (ss in 1:N.sim){
phat.weighted[,ss] <- pnorm(y.lat[,ss]) * weight.48.26
}
# prediction per canton:
phat.canton <- matrix(NA,26,N.sim)
for (nn in 1:N.sim){
for (tt in 1:26){
aa <- (tt-1)*48 + 1
bb <- tt*48
phat.canton[tt, nn] <- sum(phat.weighted[aa:bb, nn])
}
}
# Estimate is mean of all simulations
preport <- 100*round(rowMeans(phat.canton),3)
preport <- paste(preport,rep("%",26), sep=" ")
testi <- rowMeans(phat.canton)
if (NaN %in% testi) {
alarm <- 1
missing.canton <- which(is.na(testi))
}
# end of uncertainty==TRUE
}
if (uncertainty==FALSE){
if (output>0){
cat("----------------------------------------------------\n No uncertainty measures requested \n")
cat("----------------------------------------------------\n")
}
# Read out coefs
educ.a.block <- alpha.re$education[,,1]
age.a.block <- alpha.re$age[,,1]
woman.a.block <- alpha.re$woman[,,1]
ctnr.a.block0 <- alpha.re$cantonnr[,,1]
ctnr.a.block <- rep(0,26)
keeper <- which(as.numeric(c(1:26) %in% response.model@frame$cantonnr)==1)
a <- 0
for (r in 1:26){
if (r %in% keeper){
ctnr.a.block[r] <- ctnr.a.block0[r-a]
}
if ((r %in% keeper)==FALSE){
ctnr.a.block[r] <- 0
a <- a+1
}
}
region.a.block <- alpha.re$region[,,1]
beta.a.block <- beta.fe
# Help objects for the census information (educ6, age4, woman2)
education.pointer <- c(rep(1,8),rep(2,8),rep(3,8),rep(4,8),rep(5,8),rep(6,8)) # for each of the 48 SV
age.pointer <- rep(c(1,1,2,2,3,3,4,4),6) # ditto
woman.pointer <- rep(c(1,2),24)
check.vector <- rep(NA,48) # remember structure of census: educ (6), age (4), woman (2)
for (i in 1:48){
check.vector[i] <- educ.a.block[education.pointer[i]] + age.a.block[age.pointer[i]] + woman.a.block[woman.pointer[i]]
}
# create individual level contribution on latent, 48*26
ylat.lev1.48.26 <- rep(check.vector, 26)
# create ceontext level contribution on latent, 48*24
#print(region.a.block)
ylat.lev2 <- X2.26 %*% as.matrix(beta.a.block) + ctnr.a.block + region.a.block[region.pointer]
# ylat.lev2 <- X2.26 %*% beta.a.block
ylat.lev2.48.26 <- rep(NA, 1248)
for (zz in 1:26){
aa <- (zz-1)*48 + 1
bb <- zz*48
for (jj in aa:bb){
ylat.lev2.48.26[jj] <- ylat.lev2[zz]
}
}
y.lat <- ylat.lev2.48.26 + ylat.lev1.48.26
# read-in census data
#data("Census48.rda")
census48 <- census$census48.MAZH2013
colStandard <- function(x){
x <- x/sum(x)
}
weight.48.26 <- apply(census48,2,colStandard)
# generate predictions
phat.weighted <- pnorm(y.lat) * weight.48.26
# prediction per canton:
phat.canton <- rep(NA,26)
for (tt in 1:26){
aa <- (tt-1)*48 + 1
bb <- tt*48
phat.canton[tt] <- sum(phat.weighted[aa:bb])
}
preport <- 100*round(phat.canton,3)
preport <- paste(preport,rep("%",26), sep=" ")
testi <- phat.canton
if (NaN %in% testi) {
alarm <- 1
missing.canton <- which(is.na(testi))
}
# end of uncertainty==TRUE
}
cantons.name <- c("ZH", "BE", "LU", "UR", "SZ", "OW", "NW", "GL", "ZG", "FR", "SO", "BS", "BL", "SH", "AR", "AI", "SG", "GR", "AG", "TG", "TI", "VD", "VS", "NE", "GE", "JU")
reportL <- data.frame(cbind(cantons.name,preport))
colnames(reportL) <- c("Canton", "Prediction")
if (output>1) print(reportL)
if (alarm==1){
N.missing <- length(as.vector(missing.canton))
TEXT <- paste(cantons.name[missing.canton][1:N.missing],sep=" ", collapse=" / ")
print(TEXT)
writeLines(paste("Your survey model's data not contain any respondents from canton(s):", TEXT, "\n(see option 'augment.data')", sep=" "))
#print(cantons.name)
#print(missing.canton)
}
time1 <- proc.time() - ptm
if (output>0){
cat("\n----------------------------------------------------\n ")
writeLines((paste(" Estimation time:",round(time1[3],3),"seconds.", sep=" ")))
# cat("----------------------------------------------------\n (® Leemann, 2013) \n")
writeLines(paste("----------------------------------------------------"))
}
return(phat.canton)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.