R/mrp.zurich.R

In MrPzurich2013: MrP functions for MA class Zurich 2013

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)

}