Linear Network Autocorrelation Model
WEN, Tzai-Hung (NTU Geography)
Source: McFarland, Daniel, et.al. 2010. “Social Network Analysis Labs in R.” Stanford University.
library(reshape)
library(igraph)
library(sna)
library(numDeriv)
1. Load Students Data from NetData
data(studentnets.peerinfl, package="NetData")
2. View each student’s attributes
Student’s Preference: e.g. “sub” (liking of subject), “cmt” (liking of classmates), and “tch” (liking of teacher)
head(attitudes, n=10)
## std_id sem_id cls_id tlks tlkt like_c imp egrd jgrd sub tot frn cmt tch
## 1 16144 1 851 4 2 NA NA 2 -1 2 2 3 2 2
## 2 16181 1 851 4 2 2 NA 3 0 2 3 4 3 2
## 3 16247 1 851 4 4 2 NA 4 0 2 2 4 3 2
## 4 16399 1 851 4 2 3 NA 4 1 3 2 3 1 2
## 5 15469 1 901 2 2 3 NA 2 0 3 4 4 4 4
## 6 15947 1 901 3 2 3 NA 3 0 3 4 4 3 4
## 7 15973 1 901 3 2 3 NA 4 0 3 4 4 3 4
## 8 16073 1 901 4 3 3 NA 3 0 3 3 4 3 4
## 9 16177 1 901 2 2 3 NA 3 0 3 3 3 2 4
## 10 16034 1 901 3 3 4 NA 3 0 3 4 4 3 4
## voln call misb sanc sftch tfstd tfros chal prev
## 1 NA NA NA NA NA NA NA NA 0
## 2 NA NA NA NA NA NA NA 3 0
## 3 NA NA NA NA NA NA NA 2 0
## 4 NA NA NA NA NA NA NA 2 0
## 5 NA NA NA NA NA NA NA 1 0
## 6 NA NA NA NA NA NA NA 3 0
## 7 NA NA NA NA NA NA NA 3 0
## 8 NA NA NA NA NA NA NA 3 0
## 9 NA NA NA NA NA NA NA 3 0
## 10 NA NA NA NA NA NA NA 3 1
3. View each student’s friend relationships in semester 1 (edge-list)
head(sem1, n=10)
## std_id alter_id timea ona oca acta nacta wrka nbra otha wkna bfa lva
## 2 149824 119516 5 1 0 0 0 0 0 0 0 0 0
## 3 149824 122634 1 1 0 0 0 0 0 0 0 0 0
## 4 149824 114679 1 0 1 0 0 0 0 0 0 0 0
## 7 16868 17142 0 1 0 1 0 0 0 0 0 0 0
## 8 16877 16681 1 0 1 1 0 0 0 0 0 0 0
## 9 16877 16769 1 0 1 1 0 0 0 0 0 0 0
## 10 16877 16851 3 1 0 0 0 0 1 0 0 0 0
## 11 16877 16868 7 0 1 0 0 0 1 0 0 0 0
## 13 16902 16690 0 0 1 0 0 0 0 0 0 0 0
## 14 16902 16868 0 0 1 1 0 0 0 0 0 0 0
## cls_id
## 2 251
## 3 251
## 4 251
## 7 851
## 8 851
## 9 851
## 10 851
## 11 851
## 13 851
## 14 851
4. Reshape
attitudesw = reshape( attitudes, idvar="std_id",timevar="sem_id", direction="wide" )
head(attitudesw, n=10)
## std_id cls_id.1 tlks.1 tlkt.1 like_c.1 imp.1 egrd.1 jgrd.1 sub.1 tot.1
## 1 16144 851 4 2 NA NA 2 -1 2 2
## 2 16181 851 4 2 2 NA 3 0 2 3
## 3 16247 851 4 4 2 NA 4 0 2 2
## 4 16399 851 4 2 3 NA 4 1 3 2
## 5 15469 901 2 2 3 NA 2 0 3 4
## 6 15947 901 3 2 3 NA 3 0 3 4
## 7 15973 901 3 2 3 NA 4 0 3 4
## 8 16073 901 4 3 3 NA 3 0 3 3
## 9 16177 901 2 2 3 NA 3 0 3 3
## 10 16034 901 3 3 4 NA 3 0 3 4
## frn.1 cmt.1 tch.1 voln.1 call.1 misb.1 sanc.1 sftch.1 tfstd.1 tfros.1
## 1 3 2 2 NA NA NA NA NA NA NA
## 2 4 3 2 NA NA NA NA NA NA NA
## 3 4 3 2 NA NA NA NA NA NA NA
## 4 3 1 2 NA NA NA NA NA NA NA
## 5 4 4 4 NA NA NA NA NA NA NA
## 6 4 3 4 NA NA NA NA NA NA NA
## 7 4 3 4 NA NA NA NA NA NA NA
## 8 4 3 4 NA NA NA NA NA NA NA
## 9 3 2 4 NA NA NA NA NA NA NA
## 10 4 3 4 NA NA NA NA NA NA NA
## chal.1 prev.1 cls_id.2 tlks.2 tlkt.2 like_c.2 imp.2 egrd.2 jgrd.2 sub.2
## 1 NA 0 851 4 2 3 3 1 -1 2
## 2 3 0 851 3 2 2 2 3 1 1
## 3 2 0 726 4 1 3 3 4 1 4
## 4 2 0 851 4 2 2 3 2 -1 3
## 5 1 0 NA NA NA NA NA NA NA NA
## 6 3 0 NA NA NA NA NA NA NA NA
## 7 3 0 NA NA NA NA NA NA NA NA
## 8 3 0 692 2 3 2 3 2 0 2
## 9 3 0 692 3 2 2 3 3 0 2
## 10 3 1 NA NA NA NA NA NA NA NA
## tot.2 frn.2 cmt.2 tch.2 voln.2 call.2 misb.2 sanc.2 sftch.2 tfstd.2
## 1 1 2 1 2 1 1 1 1 0 0
## 2 3 3 2 2 3 2 1 1 0 0
## 3 2 4 3 3 2 2 3 1 2 2
## 4 1 2 2 1 1 1 1 1 0 -1
## 5 NA NA NA NA NA NA NA NA NA NA
## 6 NA NA NA NA NA NA NA NA NA NA
## 7 NA NA NA NA NA NA NA NA NA NA
## 8 2 4 3 2 2 3 2 2 0 1
## 9 3 4 3 2 1 2 1 1 0 1
## 10 NA NA NA NA NA NA NA NA NA NA
## tfros.2 chal.2 prev.2
## 1 1 NA NA
## 2 0 NA NA
## 3 2 NA NA
## 4 0 NA NA
## 5 NA NA NA
## 6 NA NA NA
## 7 NA NA NA
## 8 1 NA NA
## 9 1 NA NA
## 10 NA NA NA
5. Calculate Peer Influences
attitudesw$mfrsub.1 = 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)
}
head(attitudesw, n=10)
## std_id cls_id.1 tlks.1 tlkt.1 like_c.1 imp.1 egrd.1 jgrd.1 sub.1 tot.1
## 1 16144 851 4 2 NA NA 2 -1 2 2
## 2 16181 851 4 2 2 NA 3 0 2 3
## 3 16247 851 4 4 2 NA 4 0 2 2
## 4 16399 851 4 2 3 NA 4 1 3 2
## 5 15469 901 2 2 3 NA 2 0 3 4
## 6 15947 901 3 2 3 NA 3 0 3 4
## 7 15973 901 3 2 3 NA 4 0 3 4
## 8 16073 901 4 3 3 NA 3 0 3 3
## 9 16177 901 2 2 3 NA 3 0 3 3
## 10 16034 901 3 3 4 NA 3 0 3 4
## frn.1 cmt.1 tch.1 voln.1 call.1 misb.1 sanc.1 sftch.1 tfstd.1 tfros.1
## 1 3 2 2 NA NA NA NA NA NA NA
## 2 4 3 2 NA NA NA NA NA NA NA
## 3 4 3 2 NA NA NA NA NA NA NA
## 4 3 1 2 NA NA NA NA NA NA NA
## 5 4 4 4 NA NA NA NA NA NA NA
## 6 4 3 4 NA NA NA NA NA NA NA
## 7 4 3 4 NA NA NA NA NA NA NA
## 8 4 3 4 NA NA NA NA NA NA NA
## 9 3 2 4 NA NA NA NA NA NA NA
## 10 4 3 4 NA NA NA NA NA NA NA
## chal.1 prev.1 cls_id.2 tlks.2 tlkt.2 like_c.2 imp.2 egrd.2 jgrd.2 sub.2
## 1 NA 0 851 4 2 3 3 1 -1 2
## 2 3 0 851 3 2 2 2 3 1 1
## 3 2 0 726 4 1 3 3 4 1 4
## 4 2 0 851 4 2 2 3 2 -1 3
## 5 1 0 NA NA NA NA NA NA NA NA
## 6 3 0 NA NA NA NA NA NA NA NA
## 7 3 0 NA NA NA NA NA NA NA NA
## 8 3 0 692 2 3 2 3 2 0 2
## 9 3 0 692 3 2 2 3 3 0 2
## 10 3 1 NA NA NA NA NA NA NA NA
## tot.2 frn.2 cmt.2 tch.2 voln.2 call.2 misb.2 sanc.2 sftch.2 tfstd.2
## 1 1 2 1 2 1 1 1 1 0 0
## 2 3 3 2 2 3 2 1 1 0 0
## 3 2 4 3 3 2 2 3 1 2 2
## 4 1 2 2 1 1 1 1 1 0 -1
## 5 NA NA NA NA NA NA NA NA NA NA
## 6 NA NA NA NA NA NA NA NA NA NA
## 7 NA NA NA NA NA NA NA NA NA NA
## 8 2 4 3 2 2 3 2 2 0 1
## 9 3 4 3 2 1 2 1 1 0 1
## 10 NA NA NA NA NA NA NA NA NA NA
## tfros.2 chal.2 prev.2 mfrsub.1 nfrsubgt3
## 1 1 NA NA 2.400000 0
## 2 0 NA NA 2.916667 4
## 3 2 NA NA 2.600000 1
## 4 0 NA NA NaN 0
## 5 NA NA NA 3.000000 0
## 6 NA NA NA 3.000000 2
## 7 NA NA NA 3.333333 1
## 8 1 NA NA 2.428571 0
## 9 1 NA NA 3.125000 2
## 10 NA NA NA 3.500000 1
6. Deal with friendship network (W)
sem1graph = graph.data.frame(d = sem1[1:2], directed=TRUE)
sem1matrix = get.adjacency(sem1graph)
sem1matrix = as.matrix(sem1matrix )
attitudesw = attitudesw[match(row.names(sem1matrix), attitudesw$std_id),]
atts = attitudesw[!is.na(attitudesw$sub.2),]
atts = atts[!is.na(atts$sub.1),]
atts = atts[!is.na(atts$mfrsub.1),]
W = sem1matrix
# 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))]
# make sure (this will return 0 if we did it right):
which(rownames(W) != atts$std_id)
## integer(0)
7. Fit a Linear Network Autocorrelation Model
# Linear Network Autocorrelation Model.
pim1<-lnam(atts$sub.2,cbind(atts$sub.1, atts$mfrsub.1), W2 = W)
summary(pim1)
##
## Call:
## lnam(y = atts$sub.2, x = cbind(atts$sub.1, atts$mfrsub.1), W2 = W)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.4678 -0.4678 0.1391 0.4455 2.0064
##
## Coefficients:
## Estimate Std. Error Z value Pr(>|z|)
## X1 0.60696 0.04019 15.103 < 2e-16 ***
## X2 0.34667 0.04232 8.191 2.22e-16 ***
## rho2.1 0.02887 0.01691 1.707 0.0878 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Estimate Std. Error
## Sigma 0.7247 0.001
##
## Goodness-of-Fit:
## Residual standard error: 0.7322 on 342 degrees of freedom (w/o Sigma)
## Multiple R-Squared: 0.4136, Adjusted R-Squared: 0.4085
## Model log likelihood: -379 on 341 degrees of freedom (w/Sigma)
## AIC: 766.1 BIC: 781.4
##
## Null model: meanstd
## Null log likelihood: -476.3 on 343 degrees of freedom
## AIC: 956.5 BIC: 964.2
## AIC difference (model versus null): 190.4
## Heuristic Log Bayes Factor (model versus null): 182.8