Building a Schelling Segregation Model in R

Happy New Year! Last year, our good friend Shawn over at Electric Archaeology introduced us to an excellent animated, interactive representation of Thomas Schelling’s (1969) model of segregation called “Parable of the Polygons”. I’ve always liked Schelling’s model because I think it illustrates the concepts of self-organization and emergence, and is also easy to explain, so it works as a useful example of a complex system. In the model, individuals situated in a gridded space decide whether to stay put or move based on a preference for neighbours like them. The model demonstrates how features of segregated neighborhoods can emerge even when groups are relatively ‘tolerant’ in their preferences for neighbors.

Here, I’ve created a simple version of Schelling’s model using R (building on Marco Smolla’s excellent work on creating agent-based models in R). Schelling’s model is situated on a grid, and in its simplest form, the cells of the grid will be in one of three states: uninhabited, inhabited by a member of one group, or inhabited by a member of a second group. This could be represented as a matrix of numbers, with each element being either 0, 1, or 2. So we’ll bring together these components as follows:

number<-2000
group<-c(rep(0,(51*51)-number),rep(1,number/2),rep(2,number/2))
grid<-matrix(sample(group,2601,replace=F), ncol=51)

par(mfrow=c(1,2))
image(grid,col=c("black","red","green"),axes=F)
plot(runif(100,0,1),ylab="percent happy",xlab="time",col="white",ylim=c(0,1))

Here, we start with a 51 x 51 grid of 2000 occupied cells.  To create this, a vector called group is generated that contains 1000 1s, 1000 2s, and the remainder are 0s. These are collated into a matrix called grid through random sampling of the group vector. Finally, the matrix is plotted as an image where occupied cells are colored green or red depending on their number while unoccupied cells are colored black, like so:

Rplot01
The next step is to establish the common preference for like neighbors, and to setup a variable which tracks the overall happiness of the population.

alike_preference<-0.60
happiness_tracker<-c()

Finally, we’ll need a function which we’ll use later to establish who the neighbors are for a given patch, which we will call get_neighbors. To do this, we’ll feed the function a set of two xy-coordinates as a vector (e.g. 2 13 ), and using a for-loop, pull each neighbor with the Moore neighborhood (8 surrounding patches) in order, counterclockwise from the right. Then we’ll need to ensure that if a neighboring cell goes beyond the bounds of the grid (>51 or <1), we account for this by grabbing the cell at the opposite end of the grid. This function will return eight pairs of coordinates as a matrix.

get_neighbors<-function(coords) {
  n<-c()
  for (i in c(1:8)) {
 
    if (i == 1) {
      x<-coords[1] + 1
      y<-coords[2]
    }

    if (i == 2) {
      x<-coords[1] + 1
      y<-coords[2] + 1
    }
  
    if (i == 3) {
      x<-coords[1]
      y<-coords[2] + 1
    }
    
    if (i == 4) {
      x<-coords[1] - 1
      y<-coords[2] + 1
    }
    
    if (i == 5) {
      x<-coords[1] - 1
      y<-coords[2]
    }
    
    if (i == 6) {
      x<-coords[1] - 1
      y<-coords[2] - 1
    }
   
    if (i == 7) {
      x<-coords[1]
      y<-coords[2] - 1
    }
    
    if (i == 8) {
      x<-coords[1] + 1
      y<-coords[2] - 1
    }
   
    if (x < 1) {
      x<-51
    }
    if (x > 51) {
      x<-1
    }
    if (y < 1) {
      y<-51
    }
    if (y > 51) {
      y<-1
    }
    n<-rbind(n,c(x,y))
  }
  n
}

Now to get into the program. We’ll run the process 1000 times to get output, so the whole thing will will be embedded in a for loop. Then, we’ll set up some variables which keep track of happy versus unhappy cells:

