Simulating Social Network Hotspots

# This simulation will generate a simulated data set with social network connections.

# Unlike the previous post which made all connections randomly this post will have hotspots such as schools and common areas that networkers access which make it more likely for them to become part of the same social network.

# For reference on how to us R to analyze social networks, check out: Mike Nowak and Sean Westwood. 2010. "Social Network Analysis Labs in R." Stanford University.

# First lets decide how many people we would like to simulate in our social network.

npeople <- 50

# We will create an edgelist which defines all of the relationships between everybody in the network.

# First let's create a pure random network.

# First we will populate it by a list of all of the edges.
rnet <- data.frame(ego=rep(1:npeople,each=npeople), alter=rep(1:npeople,times=npeople), friendship=0, friendships=0)

# Now let's add three hotspots.  These hotspots will be randomly assigned to inviduals.

rnet$hot1 <- runif(npeople)
# Perhaps, physical location.
rnet$hot2 <- runif(npeople)
# Age group perhaps.
rnet$hot3 <- runif(npeople)
# Soioeconomic background.

# I will use a scaled version of the equclidean distance of the three hotspots to calculate the probability that any two people are connected.
euc.scale = 4

# Unlike the previous post which allowed anybody to be within a network with anybody else the hotspots will restrict social networks to people who are near each other in all three hotspots.
for (i in 1:npeople) for (ii in (i+1):npeople) if (i!=ii) {
  # Calculate Euclidean distance between person i and ii.
  distance <- ((rnet$hot1[i]-rnet$hot1[ii])^2+
               (rnet$hot2[i]-rnet$hot2[ii])^2+
               (rnet$hot2[i]-rnet$hot2[ii])^2)^.5
  # Make sure the probability of two people being friends is always positive.
  p.friend <- max(1-distance*euc.scale,0)
 
  # Now draw one draw using the probabily.
  if (rbinom(1,1,p.friend)==1) {
    rnet$friendship[(rnet$ego==i & rnet$alter==ii)] <- 1
    rnet$friendship[(rnet$ego==ii & rnet$alter==i)] <- 1
   
    rnet$friendships[rnet$ego==i] <- rnet$friendships[rnet$ego==i] + 1
   # rnet$friendships[rnet$ego==ii] <- rnet$friendships[rnet$ego==ii] + 1
   
  }
}

head(rnet)

# We can see that as we had planned there is a specific number of connections.
summary(rnet)

summary(rnet$friendships)

# For graphing our connections we will use the package igraph
require(igraph)

# Let's try generating our first graph:
g = graph.data.frame(rnet[rnet$friendship==1,], directed=F)
plot(g, vertex.size=(rnet$friendships[rnet$friendship==1]+2)*2)
# I was hoping to have the scale of the dots to be proportional to the numnber of connections.  However, clearly this is not yet working.  Something to figure out for later posts.

Simulating Social Network Data

# Code updated: Sorry about that sometimes blogger seems to go crazy with the code thinking it is html or something.

# This simulation will generate a simulated data set with social network connections.

# For reference on how to us R to analyze social networks, check out: Mike Nowak and Sean Westwood. 2010. "Social Network Analysis Labs in R." Stanford University.

# First lets decide how many people we would like to simulate in our social network.

npeople = 50
# We will create an edgelist which defines all of the relationships between everybody in the network.

# First let's create a pure random network.

# First we will populate it by a list of all of the edges.

rnet <- data.frame(ego=rep(1:npeople,each=npeople), alter=rep(1:npeople,times=npeople), friendship=0)

# Right now we have all of the edges (npeople^2) and now we just need to populate it with connections.

# These connection will be random on our first pass at this.

# Because these connections are "friendship" connections we will assume that they go both directions.

# Let's define the density of connections:  Ie number of connections per potential connection.

conDen = .05

# This double loop should make it so that every connection is examined for a potential connection

for (i in 1:npeople) for (ii in (i+1):npeople) if ((rbinom(1,1,conDen)==1)&(i!=ii)) {
  print(paste(i,ii))
  rnet$friendship[(rnet$ego==i & rnet$alter==ii)] <- 1
  rnet$friendship[(rnet$ego==ii & rnet$alter==i)] <- 1
}

head(rnet)

# We can see that as we had planned there is a specific number of connections.
summary(rnet)

# For graphing our connections we will use the package igraph
require(igraph)

# Let's try generating our first graph:
plot(graph.data.frame(rnet[rnet$friendship==1,], directed=F), main="Purely Random Connections")

# I reran this code a number of times before I found a network that did not have any stray members who were unconnected.  However, a large portion of the time at least one member us unconnected with the larger network.

