##############################################
# Lab 7: PEER INFLUENCE & QAP REGRESSION LAB #
##############################################
# NOTE: if you have trouble because some packages are not installed,
# see lab 1 for instructions on how to install all necessary packages.
###########################################################################
#
# Lab 7
#
# Part I - Develop peer influence models that meet the high standards of
# publication in todays top journals (i.e., addressing autocorrelation
# and issues of causality).
#
# Part II - Use QAP regression to predict increases in social interaction
# and task interaction for two classrooms. Basically run dyadic level
# regressions suitable for continuous variables and count data (what
# ERGM- Exponential Random Graph's cannot do).
#
###########################################################################
##########################################################################
# PART I -- PEER INFLUENCE
###########################################################################
#
# This lab examines peer effects in classroom data. It first introduces
# the reader to the necessary data processing and manipulation, then basic
# visualization related to peer-effects. It then introduces the Linear Network
# Autocorrelation Model (lnam), which can be used to better model the
# dependencies in network data (well, better than an OLS model anyway).
#
# The next section introduces the reader to the concept of matching--treating
# the explanatory variable of interest as an experimental "treatment"
# and then matching on other covariates. Matching allows the application
# of straightforward analytical techniques used to analyze experiments,
# namely comparison of means. Code for bootstrapping is introduced to allow for
# non-parametric estimation of the standard errors for the mean. Code is also provided
# to produce a dotplot showing the means and a (bootstrapped) 95% CI.
#
# The data was collected by Dan McFarland (McFarland 2001). The key measures of interest
# for this example are the time 1 and 2 measures of how much each student liked
# the subject, and the friend network at time 1.
#
# The codebook is available here: http://stanford.edu/~messing/codebook.v12.txt
###############################################################################
#Clear Variables
rm(list=ls(all=TRUE))
gc()
# Install and load necessary packages:
install.packages("reshape", dependencies = TRUE)
library(reshape)
library(igraph)
data(studentnets.peerinfl, package="NetData")
######################################################################################
# Now we need to reshape our attitudes data according to semester. The data is
# currently in "long" format, so that student ids and measures are repeated on
# multiple lines and "semester" is a single variable. We need
# to transform the data so that measures are repeated on a single
# line, or "wide" format. We'll use the excellent "reshape" package
# to accomplish this. Take a look at the documentation using ?reshape
?reshape
attitudesw = reshape( attitudes, idvar="std_id",
timevar="sem_id", direction="wide" )
# Notice that all semester 1 variables now feature a ".1" postscript, while all
# semester 2 variables feature a ".2" postscript
# We'll first visualize the data. We want to know whether the change in
# a student's appreciation of a subject (t_2 - t_1) changes in the direction of his
# friends' appreciation at t_1.
# So, we'll first take the difference of a few dependent variables
# for which we might expect peer influence (e.g. "sub" (liking of subject),
# "cmt" (liking of classmates), and "tch" (liking of teacher)),
# and then take the mean value for an individual's friends' values at t_1.
# First create the delta variables:
attitudesw$deltasub = attitudesw$sub.2 - attitudesw$sub.1
attitudesw$deltatch = attitudesw$tch.2 - attitudesw$tch.1
attitudesw$deltacmt = attitudesw$cmt.2 - attitudesw$cmt.1
# Next we'll create the mean of friends' sub.1.
# We'll use the adjacency matrix, multiply that by the attitudesw$sub.1
# (element-wise, not matrix multiplication), and then take the mean of each row.
sem1graph = graph.data.frame(d = sem1[1:2], directed=TRUE)
#sem1graph = network(x = sem1[1:2], directed=TRUE)
#sem2graph = graph.data.frame(d = sem2, directed=TRUE)
# But first, we'll need to clean the data to make sure the rows line up!
# let's check to see whether the edge data has the same cases as the attitudes data
which(!(V(sem1graph)$name %in% as.character(attitudesw$std_id)))
which(!(as.character(attitudesw$std_id) %in% V(sem1graph)$name ))
# They are not the same... we'll need to address this below.
# Now let's get the matrix representation of the network:
sem1matrix = get.adjacency(sem1graph, binary=T)
# This is often referred to as "W" in the literature
# When you have a large square matrix like this, you can get a good idea of
# it's density/sparsity using the image function.
image(sem1matrix)
# It's also generally good to know the degree distribution of any network you
# are handling:
plot(density(degree(sem1graph)))
# Looks like degree might be distributed exponentially. We can check as follows:
simexpdist100 = rexp(n=length(V(sem1graph)), rate = .100)
simexpdist125 = rexp(n=length(V(sem1graph)), rate = .125)
simexpdist150 = rexp(n=length(V(sem1graph)), rate = .150)
lines(density(simexpdist100), col="red")
lines(density(simexpdist125), col="blue")
lines(density(simexpdist150), col="green")
# It's not a precise fit, but it might be a good approximation.
# Let's reorder the attitudes data to make sure it's in the same order as
# sem1matrix
attitudesw = attitudesw[match(row.names(sem1matrix), attitudesw$std_id),]
# Make sure it worked (this should return 0 if so):
which(row.names(sem1matrix)!=as.character(attitudesw$std_id))
which(colnames(sem1matrix)!=as.character(attitudesw$std_id))
# Let's also make sure that igraph read the graph in correctly:
# (note that igraph can behave strangely if you pass it a
# parameter like "directed = FALSE" for a directed network like this one).
sem1$alter_id[sem1$std_id == 149824]
V(sem1graph)[0]
V(sem1graph)[unique(neighbors(sem1graph, v=0, mode = 1))]
# looks good.
# Now we can compute the mean of sub.1 for a student's friends,
# by element-wise multiplying sub.1 by the matrix and then taking the
# mean of each non-zero cell in each row:
# Let's first test this out with a simple matrix to make sure we know what
# we are doing:
(M = matrix(c(1,0,1,0,1,0,1,0,1), nrow=3, ncol=3))
(R = c(1,2,3))
R * M
# This is multiplying the first row by 1, the second row by 2, and so on.
# Instead we want:
t(R * t(M))
# This is multiplying the first column by 1, the second column by 2, and so on.
# Recall that we will be analyzing this data so that rows are cases, so
# this is what we want if we are going to calculate sub.1 for each
# student's friends' sub.1. (The first option would return the
# student's sub.1 for each non-zero row, not his/her friends' sub.1).
sem1Wxsub1 = t(attitudesw$sub.1 * t(sem1matrix))
sem1Wxtch1 = t(attitudesw$tch.1 * t(sem1matrix))
sem1Wxcmt1 = t(attitudesw$cmt.1 * t(sem1matrix))
# Now we'll take the mean of cells not equal to zero:
attitudesw$frsub1mean = numeric(length(attitudesw$std_id))
attitudesw$frtch1mean = numeric(length(attitudesw$std_id))
attitudesw$frcmt1mean = numeric(length(attitudesw$std_id))
for(i in 1:length(attitudesw$std_id)){
attitudesw$frsub1mean[i] = mean(sem1Wxsub1[i,sem1Wxsub1[i,]!=0 ], na.rm=T)
attitudesw$frtch1mean[i] = mean(sem1Wxtch1[i,sem1Wxtch1[i,]!=0], na.rm=T)
attitudesw$frcmt1mean[i] = mean(sem1Wxcmt1[i,sem1Wxcmt1[i,]!=0], na.rm=T)
}
# Now let's plot the data and look at the relationship:
plot(attitudesw$frsub1mean, jitter(attitudesw$deltasub))
plot(attitudesw$frtch1mean, jitter(attitudesw$deltatch))
plot(attitudesw$frcmt1mean, jitter(attitudesw$deltacmt))
# Looks noisy, but there might be a relationship for "sub".
############################################################################
# An alternative way to compute this is iteratively. The advantage
# to the following approach is that if you are working with a large
# data set, you do not need to load an adjacency matrix into
# memory, you can keep things in edgelist format, which is MUCH
# more efficient.
attitudesw$mfrsub.1 = numeric(length(attitudesw$sub.1))
attitudesw$mfrtch.1 = numeric(length(attitudesw$sub.1))
attitudesw$mfrcmt.1 = numeric(length(attitudesw$sub.1))
# If you thought the number of alters who like the class mattered
# more than the mean, you might uncomment the code for the nfrsubgt3
# variable and incorporate it into your analysis (two occurrences):
#attitudesw$nfrsubgt3 = numeric(length(attitudesw$sub.1))
for(i in 1:length(attitudesw$std_id)){
# first get alters of student i:
altrs = sem1$alter_id[ sem1$std_id == attitudesw$std_id[i] ]
# then get alters' attitudes
altatts = attitudesw$sub.1[ attitudesw$std_id %in% altrs ]
# now count how many friends like the class more than "3"
attitudesw$nfrsubgt3[i] = length(which(altatts > 3))
# then take the mean, ignoring NAs:
attitudesw$mfrsub.1[i] = mean(altatts, na.rm = TRUE)
# Note that this can all be done in one line of code:
# attitudesw$mfrsub.1[i] = mean(attitudesw$sub.1[ attitudesw$std_id %in% sem1$alter_id[sem1$std_id %in% attitudesw$std_id[i]]], na.rm=TRUE )
attitudesw$mfrtch.1[i] = mean(attitudesw$tch.1[ attitudesw$std_id %in% sem1$alter_id[sem1$std_id %in% attitudesw$std_id[i]]], na.rm=TRUE )
attitudesw$mfrcmt.1[i] = mean(attitudesw$cmt.1[ attitudesw$std_id %in% sem1$alter_id[sem1$std_id %in% attitudesw$std_id[i]]], na.rm=TRUE )
}
# if you are going to run this with a lot of data, you may wish to include
# the following lines inside your for-loop:
# gc()
# print(i)
# this will run the garbage collector, gc(), which helps R manage memory
# and print(i) will let you know your progress as R chugs away.
# The plot is exactly the same:
plot(attitudesw$mfrsub.1, jitter(attitudesw$deltasub))
# And the correlation should be 1
cor(attitudesw$mfrsub.1,attitudesw$frsub1mean, use= "complete")
# The plots suggest that there might be a relationship for "sub", or
# how much each student likes the subject.
# Let's cheat a little bit and run linear models predicting
# change in each of our variables based on mean friend's values at t=1.
# The data violates many of the assumptions of OLS regression so the
# estimates are not particularly meaningful, other than to let us know
# that we may be on the right path.
summary(lm(deltasub~mfrsub.1, data=attitudesw))
summary(lm(deltatch~mfrtch.1, data=attitudesw))
summary(lm(deltacmt~mfrcmt.1, data=attitudesw))
# Significance is always encouraging, even though we are only predicting a tiny
# fraction of the variance (R^2 = .026)
# We'll cheat a little more and remove everyone whose attitudes did not change:
summary(lm(deltasub~mfrsub.1, data=attitudesw, subset = deltasub!=0 ))
summary(lm(deltatch~mfrtch.1, data=attitudesw, subset = deltasub!=0 ))
summary(lm(deltacmt~mfrcmt.1, data=attitudesw, subset = deltasub!=0 ))
# Very nice, for "sub" the effect grows in strength, signficance,
# and R^2 goes up to .067.
# Let's also take a look at the plot:
pdf("7.1_Peer_effect_on_opinion_of_subject.pdf")
plot( x = attitudesw$mfrsub.1,
y = jitter(attitudesw$deltasub),
main = "Peer effects on opinion of subject",
xlab = "mean friends' opinion",
ylab = expression(Delta~"opinion from "~t[1]*~to~t[2])
)
# We can draw a regression line on the plot based on the regression
# we estimated above:
abline(lm(deltasub~mfrsub.1, data=attitudesw ))
dev.off()
# We've got some evidence of a peer effect on how much a student
# likes the subject. Our evidence is based on temporal precendence:
# your friends opinion of the subject predicts how much
# your opinion the subject changes from t_1 to t_2.
############################################################################
# Now let's run some "real" models.
# We'll first estimate the model Y_2 = Y_1 + W*Y_alter, using a network
# autocorrelation model. Read up on it here:
# Roger Th. A. J. Leenders, Modeling social influence through network
# autocorrelation: constructing the weight matrix, Social Networks, Volume 24, Issue 1, January 2002, Pages 21-47, ISSN 0378-8733, DOI: 10.1016/S0378-8733(01)00049-1.
# (http://www.sciencedirect.com/science/article/B6VD1-44SKCC2-1/2/cae2d6b4cf4c1c21f4f870fd2d58b5cc)
# This model arguably takes into accounts the correlation between friends in
# the disturbances (error term). Yet there are problems with this kind of model.
# See:
# Eric J. Neuman, Mark S. Mizruchi, Structure and bias in the network autocorrelation
# model, Social Networks, In Press, Corrected Proof, Available online 31 May 2010, ISSN 0378-8733, DOI: 10.1016/j.socnet.2010.04.003.
# (http://www.sciencedirect.com/science/article/B6VD1-506H0SN-1/2/73a16c627271d29d9283b6d69b873a07)
# With that in mind, let's proceed:
library(sna)
library(numDeriv)
# We'll use the mfrsub.1 variable to approximate the W*Y_alter term. (If
# we actually use the matrix, we cannot estimate the model--there
# would be the same number of variables as cases).
# Let's remove the NA values from the attitudes measures
# (otherwise lnam will crash). We'll remove rows with NAs for
# sub.1, sub.2, or mfrsub.1. Note that we could remove any row with an
# NA value with the syntax: na.omit(attitudesw), but if there are
# rows that contain our variables of interest but are NA for some other
# variable, na.omit would also drop these rows and we would lose valuable
# data.
atts = attitudesw[!is.na(attitudesw$sub.2),]
atts = atts[!is.na(atts$sub.1),]
atts = atts[!is.na(atts$mfrsub.1),]
W = sem1matrix
# Now we'll make sure the rows and columns in W are in the same
# order as atts:
W = W[match(atts$std_id,row.names(W)), match(atts$std_id,colnames(W))]
# Let's make sure (this will return 0 if we did it right):
which(rownames(W) != atts$std_id)
# Now we'll estimate the peer influence model using lnam, which stands for
# Linear Network Autocorrelation Model.
pim1<-lnam(atts$sub.2,cbind(atts$sub.1, atts$mfrsub.1), W2 = W)
summary(pim1)
# This shows a pretty strong peer effect on subjects at time2. Notice that
# the adjusted R^2 is .41. This is not good news. Even after controlling
# for t1 values of attitude for the subject, we are predicting less than
# half of the variance in t_2 values. That means that there is a good
# chance of omitted variable bias.
# Below is an attempt to run the model with the actual WYalt as the X matrix.
# For the reasons stated above, it doesn't work. It is included in case the
# reader wants to uncomment this code and see for him/herself.
#WYalt = sem1Wxsub1
#WYalt = WYalt[match(atts$std_id,row.names(WYalt)), match(atts$std_id,colnames(WYalt))]
#image(WYalt)
#
##pim2<-lnam(y = atts$sub.2, x = WYalt, W2 = W)
#pim2<-lnam(y = atts$sub.2, x = WYalt)
#summary(pim2)
############################################################################
# Matching
#
# One method that is less vulnerable to some problems related to regression is
# matching. In matching, one divides the key explanatory variable into discrete
# "treatments" and then matches participants so that they are as similar on the
# remaining set of covariates as possible. A simple difference in means/medians/
# quantiles can be used. From there we can use a simple t-test to estimate the
# "treatment" effect without having to worry about multicolinearity.
# If the data violates the normality assumption (as will be the
# case here), non-parametric methods such as a permutation test or bootstrapped
# standard errors can be used. This approach also also violates less statistical
# assumptions and therefore should be expected to be less biased than running a
# simple OLS model.
# For the original article on matching, see:
# Donald B. Rubin, Matching to Remove Bias in Observational Studies
# Biometrics, Vol. 29, No. 1 (Mar., 1973), pp. 159-183
# Stable URL: http://www.jstor.org/stable/2529684
# Also worthy of note is the Rubin Causal Model, which uses matching.
# Read up on it here: http://en.wikipedia.org/wiki/Rubin_Causal_Model
# However, matching is not a cure-all and is obviously sensitive to the variables
# on which one matches. For an excellent article on this, see:
# Jeffrey A. Smith, Petra E. Todd, Does matching overcome LaLonde's critique of
# nonexperimental estimators?, Journal of Econometrics, Volume 125, Issues 1-2,
# Experimental and non-experimental evaluation of economic policy and models,
# March-April 2005, Pages 305-353.
# (http://www.sciencedirect.com/science/article/B6VC0-4CXTY59-1/2/06728bd79899fd9a5b814dfea9fd1560)
# A good introduction to matching using the MatchIt package can be found courtesy
# of Gary King here: http://gking.harvard.edu/matchit/docs/matchit.pdf
library(MatchIt)
# First, we dichotomize the independent variable of interest, mean friend opinion regarding
# subject: mfrsub.1
atts$mftreat = ifelse(atts$mfrsub.1 > 3, 1,0)
atts = atts[c("mftreat", "deltasub", "mfrsub.1","egrd.1", "sub.1", "tot.1", "frn.1", "cmt.1", "prev.1")]
# Before we do any matching, let's take a look at our new treatment variable.
# We'll make this into a factor, which will help interpret our results. A factor
# is a special type of object in R that is used for qualitative (nominal or ordinal)
# data/variables. A factor is also generally more efficient with respect to memory,
# because it stores an integer value corresponding to a label (e.g. 1 for "low" and 2
# for "high") rather than storing the entire string.
# Factors have an order, which you can set by ordering the labels in ascending order
# when you create the factor, or by using the reorder() function. We'll talk about this
# more in the section on dotplots below.
atts$mftreatf = factor(atts$mftreat, labels = c("low", "high"))
# Let's compute the means:
tapply(X = atts$deltasub, INDEX = list(atts$mftreatf), FUN = mean)
# Note, the tapply function is useful if you have more than one factor for which
# you want to calcuate some function over. For example, if there were another
# factor, say exptreat, we could get the mean for each cell using this syntax:
# tapply(X = atts$deltasub, INDEX = list(atts$mftreat, atts$exptreat), FUN = mean)
# A t-test is not totally kosher here because the DV is ordinal, and so it violates
# assumptions behind the t-test.
# But let's take a look anyway:
t.test(deltasub~mftreatf, data=atts)
# A better alternative is the Wilcoxon test (a.k.a. Mann-Whitney test),
# suitable for ordinal variables.
wilcox.test(x=atts$deltasub[atts$mftreatf=="low"],
y=atts$deltasub[atts$mftreatf=="high"],
conf.int = TRUE)
# Now let's match and redo the analyses.
# We normally would assign each subject/case to treatment or control groups based on
# whatever corresponds more closely to a treatment or a control. In this case, it
# is debatable whether having friends who like the class is more of a "treatment"
# than having friends who do not like the class. However, there are 246 students
# whose friends do not like the class versus 99 whose friends do like the class,
# so we can say that having friends who like the class is certainly more unusual
# than the opposite. Accordingly, we'll conceptualize having friends who like the
# class as the treatment.
# In the actual matching proceedure, what we'll do is try to find individuals in
# the "control" condition (having friends who do not like the subject) who look
# most similar to each treated individual based on a variety of variables. Ideally,
# we will achieve "balance" on the covariates between those in the treatment condition,
# and those we select from the control condition, so that the individuals are as similar
# as possible.
# Here are the variables we'll use (all are at time 1):
# sub.1 - how much the student likes the subject
# egrd.1 - expected grade
# tot.1 - how much the student likes the way the course is taught
# frn.1 - how interesting the subject finds his/her friends
# cmt.1 - how interesting the student finds his/her classmates
# prev.1 - has the student had the teacher previously)
# We'll use "nearest neighbor" matching, which uses a distance measure to
# select the individual in the "control" group that is closest to each
# individual in the "treatment" group. The default distance measure is the
# logistic regression propensity score.
# The codebook for this data is available here: http://stanford.edu/~messing/codebook.v12.txt
m.out = matchit(mftreat ~ egrd.1 + sub.1 + tot.1 + frn.1 + cmt.1 + prev.1,
data = atts, method = "nearest")
summary(m.out)
plot(m.out)
# Now let's assess the extent to which our matching algorithm effectively matched each
# treatment subject to a control subject, e.g., the extent to which we achieved balance.
# The plot function displays Q-Q plots of the control (X-axis) versus treatment groups
# (Y-axis) for the original sample and the matched sample. The Q-Q plots indicate perfect
# balance on the covariate in question if all dots are perfectly aligned on the
# 45 degree line. The further the deviation from this line, the worse the samples are
# balanced on this variable.
# Based on the plots, matching appears to have improved balance on our covariates a
# little bit, but not perfectly.
# Ideally we want to experiment with other covariates and try to find the combination
# that best and most parsimoniously captures the key phenomena we want to control for.
matchatts = match.data(m.out)
# We'll make our treatment variable into a factor:
matchatts$mftreatf = factor(matchatts$mftreat, labels = c("low", "high"))
# Let's compute the means:
tapply(X = matchatts$deltasub, INDEX = list(matchatts$mftreatf), FUN = mean)
t.test(deltasub~mftreatf, data=matchatts)
wilcox.test(x=matchatts$deltasub[matchatts$mftreatf=="low"],
y=matchatts$deltasub[matchatts$mftreatf=="high"],
conf.int = TRUE)
# We can also perform a full permuation test. A permutation test is an exact test of
# whether it is possible to reject the null hypothesis that two distributions are
# the same. It works by first calulating the mean difference, which fuctions as the
# test statistic. The difference in sample means is calculated and recorded for every
# possible way of dividing these pooled values into two groups of size n_A and n_B
# for every permutation of the group labels A and B. The set of these calculated
# differences is the exact distribution of possible differences under the null
# hypothesis that group label does not matter. This test fuctions like a t.test but
# does not rely on any distributional assumptions about the data.
# For additional information, see:
# http://en.wikipedia.org/wiki/Resampling_%28statistics%29#Permutation_tests
library(coin)
independence_test(deltasub ~ mftreatf, data = matchatts,
distribution = exact())
# Another alternative is to look at bootstrapped estimates of the mean and standard
# error. This is often an attractive way to estimate statistical significance
# because our data violates the distributional assumptions involved in estimating
# standard errors in this way (namely, that our dependent variable resembles anything
# like a normal or t distribution--it cannot because it is ordinal with only
# 7 levels). Because boostrapping relies on resampling instead of distributional
# assumptions to estimate variance, it is more robust.
# There are packages available to generate bootstrapped estimates, but seeing the code
# is a valuable way to get a sense of exactly how bootstrapping works. Here is the code
# for a bootstrapping function designed to estimate the mean and standard error:
b.mean <- function(data, nsim) {
# Create a list object to store the sets of resampled data (in this case nsim = 1000).
resamples = list()
# Create a vector of the same length to store the resampled means
r.mean = numeric(nsim)
# Generate a sample with replacement
for(i in 1:nsim){
# generates a random sample of our data with replacement
resamples[[i]] = sample(x=data, size=length(data), replace=T)
# Now calcuate the mean for this iteration
r.mean[i] = mean(resamples[[i]], na.rm=T)
}
# Calculate the mean of the mean of the simulated estimates above:
mean.adj = mean(r.mean, na.rm=T)
# Calculate how this differs from the arithmatic mean estimate of the original
# sample:
bias = mean(data, na.rm=T) - mean.adj
# Generate the standard error of the estimate
std.err <- sqrt(var(r.mean))
# Return results
return( data.frame(mean = mean(data, na.rm=T), #the mean estimate
mean.adj = mean.adj, # the adjusted estimate based on simulations
bias = bias, # the mean minus the adjusted estimate
std.err=std.err # the standard error of the estimate (based on simulations)
)
)
}
# Before we use it on our data, let's make sure it works. We will simulate some
# data for which we know the parameters and then estimate the parameters using
# the bootstrap:
simdata = rnorm(n = 1000, mean = 5, sd = 1)
b.mean(simdata, nsim=1000)
# Calculate the theoretical mean and standard error
mean(simdata)
(se = sd(simdata)/sqrt(length(simdata)-1))
# here we use the formula SE = SD/sqrt(N - 1)
# We have demonstrated that our bootstraping function is a good estimator
# for "perfect" data that meets the assumptions behind traditional
# statistical tests. Now let's move on to our data.
# For those with friends with on average postive attitudes:
b.mean(data=matchatts$deltasub[matchatts$mftreatf=="high"],
nsim=1000)
# For those with friends with on average negative attitudes:
b.mean(data=matchatts$deltasub[matchatts$mftreatf=="low"],
nsim=1000)
# Here's how to do it using the boot package, which is faster and more
# extensible than the R-code above:
library(boot)
# You have to write a special function for the boot package, which
# can be confusing at first, but is useful for more complicated
# types of bootstrapping. The boot package is also orders of magnitude
# faster than doing things in R, which is useful if you start to do
# more complicated bootstrapping or use a higher number of simulations.
# First, write a basic function that will return the estimate in question
# for x data using d cases:
samplemean <- function(x, d) {
return(mean(x[d]))
}
(bootlow = boot(data = matchatts$deltasub[matchatts$mftreatf=="low"],
statistic = samplemean,
R = 1000))
(boothigh =boot(data = matchatts$deltasub[matchatts$mftreatf=="high"],
statistic = samplemean,
R = 1000))
# Note that estimates of the standard errors will be slightly different
# each time due to randomness involved in estimation. As you might expect,
# as the number of simulations increases, variance will decrease.
# Based on our bootstrapped estimates and standard error calculations,
# it's clear that the two means are quite different. So, based on our
# matched samples, it seems there is a significant peer effect.
##########################################################################
# Let's try full matching so that we are not throwing away any of our data:
m.out = matchit(mftreat ~ egrd.1 + sub.1 + tot.1 + frn.1 + cmt.1 + prev.1,
data = atts, method = "full")
matchatts = match.data(m.out)
# We'll make this into a factor, which will help inturpret our results:
matchatts$mftreatf = factor(matchatts$mftreat, labels = c("low", "high"))
# Note that we'll now have to use the weights that the output provided because
# we told it to use the full sample.
# Let's compute the means:
tapply(X = matchatts$deltasub,
INDEX = list(matchatts$mftreatf),
FUN = weighted.mean, weights = matchatts$weights )
# There is no straightforward way to compute a t-test in R with weighted data.
# We can however compute weighted standard errors and a corresponding
# 95% CI:
# Recall that a 95% CI of an estimate is calculated by taking the estimate
# +/- the standard error * 1.96. Actually, +/- 1.96 is just an estimate of the
# the 0.05 and 0.95 quantiles of a statistical distribution. We'll use
# qt(quantile, N-1) in this case.
library(Hmisc)
meanlow = wtd.mean(matchatts$deltasub[matchatts$mftreatf=="low"],
weights = matchatts$weights[matchatts$mftreatf=="low"])
varlow = wtd.var(matchatts$deltasub[matchatts$mftreatf=="low"],
weights = matchatts$weights[matchatts$mftreatf=="low"])
selow = sqrt(varlow)/sqrt(sum(matchatts$weights[matchatts$mftreatf=="low"]))
meanlow + selow * qt(.025, ( sum(matchatts$weights[matchatts$mftreatf=="low"]) - 1))
meanlow + selow * qt(.975, ( sum(matchatts$weights[matchatts$mftreatf=="low"]) - 1))
meanhigh = wtd.mean(matchatts$deltasub[matchatts$mftreatf=="high"],
weights = matchatts$weights[matchatts$mftreatf=="high"])
varhigh = wtd.var(matchatts$deltasub[matchatts$mftreatf=="high"],
weights = matchatts$weights[matchatts$mftreatf=="high"])
sehigh = sqrt(varhigh)/sqrt(sum(matchatts$weights[matchatts$mftreatf=="high"]))
meanhigh + sehigh * qt(.025, ( sum(matchatts$weights[matchatts$mftreatf=="high"]) - 1))
meanhigh + sehigh * qt(.975, ( sum(matchatts$weights[matchatts$mftreatf=="high"]) - 1))
# Here's how to do it with the boot package:
sample.wtd.mean <- function(x, d, wts) {
return(weighted.mean(x[d], w = wts))
}
(bootlow = boot(data = matchatts$deltasub[matchatts$mftreatf=="low"],
wts = matchatts$weights[matchatts$mftreatf=="low"],
statistic = sample.wtd.mean,
R = 1000))
(boothigh =boot(data = matchatts$deltasub[matchatts$mftreatf=="high"],
wts = matchatts$weights[matchatts$mftreatf=="high"],
statistic = sample.wtd.mean,
R = 1000))
# We can make a nice plot of the results with the lattice package:
library(lattice)
# Dot plots:
visual = list()
visual$var = c("Friends do not like subj.", "Friends like subj.")
visual$M = c(bootlow$t0, boothigh$t0)
visual$se = c(sd(bootlow$t), sd(boothigh$t))
visual$N = c(length(matchatts$weights[matchatts$mftreatf=="high"]),
length(matchatts$weights[matchatts$mftreatf=="low"]))
visual$lo = visual$M - visual$se
visual$up = visual$M + visual$se
visual = as.data.frame(visual)
#print it out:
pdf(file="7.2_dotplot_matchoutput_bootstrapped_SE.pdf", width = 6, height = 3)
dotplot( reorder(var,M) ~ M,
data = visual,
main="Effect of friends' attitude on opinion change,\nmatched samples with bootstrapped SE",
panel = function (x, y, subscripts, groups, ...) {
panel.dotplot(x, y)
panel.segments(x0=visual$lo[subscripts],y0=as.numeric(y),
x1=visual$up[subscripts],y1=as.numeric(y),
lty = 1 )
},
xlab=expression(Delta~"opinion from "~t[1]*~to~t[2]),
xlim= c(-.5, .2)
)
dev.off()
# Very nice!
############################################
# QUESTION #1 - What do these results show?
# What assumptions are being made?
# What might lead you to question the results? What would improve them?
############################################
#########################################################################
# Extra-credit:
#
# Use the McFarland dataset to acquire other individual level variables
# and develop better matching. Then follow this lab and explore the
# variety of outcomes suitable for illustrating peer effects in
# classrooms (e.g., PSAT scores, engagement, conflict, etc). Results
# are likely suitable for an A-journal publication.
#
# Again, the codebook is available here: http://stanford.edu/~messing/codebook.v12.txt
#
#########################################################################
#########################################################################
# PART II -- QAP REGRESSION
#########################################################################
#
# Here we want to compare different classrooms and discern what leads
# participants to become more or less sociable and academically engaged
# with one another.
#
# Class m173 is an accelerated math class of all 10th graders taught
# by a rigid teacher who used teacher-centered instructional formats.
# Class m 182 is a regular math class of 10th and 11th grade students
# taught by a casual teacher who used student-centered instructional formats.
#
#########################################################################
# igraph and sna don't mix well. Before we load sna, we will detach igraph.
detach(package:igraph)
# Load the "sna" library
library(sna)
#Clear Variables
rm(list=ls(all=TRUE))
gc()
#(2) QAP Regressions for m l73
# Note that each matrix must be the same size (n x n). You'll want to make
# sure all input and output matrices are the same size. Predictor matrices
# thus are not individual level variables; they are differences, matches,
# etc. You'd use a matrix of differences in of quantitative (node level)
# variables, and matches in the case of categorical/dichotomous (node
# level) variables. There's a really easy way to do it in R
# (say your node-level variables are x and y)
# x <- seq(1,4)
# y <- seq(2,5)
# outer(x,y,FUN="-") # find the difference between x and y
# outer(x,y,FUN="==") # produce 0/1 similairty matrix
# For this lab we will use matrices that were previously generated in UCINet
# Loading predictor matrices. Each matrix is a n x n matrix
# data is saved in a CSV format. You can open the files in Excel
# to see the structure.
data(studentnets.mrqap173, package="NetData")
# Look at what we loaded via
ls()
# We need the data in matrix format
# predictor matrices
m173_sem1_SSL <- as.matrix(m173_sem1_SSL)
m173_sem1_TSL <- as.matrix(m173_sem1_TSL)
m173_sem1_FRN <- as.matrix(m173_sem1_FRN)
m173_sem1_SEAT <- as.matrix(m173_sem1_SEAT)
m173_sem1_RCE <- as.matrix(m173_sem1_RCE)
m173_sem1_GND <- as.matrix(m173_sem1_GND)
# Load response matrices
m173_sem2_SSL <- as.matrix(m173_sem2_SSL)
m173_sem2_TSL <- as.matrix(m173_sem2_TSL)
# In order to run the QAP regression we must create an array of matrices
# containing all the predcitor matrices. We are, in effect, creating a
# 3-d matrix (predcitor x n x n).
# Important: All matrices must be the same size!
response_matrices <- array(NA, c(6, length(m173_sem1_SSL[1,]),length(m173_sem1_SSL[1,])))
response_matrices[1,,] <- m173_sem1_SSL
response_matrices[2,,] <- m173_sem1_TSL
response_matrices[3,,] <- m173_sem1_FRN
response_matrices[4,,] <- m173_sem1_SEAT
response_matrices[5,,] <- m173_sem1_RCE
response_matrices[6,,] <- m173_sem1_GND
##############################
#(2a) SSL2 <- SSL1 + TSL1 + FRN1 + SEAT1 + RCE + GND
##############################
# Fit a netlm model by using netlm, the response matrix and the array of predictor matrices
# This may take a LONG time.
nl<-netlm(m173_sem2_SSL,response_matrices)
# Make the model easier to read by adding lables for each predictor matrix.
nlLabeled <- list()
nlLabeled <- summary(nl)
# Labels are provided in the same order as they were assigned in the response_matrices array
nlLabeled$names <- c("Intercept", "Social normalized and labeled (SSL1)", "Task normalized and labeled (TSL1)", "Friends 1=friend, 2=best friend(FRN1)", "Seat in first semester (Seat1)","Race (RCE)","Gender (GND)")
# Round the ocefficients to two decimals
nlLabeled$coefficients = round(nlLabeled$coefficients, 2)
nlLabeled
##############################
#(2b) TSL2 <- TSL1 + SSL1 + FRN1 + SEAT1 + RCE + GND
##############################
# Fit a netlm model by using netlm, the response matrix and the array of predictor matrices
nl<-netlm(m173_sem2_TSL,response_matrices)
#make the model easier to read
nlLabeled <- list()
nlLabeled <- summary(nl)
nlLabeled$names <- c("Intercept", "Social normalized and labeled (SSL1)", "Task normalized and labeled (TSL1)", "Friends 1=friend, 2=best friend(FRN1)", "Seat in first semester (Seat1)","Race (RCE)","Gender (GND)")
nlLabeled$coefficients = round(nlLabeled$coefficients, 2)
nlLabeled
##############################
#(3) QAP Regressions for m 182
##############################
# Repeat for class m 182
# Clear Variables
rm(list=ls(all=TRUE))
gc()
data(studentnets.mrqap182, package = "NetData")
# Look at what we loaded via
ls()
# Again, we need the data in matrix format
# predictor matrices
m182_sem1_SSL <- as.matrix(m182_sem1_SSL)
m182_sem1_TSL <- as.matrix(m182_sem1_TSL)
m182_sem1_FRN <- as.matrix(m182_sem1_FRN)
m182_sem1_SEAT <- as.matrix(m182_sem1_SEAT)
m182_sem1_RCE <- as.matrix(m182_sem1_RCE)
m182_sem1_GND <- as.matrix(m182_sem1_GND)
#response matrices
m182_sem2_SSL <- as.matrix(m182_sem2_SSL)
m182_sem2_TSL <- as.matrix(m182_sem2_TSL)
#This class will require you to make multiple extraction operations.
#A student exits the class at the end of first semester and another
#enters at the start of second semester. These students must be removed
#before you can conduct QAP regression (you need the same row-column orderings).
#Extract actor 15 for TSL and SSL of second semester (# 15 <- new student).
response_matrices <- array(NA, c(6, length(m182_sem1_SSL[1,]),length(m182_sem1_SSL[1,])))
response_matrices[1,,] <- m182_sem1_SSL
response_matrices[2,,] <- m182_sem1_TSL
response_matrices[3,,] <- m182_sem1_FRN
response_matrices[4,,] <- m182_sem1_SEAT
response_matrices[5,,] <- m182_sem1_RCE
response_matrices[6,,] <- m182_sem1_GND
##############################
#(3a) SSL2 <- SSL1 + TSL1 + FRN1 + SEAT1 + RCE + GND
##############################
#Fit a netlm model
nl<-netlm(m182_sem2_SSL,response_matrices)
#make the model easier to read
nlLabeled <- list()
nlLabeled <- summary(nl)
nlLabeled$names <- c("Intercept", "Social normalized and labeled (SSL1)", "Task normalized and labeled (TSL1)", "Friends 1=friend, 2=best friend(FRN1)", "Seat in first semester (Seat1)","Race (RCE)","Gender (GND)")
nlLabeled$coefficients = round(nlLabeled$coefficients, 2)
nlLabeled
##############################
#(3b) TSL2 <- SSL1 + TSL1 + FRN1 + SEAT1 + RCE + GND
##############################
#Fit a netlm model
nl<-netlm(m182_sem2_TSL,response_matrices)
#make the model easier to read
nlLabeled <- list()
nlLabeled <- summary(nl)
nlLabeled$names <- c("Intercept", "Social normalized and labeled (SSL1)", "Task normalized and labeled (TSL1)", "Friends 1=friend, 2=best friend(FRN1)", "Seat in first semester (Seat1)","Race (RCE)","Gender (GND)")
nlLabeled$coefficients = round(nlLabeled$coefficients, 2)
nlLabeled
#Report your results. Describe what they mean? Repeat for ETSL of second semester.
#######################################################################
# QUESTION #2 - Compare your results for m 173 and m 182.
# Do the classes have different results?
# If not, why not?
# If so, why so?
# What sort of story can you derive about the change in
# task and social interactions within these classrooms?
#
# Remember you are predicting changes in relations of sociable or task
# interaction using other types of relations (i.e., friendship,
# proximate seating, same race, etc).
#########################################################################
#
# Extra-credit - run the QAP regression part of the lab on the
# entire McFarland dataset of classrooms, and perform meta-analyses
# on the results using multi-level modeling - then presto -
# you will have results on academic and social relations in classrooms
# which is suitable for publication in an A-journal.
#
#########################################################################