for (t in c(1:1000)) {
happy_cells<-c()
unhappy_cells<-c()  

Each of these tracker vectors (happy_cells and unhappy_cells) will keep track of additional vectors that contain coordinates of cells that are happy or unhappy.

Next, we’ll use two for loops to iterate through each row and column in the matrix. For each cell (here called current), we’ll take its value (0, 1, or 2). If the cell is not empty (that is, does not have a value of 0), then we’ll create variables that keep track of the number of like neighbors and the total number of neighbors (that is, neighboring cells which are inhabited), and then we’ll use our get_neighbors function to generate a vector called neighbors . Then we’ll use a for loop to iterate through each of those neighbors, and compare their values to the value of the current patch. If it is a match, we add 1 to the number of like neighbors. If it is inhabited, we add 1 to the total number of neighbors (a varable called all_neighbors). Then we divide the number of like neighbors by the total number of neighbors, and compare that number to the like-neighbor preference to determine whether the current patch is happy or not (The is.nan function is used here to escape situations where a cell is completely isolated and thus would involve division by 0). Happy patches are added to our happy patches variable, while unhappy patches are added to our unhappy patches variable, both as matrices of coordinates.

for (j in c(1:51)) {
  for (k in c(1:51)) {
    current<-c(j,k)
    value<-grid[j,k] 
    if (value > 0) {
      like_neighbors<-0
      all_neighbors<-0
      neighbors<-get_neighbors(current)
      for (i in c(1:nrow(neighbors))){
        x<-neighbors[i,1]
        y<-neighbors[i,2]
        if (grid[x,y] > 0) {
          all_neighbors<-all_neighbors + 1
        }
        if (grid[x,y] == value) {
          like_neighbors<-like_neighbors + 1
        }
      }
      if (is.nan(like_neighbors / all_neighbors)==FALSE) {
        if ((like_neighbors / all_neighbors) < alike_preference) {
            unhappy_cells<-rbind(unhappy_cells,c(current[1],current[2]))
        }
          else {
            happy_cells<-rbind(happy_cells,c(current[1],current[2]))
          }
        }
   
      else {
        happy_cells<-rbind(happy_cells,c(current[1],current[2]))
      }
    }
  }
}

Next, we’ll get our overall happiness by dividing the number of happy cells by the total number of occupied cells, and update our happiness tracker by appending that value to the end of the vector.

happiness_tracker<-append(happiness_tracker,length(happy_cells)/(length(happy_cells) + length(unhappy_cells)))

Next, we’ll get our unhappy patches to move to unoccupied spaces. To do this, we’ll randomly sample unhappy cells so we’re not introducing a spatial bias. Then, we’ll iterate through that sample, calling each patch in that group a mover, and picking a random spot in the grid as a moveto. A while loop will continue to pick a new random moveto while the current moveto is inhabited. Once an uninhabited moveto has been found, the mover’s value is applied to that patch, and removed from the original mover patch.

rand<-sample(nrow(unhappy_cells))
for (i in rand) {
  mover<-unhappy_cells[i,]
  mover_val<-grid[mover[1],mover[2]]
  move_to<-c(sample(1:51,1),sample(1:51,1))
  move_to_val<-grid[move_to[1],move_to[2]]
  while (move_to_val > 0 ){
    move_to<-c(sample(1:51,1),sample(1:51,1))
    move_to_val<-grid[move_to[1],move_to[2]]
  }
  grid[mover[1],mover[2]]<-0
  grid[move_to[1],move_to[2]]<-mover_val
}

Finally, we’ll check the output.

par(mfrow=c(1,2))
image(grid,col=c("black","red","green"),axes=F)
plot(runif(100,0,1),ylab="percent happy",xlab="time",col="white",ylim=c(0,1))
lines(happiness_tracker,oma = c(0, 0, 2, 0),col="red")
}

With the for loop we created around the whole program, we get animation in our graphical display. Here’s what we get when we set the alike_preference value to 70 percent:animation1

And here’s what happens when it’s set to 72 percent:

animation2

Finally, here’s what happens when it’s set to 80:

animation3

For comparison, check out this version in NetLogo or this HTML version.

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s