# In later posts I will enrich this simulation to allow for networks to be generated dynamically.  That is friends bridge friendships between friends.

Generalized Graded Response Data Generating Command

# The following code will draw random test items from a graded response model.

# The command mat.binom will be useful in drawing a matrix of binomial results. I will use a command I programmed in a previous post. (http://www.econometricsbysimulation.com/2012/09/item-response-theory-estimation.html)
mat.binom <- function="function" n="n" p="p">  bin.mat <- p="p">  for (i in 1:nrow(p)) {
    for (ii in 1:ncol(p)) {
    # This will draw a random binomial for each of the probabilities.
    bin.mat[i,ii] <- i="i" ii="ii" n="n" p="p" rbinom="rbinom">  }
  }
 return(bin.mat)
}


rgrm <- d="1.7)" function="function" p="p" theta="0,a=cbind(rep(1,5),rep(1,5)),b=cbind(rep(0,5),rep(1,5)),">  # b is now an input matrix with each row representing a different item and each column representing a different grade for that item.  The number of columns should be equal to the maximum number of grades for all of the items.  Items may have less than the full number of grades indicated by NA values at upper levels.

  # We will use the grm Cululative Grade Function defined in a previous post as what we will use to generate our probabilities.
  # http://www.econometricsbysimulation.com/2012/09/generalize-graded-response-model.html
  grm <- a="1," b="c(0,1,2)," cplot="T," function="function" p="p" pplot="T," stackpoly="F" theta="1,">    ngrade = max(length(b),length(a))
    CGF = matrix(NA, ncol=ngrade  , nrow=length(theta))
    for (i in 1:length(theta)) CGF[i,] = exp(a*(theta[i]-b))/(1+exp(a*(theta[i]-b)))
    return(CGF)
  }

  # This will be the number of items to generate.
  if (sum(dim(a)!=dim(b))>0) warning()
  nitems = nrow(b)

  # This matrix will hold the results of the items.
  Y <- matrix="matrix" ncol="nitems)</p" nrow="length(theta),">
  rownames(Y) = paste("Stud", 1:length(theta))
  colnames(Y) = paste("Item", 1:nitems)

  for (i in 1:nitems) {
    # Draw the submatrix of a and b not equal to NA
    a.sub = !is.na(a[i,])
    b.sub = !is.na(b[i,])
    # Draw a uniform draw for every individual
    unif.draw = runif(length(theta))
    # Spread those draws out into a matrix for each grade of each item for each individual
    unif.draws = matrix(unif.draw, ncol=length(a.sub), nrow=length(theta))
    # Draw the probability of getting each grade for each item
    itemi = grm(theta=theta, a=a[i,a.sub], b=b[i,b.sub])
    # Calculate the grade for each item.
    Y[,i]<-apply itemi="itemi">unif.draw,1,sum)
    # Since the grm function creates a matrix with the probability of getting that value or more for each grade we can simply count the number of times that the random uniform draw got below that value.  If the draw is low enough then the student gets full credit.  If it is high enough then the student gets no credit.
  }
  # Specify the command's return value
  return(Y)
}

# Notice, theta can be any length but a and b must have the same dimensions
a=cbind(rep(1,6),rep(1,6),rep(c(NA,1),3))
b=cbind(rep(1,6),rep(2,6),rep(c(NA,3),3))
# This will create six items with grades from 1 to 2 and 1 to 2 to 3.  The NA's mean that grade is missing (undefined).

# The length of theta is equal to the number of students
theta=seq(0,7,.5)

rgrm(theta=theta,a=a , b=b)

Value-added modelling - Stata simulation - iPad example

* Value-added modelling is a common approach to use to try to infer the "quality" or "value" of inputs.

* In education, these methods have become quite popular politically and academically.

* But in a sense using these methods in education is an abstraction.

* Let us first in order to grasp how value added methods were first developed think of a production example.

* Imagine that you are manufacturing tablet computers (iPad).

* On each step of the production process there is a different "value" that is added to the process as a result the particular inputs.

* In order to conceptualize this think of the product at each stage being worth a certain value that is then bought buy the next person in the manufacturing chain.

* However, you do not want to use the market to assemble the tablet computers.

* You want to assemble them in house.

* Therefore you want to figure out a way of inferring how much value each input has.

* To do this, imagine that run a series of non-market valuations after each stage of the production process to infer a current price.

* Then you take the change in that market price as a measure of the value of those inputs.

set seed 11

* Let's do this:

clear
set obs 10

* 10 different companies

gen comp_id = mod(_n-1,10)+1
  label var comp_id "Company ID"

