Model Thinking- Segregation and Peer Effects- The Standing Ovation Model (basic)- Experiments in R

StandingOvation.R

# Author: Tim Riffe: tim.riffe@gmail.com

# A basic implementation of the standing ovation model in R.
# The first part walks through the components, explaining the 
# inputs and what's going on. It's a full example, start to finish
# this is follwed by a function StandingOvationSimple(), that just takes 3 inputs 
# (with some other optional parameters) and will animate the process for you
# you may wish to modify the function to return a series of statistics taken
# along the progression of the model. Enjoy!


# Part I
# walk through an example:
# ------------------------------------------------------------------------------
# standing ovation with perfect visibility and no weights:

# for this code, you need to use the fields() package, run:
# install.packages("fields") # if necessary

# first define a matrix of thresholds for minimum quality of show that it would take to stand.
# assume rectangular theater.
set.seed(1)

# each person has a minimum show quality required in order to stand
# here, we randomly assign a threshold to everyone in the audience
# using a uniform distribution between .5 and 1.
Thresholds <- matrix(runif(200,min=.5,max=1),ncol=20,nrow=10)
par(mar=c(2,2,2,2))
fields:::image.plot(x=1:20,y=1:10,t(Thresholds),asp=1,axes=FALSE,xlab="",ylab="",zlim=c(0,1))
rect(5,13,15,18,col="black")
text(10,15,"STAGE",col="white",xpd=TRUE,cex=2)
text(10,-3,"Each person's threshold for standing",xpd=TRUE)
rect(.5,.5,20.5,10.5)

# likewise, not everyone perceives the show with the same precision
# (sophistication). There is some error. I'm going to assume some error
# using a normal distribution with sd of .1
Errors <- matrix(rnorm(200,sd=.1),ncol=20,nrow=10)
par(mar=c(2,2,2,2))
fields:::image.plot(x=1:20,y=1:10,t(Errors),asp=1,axes=FALSE,xlab="",ylab="",zlim=c(-.3,.3))
rect(5,13,15,18,col="black")
text(10,15,"STAGE",col="white",xpd=TRUE,cex=2)
text(10,-3,"Each person's error for judging show",xpd=TRUE)
rect(.5,.5,20.5,10.5)

# so really, the show will have a particular quality, which we treat as
# a fixed deterministic value. Let's say it's a quality of .7
Quality <- .4
# so therefore the perceived quality will be .7, plus each person's error:
Perception <- Quality + Errors
par(mar=c(2,2,2,2))
fields:::image.plot(x=1:20,y=1:10,t(Perception),asp=1,axes=FALSE,xlab="",ylab="",zlim=c(0,1))
rect(5,13,15,18,col="black")
text(10,15,"STAGE",col="white",xpd=TRUE,cex=2)
text(10,-3,"Perceived Quality of Show",xpd=TRUE)
rect(.5,.5,20.5,10.5)

# now, if perceived quality is higher than one's threshold for standing, one will stand
# without anyone else standing. These are potential 'first-standers'
# StandFirst is a logical matrix. Either you're standing or you're not...
StandFirst <- Perception > Thresholds
par(mar=c(2,2,2,2))
fields:::image.plot(x=1:20,y=1:10,t(StandFirst),asp=1,axes=FALSE,xlab="",ylab="",col=c("white","red"))
rect(5,13,15,18,col="black")
text(10,15,"STAGE",col="white",xpd=TRUE,cex=2)
text(10,-3,"The people in seats marked red stand first",xpd=TRUE)
rect(.5,.5,20.5,10.5)

# However, everyone is influenced by everyone else as well. If the ball gets rolling
# then everyone might stand. Each person will also have some minimum threshold for the
# proportion of others standing that it would take in order to also join and stand
# right now the proportion standing is:
sum(StandFirst)/length(StandFirst) # possibly not that high
# assign everyone a threshold, where if surpassed they will stand. Let's call it
# banwagonry, distributed between 1/size of public to 1 (perfect stubbornness)
Bandwagonry <- matrix(runif(200,min=(1/200),max=1),ncol=20,nrow=10)
# and really if you're already standing then we're not interested in you (seats in white)
Bandwagonry[StandFirst] <- NA
# take a look:
fields:::image.plot(x=1:20,y=1:10,t(Bandwagonry),asp=1,axes=FALSE,xlab="",ylab="",zlim=c(0,1))
rect(5,13,15,18,col="black")
text(10,15,"STAGE",col="white",xpd=TRUE,cex=2)
text(10,-3,"low numbers mean you're quick to follow\nhigh numbers mean you're self-assured",xpd=TRUE)
rect(.5,.5,20.5,10.5)

