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:

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:

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

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

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