* Each of the companies has a different assembly line structure.
* Let's imagine that at each stage that company's structure ads a constant amount of value to all products.

gen comp_fe = runiform()
  label var comp_fe "The fixed effect (specific to that company) added to each product each stage"

* Now let's imagine that each company produces 10 different products (over the sample time)
expand 10

sort comp_id

* Generate a list of product ids
gen prod_id = _n
  label var prod_id "Product ID"

* Each product line has some inherent design component that makes it more or less valued at each progressive stage than other products.
gen prod_fe = rnormal()/2 + 1/4
  label var prod_fe "Product Fixed Effect"

* Now imagine that there are 5 stages of production for each product
expand 5
bysort prod_id: gen prod_stage = _n

* You have different stages with the initial stage product idea.
* Stage 1 gather raw materials
* Stage 2 manufacture components
* Stage 3 assemble components (probably in China)
* Stage 4 ship product
* Stage 5 sell product at the retail locations



label var prod_stage "Production stage"

* Now imagine also that there are 100 contractors in tablet computer production market.

* These contractors get random contracts as to which product to work on.

* This is what we really want to know.

* How good are these contractors at "adding value" to the product.

* This is the trickiest part of the code so far.

* First we need to generate the contractors.

* Keep the data that we have generated so far:
preserve

* There are several ways of doing this.

* I will use the many:1 merge command to accomplish this task.
clear

* Imagine that our 100 contractors are subsidiaries of 5 different umbrella companies.

set obs 5

gen cont_company_id=_n
  label var cont_company_id "Contracting company ID"

* Each of these companies has a different work ethic
gen cont_company_fe = rnormal()*.25
  label var cont_company_fe "Contracting company effectiveness"

* Each of the contracting companies has 20 contractors they manage.

expand 20

gen cont_id = _n
  label var cont_id "Contractor ID"

gen cont_fe = rnormal() + 1
  label var cont_fe "Contractor effectiveness"

* Now we will save the contractor information to a temporary data file
save "contractor.dta", replace

restore

* Let us first assign contractor IDs randomly:

gen cont_id=int(runiform()*100+1)

* Now let's merger in the contractor data

merge m:1 cont_id using "contractor.dta"
drop _merge

* Let us think that whenever a different product is developed there is some unobserved component that adds or subtracts random value from a product line independent of all other inputs at each stage.

gen rand_effect = rnormal()
  label var rand_effect "Random idiosyncratic production effect unique to each product at each stage"
* In other words the error component

* Finally, imagine that at each stage there is an "average" amount of value added at that stage.

* We will use another merge command to do this:

preserve

clear
set obs 5

gen prod_stage=_n

gen stage_fe = runiform()*2
  label var stage_fe "Production stage fixed effect"

save "Stage.dta", replace

restore

merge m:1 prod_stage using "Stage.dta"
drop _merge

* Let us first add on a production stage zero representing the initial "value" of the product idea.

* This is a little tricky to do.

expand 2 if prod_stage == 1, generate(expand_indicator)

* I will expand the data for production stage 1 and indicate the created data with expand indicator.
replace prod_stage = 0 if expand_indicator==1
drop expand_indicator

* Now I want to make sure that stage 0 production is not done by any contractors so there is only the company effect and product effect.

foreach v in cont_id cont_company_id cont_company_fe cont_fe rand_effect stage_fe {
  replace `v' = 0 if prod_stage == 0
}

*****
* Now let us start to generate the values of the products at each stage.

* This is a cumulative model ie. Value Added

* So in effect y=lambda*y[t-1] + XB + v

* To begin with we will generate the initial value of y.

gen value=abs(rnormal()) + prod_fe +  comp_fe if stage==0

* First let us double check to make sure our data is sorted properly.
sort prod_id prod_stage

* This is the retention value of a product from each previous stage.

* If lambda is low then it means that once the product is processed it cannot be used in the previous production stage.

* If lambda is high then it means that the previous value is retained plus any value added of progressive stages.
gen lambda=.95

* Now let us generate the cumulative value added data.
replace value=lambda*value[_n-1] + prod_fe +  comp_fe + cont_fe +cont_company_fe + stage_fe + rand_effect if stage>0

**** Simulation END

bysort prod_stage: sum value
* We can see that on average at each production stage there is an increasing value of the product.
* However, we can also see that the variance in values increases as the value increases.
* This is because of the cumulative variance effect of the various components combined with the high retention value of each previous stage (lambda).

* Now that we have data generated through a Value-added simulation we can start testing different value added estimators.

* That will be for a later post!