# at this point, let's assume everyone sees everyone else- it was a high energy and interactive show...
# we iterate once:
StandingNow <- StandFirst
# copy standers:
CopyStanders <- (sum(StandFirst)/length(StandFirst)) > Bandwagonry
fields:::image.plot(x=1:20,y=1:10,t(CopyStanders),asp=1,axes=FALSE,xlab="",ylab="",col=c("white","red"))
rect(5,13,15,18,col="black")
text(10,15,"STAGE",col="white",xpd=TRUE,cex=2)
text(10,-3,"the first to stand 'due to' others standing (red)",xpd=TRUE)
rect(.5,.5,20.5,10.5)

# which means that now the following people are standing:
Standers <- StandFirst
Standers[CopyStanders] <- TRUE
fields:::image.plot(x=1:20,y=1:10,t(Standers),asp=1,axes=FALSE,xlab="",ylab="",col=c("white","red"))
rect(5,13,15,18,col="black")
text(10,15,"STAGE",col="white",xpd=TRUE,cex=2)
text(10,-3,"these people are standing at one iteration in",xpd=TRUE)
rect(.5,.5,20.5,10.5)

# now we iterate, to see how far this will go
for (i in 1:50){
	p <- sum(Standers)/length(Standers)
	Bandwagonry[Standers] <- NA
	CopyStanders <- p > Bandwagonry
	if (sum(CopyStanders,na.rm=TRUE)==0){
					cat("whoever isn't standing at this point stays sitting\n")
					break
	}
	Standers[CopyStanders] <- TRUE
	fields:::image.plot(x=1:20,y=1:10,t(Standers),asp=1,axes=FALSE,xlab="",ylab="",col=c("white","red"))
	rect(5,13,15,18,col="black")
	text(10,15,"STAGE",col="white",xpd=TRUE,cex=2)
	text(10,-3,paste("people standing at time", i),xpd=TRUE)
	rect(.5,.5,20.5,10.5)
	Sys.sleep(2)
}
# final proportion standing:
p

# possible extensions:
# 1) assign people to continguous groups that have within-group influence
# 2) assign people visibility weights

# Part II
# a function to simply animate a process, given some inputs
# ------------------------------------------------------------------------------
# Now Here's a function to animate an entire standing ovation sequence, using the above-described
# objects. As inputs, it needs
# 1) a matrix of Thresholds
# 2) a matrix of Errors of the same dimensions
# 3) a matrix of Bandwagonry
# 4) a 'true' show quality.

# that way you can do experiments by simply changine the inputs

StandingOvationSimple <- function(Thresholds, Errors, Bandwagonry, True.Quality = .5, frame.pause = 2, maxit = 50){
	Perception <- True.Quality + Errors
	Standers <- Perception > Thresholds
	
	# initial plot
	fields:::image.plot(x = 1:ncol(Standers), y = 1:nrow(Standers), t(Standers), asp = 1, 
			axes = FALSE, xlab = "", ylab = "", col = c("white", "red"))
	rect(.5, .5, ncol(Standers) + .5, nrow(Standers) + .5)
	p <- sum(Standers)/length(Standers)
	title(paste(round(p*100),"% standing at time 0",sep=""))
	Sys.sleep(frame.pause)
	for (i in 1:maxit){
		p <- sum(Standers)/length(Standers)
		Bandwagonry[Standers] <- NA
		CopyStanders <- p > Bandwagonry
		if (sum(CopyStanders,na.rm=TRUE)==0){
			cat("whoever isn't standing at this point stays sitting\n")
			break
		}
		Standers[CopyStanders] <- TRUE
		fields:::image.plot(x = 1:ncol(Standers), y = 1:nrow(Standers), t(Standers), asp = 1, 
				axes = FALSE, xlab = "", ylab = "", col = c("white", "red"))
		rect(.5, .5, ncol(Standers) + .5, nrow(Standers) + .5)
		p <- sum(Standers)/length(Standers)
		title(paste(round(p*100),"% standing at time ",i,sep=""))
		Sys.sleep(frame.pause)
	}
}

# so, another test, just make the basic inputs:
set.seed(123)
# lets make higher avg threshold, but greater error too
Thresholds <- matrix(runif(200,min=.6,max=1),ncol=20,nrow=10)
Errors <- matrix(rnorm(200,sd=.15),ncol=20,nrow=10)
# Bandwagonry same distribution:
Bandwagonry <- matrix(runif(200,min=(1/200),max=1),ncol=20,nrow=10)

# this function does all the work:
StandingOvationSimple(Thresholds,Errors,Bandwagonry,True.Quality = .55, frame.pause = 1)

# you can change the size of the matrix too and it should still work the same.
# the dimensions won't make much difference unless you incorporate space/visibility/influence
# into the model