Model Thinking- Segregation and Peer Effects- Schelling Segregation- Experiments in R

Schelling.R

# here's my hacky code for a Schelling segregation model. 

# assumptions I've made in order to make life easier:
# 1) vacant cells have no influence. If a person happens to be totally surrounded by vacant cells, 
#    they don't move
# 2) when a person moves, they have no info about where they're going. This saves tons of computations
#    about the preferability of all the vacant plots
# 3) a person knows perfectly well if their threshold is met or not (threshold is here is the 
#    minimum proportion same)
# 4) all movers have the same chances of landing in any of the vacant plots.
# 5) if there are more people wanting to move than there are vacant plots, a random sample of 
#    potential movers is taken.
# 6) there are only 2 kinds of people, indicated with default colors blue and green, but that can be 
#    changed in the plot() function. In the matrix, these are indicated with TRUE and FALSE, for no 
#    reason other than computational efficiency. Obviously, if there were more groups moving around one 
#    would have to use a different data type (integer, character, factor), and these might slow things down.

# first a bunch of functions are defined. Lightly annotated. 

# There is room for improvement here. 
# [bug mentioned in previous version has been fixed- at this point I'm not getting any errors, but if 
# someone tries this out and notices one, please email me: tim.riffe@gmail.com]

# first internal functions ------------------------------------------------------------#
# given individual of position index 'ind', what are statuses of all neighbors?
queenNeighbors <- function(ind, M, n){
	inds <- ind+outer(c(-1, 0, 1), c(-n, 0, n), "+")
	M[inds[inds > 0 & inds <= length(M) & inds != ind]]
}

# determine if an individual of position index 'ind' will move
MoveYes <- function(ind, M, n, threshold = threshold){
	status <- M[ind]
	neighbors <- queenNeighbors(ind = ind, M = M, n = n)
	psame <- sum(neighbors == status, na.rm = TRUE) / length(neighbors[!is.na(neighbors)])
	if (!(is.null(psame)|is.na(psame))){
		if (psame < threshold) return(ind)
	}
}

# determine all who will move
Movers <- function(M, n, threshold = threshold){
	unlist(sapply((1:length(M))[!is.na(M)], MoveYes, M = M, threshold = threshold, n = 10))
}

# carries out a single iteration:
# 1) determine who will move
# 2) what are their statuses
# 3) where do they go (random assignment)
# 4) empty previous occupied spaces
# 5) if at some point no one wants to move, pass on a message to next function to
#    it knows to break the cycle
MoveCycle <- function(M, n = 10, threshold = threshold, nvacants = nvacants){
	movers.all.ind <- Movers(M, n = n, threshold = threshold)
	if (length(movers.all.ind) > 0){
		if (length(movers.all.ind) > nvacants){
			movers.all.ind <- sample(movers.all.ind, nvacants, replace = FALSE)
		}
		movers.all.statuses <- M[movers.all.ind]
		# random reassignment of places movers.all.statuses to available spots:
		M[sample((1:length(M))[is.na(M)], length(movers.all.ind),replace=FALSE)] <- movers.all.statuses
		M[movers.all.ind] <- NA
		return(M)
	} else {
		return(list("breaktheloop", M = M))
	}
}
#    end internal functions ------------------------------------------------------------#

# main function ------------------------------------------------------------------------#
# given starting matrix M0, a 'sameness' threshold and a maximum numer of iterations, 
# walk through the moving process. If at some point everyone is satisfied, then the loop
# breaks
SchellingProcess <- function(M0, threshold = .3, maxit = 100){
	nvacants <- sum(is.na(M0))
	M <- list()
	M[[1]] <- M0
	n <- ncol(M)
	for (i in 1:maxit){
		M[[i+1]] <- MoveCycle(M[[i]], n = n, threshold = threshold, nvacants = nvacants)
		if (!is.matrix(M[[i + 1]])){
			break
		}
	}
	out <- list(M = M, stable = !(i == maxit), threshold = threshold)
	class(out) <- "Schelling"
	return(out)
}
# end main function --------------------------------------------------------------------#

# define plot method -------------------------------------------------------------------#
# plot() method for output from SchellingProcess(), which is of class "Schelling"
plot.Schelling <- function(M, frame.pause = .1,col=c("darkblue", "green3")){
	iterations <- length(M$M)-1
	for (i in 1:iterations){
		Sys.sleep(frame.pause)
		image(t(M$M[[i]]), useRaster = TRUE, col = col, axes = FALSE)
	}
}
# end plot method ----------------------------------------------------------------------#

# function to make starting population -------------------------------------------------#
# just an easy way to get a starting matrix randomly filled in
PopulateMatrix <- function(nrow = 10, ncol = 10, prob = c(.3,.7), propvacant = .2){
	l <- nrow*ncol
	samplesize <- round(l * (1 - propvacant))
	M <- matrix(NA, nrow = nrow, ncol = ncol)
	M[sample(1:prod(dim(M)), samplesize, replace = FALSE)] <- sample(c(TRUE, FALSE), samplesize, prob = c(.3, .7), replace = TRUE)
	return(M)
}
# --------------------------------------------------------------------------------------#

# example usage ------------------------------------------------------------------------#
SchellingOutput <- SchellingProcess(PopulateMatrix(nrow=20,ncol=20, propvacant = .3),threshold = .35, maxit = 500)
SchellingOutput$stable
length(SchellingOutput$M)-1 # number of iterations
plot(SchellingOutput, frame.pause = 1)

# or to save a multipage pdf:
pdf("SchellingTest.pdf")
plot(SchellingOutput, frame.pause = 0)
dev.off()
# it's located in this directory:
getwd()