## ------------------------------------------------------------------------
#### ---- obtain trade data from top 30 countries in terms of GDP
data(IR90s)
gdp<-IR90s$nodevars[,2]
topgdp<-which(gdp>=sort(gdp,decreasing=TRUE)[30] )
Y<-log( IR90s$dyadvars[topgdp,topgdp,2] + 1 )
Y[1:5,1:5]
## ------------------------------------------------------------------------
#### ---- ANOVA for trade data
Rowcountry<-matrix(rownames(Y),nrow(Y),ncol(Y))
Colcountry<-t(Rowcountry)
anova(lm( c(Y) ~ c(Rowcountry) + c(Colcountry) ) )
## ------------------------------------------------------------------------
#### ---- comparison of countries in terms of row and column means
rmean<-rowMeans(Y,na.rm=TRUE) ; cmean<-colMeans(Y,na.rm=TRUE)
muhat<-mean(Y,na.rm=TRUE)
ahat<-rmean-muhat
bhat<-cmean-muhat
# additive "exporter" effects
head( sort(ahat,decreasing=TRUE) )
# additive "importer" effects
head( sort(bhat,decreasing=TRUE) )
## ------------------------------------------------------------------------
#### ---- covariance and correlation between row and column effects
cov( cbind(ahat,bhat) )
cor( ahat, bhat)
## ------------------------------------------------------------------------
#### ---- an estimate of dyadic covariance and correlation
R <- Y - ( muhat + outer(ahat,bhat,"+") )
cov( cbind( c(R),c(t(R)) ), use="complete")
cor( c(R),c(t(R)), use="complete")
## ----echo=FALSE,fig.height=3,fig.width=6,out.width='6in',fig.align='center'----
par(mfrow=c(1,1),mar=c(3,3,1,1),mgp=c(1.75,.75,0))
plot(rmean*1.1,cmean*1.1,type="n",xlab="average log exports",
ylab="average log imports")
abline(h=0,col="gray",lty=2) ; abline(v=0,col="gray",lty=2)
abline(0,1,col="gray",lty=2)
text(rmean,cmean,dimnames(Y)[[1]],cex=.7)
## ----IR90s_fit_srm,cache=TRUE,fig.keep='last',results='hide',fig.cap='Default plots generated by the ame command.',fig.height=6,fig.width=5,fig.align='center'----
fit_SRM<-ame(Y)
## ------------------------------------------------------------------------
muhat # empirical overall mean
mean(fit_SRM$BETA) # model-based estimate
cov( cbind(ahat,bhat)) # empirical row/column mean covariance
apply(fit_SRM$VC[,1:3],2,mean) # model-based estimate
cor( c(R), c(t(R)) ,use="complete") # empirical residual dyadic correlation
mean(fit_SRM$VC[,4] ) # model-based estimate
## ----echo=FALSE,fig.height=2.7,fig.width=5.2,out.width='5.2in',fig.align='center'----
par(mfrow=c(1,2),mar=c(3,3,1,1),mgp=c(1.75,.75,0))
plot( ahat, fit_SRM$APM,xlab=expression(paste("ANOVA ",hat(a)) ) ,
ylab=expression(paste("Bayes ",hat(a)) ) )
abline(0,1)
plot( bhat, fit_SRM$BPM,xlab=expression(paste("ANOVA ",hat(b)) ) ,
ylab=expression(paste("Bayes ",hat(b)) ) )
abline(0,1)
## ------------------------------------------------------------------------
#### ---- nodal covariates
dimnames(IR90s$nodevars)[[2]]
Xn<-IR90s$nodevars[topgdp,]
Xn[,1:2]<-log(Xn[,1:2])
#### ---- dyadic covariates
dimnames(IR90s$dyadvars)[[3]]
Xd<-IR90s$dyadvars[topgdp,topgdp,c(1,3,4,5)]
Xd[,,3]<-log(Xd[,,3])
## ----IR90s_fit_srrm,cache=TRUE,fig.keep='none',results='hide'------------
fit_srrm<-ame(Y,Xd=Xd,Xr=Xn,Xc=Xn)
## ------------------------------------------------------------------------
summary(fit_srrm)
## ----IR90s_fit_rm,cache=TRUE,fig.keep='none',results='hide'--------------
fit_rm<-ame(Y,Xd=Xd,Xr=Xn,Xc=Xn,rvar=FALSE,cvar=FALSE,dcor=FALSE)
## ------------------------------------------------------------------------
summary(fit_rm)
## ----echo=FALSE,fig.height=3,fig.width=6,out.width='6in',fig.align='center'----
ht<-c(30,30,50,30)
GOF<-gofstats(Y)
par(mfrow=c(2,2),mar=c(3,3,1,1),mgp=c(1.75,.75,0))
for(k in 1:4)
{
xlim<-range(c(fit_rm$GOF[,k],fit_srrm$GOF[,k],GOF[k]))*c(.9,1.1)
hist(fit_rm$GOF[-1,k],prob=TRUE,col="pink",xlim=xlim,main="",ylab="",
xlab=names(GOF)[k],ylim=c(0,ht[k]))
clr<-c(col2rgb("lightblue")/255,.75)
hist(fit_srrm$GOF[-1,k],prob=TRUE,add=TRUE,
col=rgb(clr[1],clr[2],clr[3],clr[4]))
abline(v=GOF[k],col="red")
}
## ----cache=TRUE,fig.keep='none',results='hide'---------------------------
fit_srrm0<-ame(Y,Xd[,,1:2],Xn,Xn)
## ----echo=FALSE,fig.height=3,fig.width=6---------------------------------
hist(fit_srrm0$GOF[-1,4],main="",yaxt="n",ylab="",prob=TRUE,col="lightgreen",breaks=12,
xlab=names(GOF)[4] ,xlim=range(fit_srrm0$GOF[,4]),ylim=c(0,17) )
hist(fit_srrm$GOF[-1,4],prob=TRUE,add=TRUE,
col=rgb(clr[1],clr[2],clr[3],clr[4]))
abline(v=GOF[4],col="red")
## ----trade_ame2,cache=TRUE,fig.keep='last',results='hide',fig.cap='Diagnostic plots for the rank-2 AME model.',fig.height=6,fig.width=5,fig.align='center'----
fit_ame2<-ame(Y,Xd,Xn,Xn,R=2)
## ------------------------------------------------------------------------
summary(fit_ame2)
## ----echo=FALSE,fig.height=4.75,fig.width=4.75,fig.align='center'--------
mu<-mean(fit_ame2$BETA[,1])
br<-apply(fit_ame2$BETA[,1+1:3] ,2,mean)
bc<-apply(fit_ame2$BETA[,1+3+1:3],2,mean)
bd<-apply(fit_ame2$BETA[,1+6+1:4] ,2,mean)
EY2<- mu+
Xbeta(Xd,bd) +
( Xn%*%br + fit_ame2$APM)%*%rep(1,nrow(Y)) +
t( ( Xn%*%bc + fit_ame2$BPM)%*%rep(1,nrow(Y)) )
par(mar=c(0,0,0,0),mgp=c(0,0,0))
circplot(1*( Y - EY2 > 0 ) , fit_ame2$U , fit_ame2$V ,pscale=1.1 )
## ----llaw_graph,fig.height=4.75,fig.width=4.75,fig.align='center',fig.cap='Graph of the friendship network between 71 lawyers. Node colors represent at which of the three offices each lawyer works.'----
data(lazegalaw)
Y<-lazegalaw$Y[,,2]
Xd<-lazegalaw$Y[,,-2]
Xn<-lazegalaw$X
dimnames(Xd)[[3]]
dimnames(Xn)[[2]]
netplot(lazegalaw$Y[,,2],ncol=Xn[,3])
## ----cache=TRUE,fig.keep='none',results='hide'---------------------------
fit_SRM<-ame(Y,model="bin")
## ----cache=TRUE,fig.keep='none',results='hide'---------------------------
fit_SRG<-ame(Y,model="bin",rvar=FALSE,cvar=FALSE,dcor=FALSE)
## ----echo=FALSE,fig.height=3,fig.width=6,fig.align='center'--------------
gY<-gofstats(Y)
clr<-c(col2rgb("lightblue")/255,.75)
par(mfrow=c(2,2),mar=c(3,3,1,1),mgp=c(1.75,.75,0))
for(k in 1:4)
{
hist(fit_SRG$GOF[-1,k],main="",yaxt="n",ylab="",prob=TRUE,col="lightgreen",
xlab=names(gY)[k] ,xlim=range(c(fit_SRM$GOF[,k],fit_SRG$GOF[,k]) ) )
hist(fit_SRM$GOF[-1,k],prob=TRUE,add=TRUE,
col=rgb(clr[1],clr[2],clr[3],clr[4]))
abline(v=gY[k],col="red")
}
## ----echo=FALSE,fig.height=3,fig.width=6,fig.align='center'--------------
par(mfrow=c(1,2),mar=c(3,3,1,1),mgp=c(1.75,.75,0))
plot(rowSums(Y,na.rm=TRUE), fit_SRM$APM,xlab="outdegrees",
ylab="additive row effects")
plot( colSums(Y,na.rm=TRUE), fit_SRM$BPM,xlab="indegrees",
ylab="additive column effects")
## ----cache=TRUE,fig.keep='none',results='hide'---------------------------
Xno<-Xn[,c(1,2,4,5,6)]
fit_SRRM<-ame(Y, Xd=Xd, Xr=Xno, Xc=Xno, model="bin")
## ------------------------------------------------------------------------
summary(fit_SRRM)
## ----cache=TRUE,fig.keep='none',results='hide'---------------------------
fit_AME<-ame(Y, Xd=Xd[,,2], R=3, model="bin")
## ----echo=FALSE,fig.height=3,fig.width=6,fig.align='center'--------------
par(mfrow=c(2,2),mar=c(3,3,1,1),mgp=c(1.75,.75,0))
for(k in 1:4)
{
hist(fit_AME$GOF[-1,k],main="",yaxt="n",ylab="",prob=TRUE,col="lightblue",
xlab=names(gY)[k] ,xlim=range(fit_AME$GOF[,k]) )
abline(v=gY[k],col="red")
}
## ------------------------------------------------------------------------
U<-fit_AME$U
V<-fit_AME$V
round(cor(U, Xno),2)
round(cor(V, Xno),2)
## ----echo=FALSE,fig.height=3,fig.width=8.25------------------------------
par(mfrow=c(1,3),mar=c(3,3,1,1),mgp=c(1.75,.75,0))
xlb<-c(expression(u[1]),expression(u[2]),expression(u[3]) )
ylb<-c(expression(v[1]),expression(v[2]),expression(v[3]) )
par(mfcol=c(1,3))
for(k in 1:3)
{
plot(U[,k],V[,k] , xlab=xlb[k],ylab=ylb[k],
pch=1+(Xn[,1]-1)*15,
col=Xn[,3]+2 ,
xlim=range(U),ylim=range(V))
abline(h=0,lty=2,col="gray") ; abline(v=0,lty=2,col="gray")
}
## ------------------------------------------------------------------------
data(sheep)
Y<-sheep$dom
gofstats(Y)
## ----echo=FALSE,fig.height=2.75,fig.width=7.5,fig.align='center'---------
par(mfrow=c(1,3),mar=c(3,3,1,1),mgp=c(1.75,.75,0))
hist(Y,main="",col="brown",xlab="number of encounters",ylab="",prob=TRUE)
plot(sheep$age,apply(Y,1,mean,na.rm=TRUE ),xlab="age",ylab="row mean")
plot(sheep$age,apply(Y,2,mean,na.rm=TRUE ),xlab="age",ylab="column mean" )
## ----cache=TRUE,fig.keep='none',results='hide'---------------------------
x<-sheep$age - mean(sheep$age)
Xd<-outer(x,x)
Xn<-cbind(x,x^2) ; colnames(Xn)<-c("age","age2")
fit<-ame(Y, Xd, Xn, Xn, model="ord")
## ------------------------------------------------------------------------
summary(fit)
## ------------------------------------------------------------------------
Y<-sampsonmonks[,,3]
apply(Y>0,1,sum,na.rm=T)
## ----cache=TRUE,fig.keep='none',results='hide'---------------------------
odmax<-rep(3,nrow(Y))
odmax[ apply(Y>0,1,sum,na.rm=T)>3 ]<-4
fit<-ame(Y,R=2,model="frn",odmax=odmax)
## ------------------------------------------------------------------------
summary(fit)
## ----echo=FALSE,fig.height=6,fig.width=5,fig.align='center'--------------
plot(fit)
## ----cache=TRUE,fig.keep='last',results='hide'---------------------------
data(dutchcollege)
Y<-1*( dutchcollege$Y[,,7] > 1 ) # indicator of positive relationship at the last timepoint
Xn<-dutchcollege$X[,1] # nodal indicator of male sex
Xd<-1*(outer(Xn,Xn,"==")) # dyadic indicator of same sex
## ----echo=FALSE----------------------------------------------------------
set.seed(3)
## ------------------------------------------------------------------------
n<-nrow(Y)
Ys<-matrix(NA,n,n) # sociomatrix for sampled data
egos<-sort(sample(n,5)) # ego sample
Ys[egos,]<-Y[egos,] # relations of egos are observed
for(i in egos)
{
ai<-which(Ys[i,]==1) # alters of i
Ys[ai,ai]<-Y[ai,ai] # relations between alters of i are observed
}
mean(is.na(Ys))
## ------------------------------------------------------------------------
egos
Ys[1:10,1:10]
## ----cache=TRUE,fig.keep='none',results='hide'---------------------------
fit_pop<-ame(Y,Xd,Xn,Xn,model="bin") # fit based on full data (population)
fit_ess<-ame(Ys,Xd,Xn,Xn,model="bin") # fit based on egocentric subsample
## ------------------------------------------------------------------------
apply(fit_pop$BETA,2,mean)
apply(fit_ess$BETA,2,mean)
## ------------------------------------------------------------------------
miss<-which(is.na(Ys))
mean( ( fit_pop$YPM[miss] - Y[miss] )^2, na.rm=TRUE )
mean( ( fit_ess$YPM[miss] - Y[miss] )^2, na.rm=TRUE )
## ----mdsimstudy,cache=TRUE,echo=FALSE,results='hide'---------------------
BETA_SS<-NULL
for(sim in 1:100)
{
set.seed(sim)
Ys<-matrix(NA,n,n) # sociomatrix for sampled data
s<-sort(sample(n,5)) # ego sample
Ys[s,]<-Y[s,] # relations of egos are observed
for(i in s)
{
ai<-which(Ys[i,]==1) # alters of i
Ys[ai,ai]<-Y[ai,ai] # relations between alters of i are observed
}
fit<-ame(Ys,Xd,Xn,Xn,nscan=5000,model="bin",print=FALSE,plot=FALSE,gof=FALSE)
BETA_SS<-rbind(BETA_SS,apply(fit$BETA,2,mean))
cat(sim,"\n")
}
## ----echo=FALSE,fig.height=4,fig.width=7,fig.align='center'--------------
cnames<-c(expression(beta[0]),expression(beta[r]), expression(beta[c]),
expression(beta[d]))
beta_pop<-apply(fit_pop$BETA,2,mean)
par(mfrow=c(2,2),mar=c(3,3,1,1),mgp=c(1.75,.75,0))
for(k in 1:4)
{
hist( BETA_SS[,k] ,col="lightgreen", main="",ylab="",xlab=cnames[k],yaxt="n",
xlim=1.2*range(c(beta_pop[k],BETA_SS[,k])) )
abline(v=beta_pop[k],lwd=2,col="blue")
}
## ----echo=FALSE,fig.height=3,fig.width=7---------------------------------
data(dutchcollege)
Y<-1*( dutchcollege$Y >= 2 )
xnode<-xnet(Y[,,dim(Y)[3]],fm=FALSE)
xnode[18,]<-.5*( xnode[18,] + apply(xnode[-18,],2,mean))
par(mar=c(1,1,1.25,1),mgp=c(1.75,75,0),mfrow=c(2,4))
for(k in 1:dim(Y)[3])
{
netplot(Y[,,k],xnode)
mtext(side=3,paste("Time",k-1))
}
## ------------------------------------------------------------------------
# outcome
Y<-1*( dutchcollege$Y >= 2 )[,,2:7]
n<-dim(Y)[1] ; t<-dim(Y)[3]
# nodal covariates
Xnode<-dutchcollege$X[,1:2] # sex and smoking status
Xnode<-array(Xnode,dim=c(n,ncol(Xnode),t))
dimnames(Xnode)[[2]]<-c("male","smoker")
# dyadic covariates
Xdyad<-array(dim=c(n,n,5,t))
Xdyad[,,1,]<-1*( dutchcollege$Y >= 2 )[,,1:6] # lagged value
Xdyad[,,2,]<-array(apply(Xdyad[,,1,],3,t),dim=c(n,n,t)) # lagged reciprocal value
Xdyad[,,3,]<-tcrossprod(Xnode[,1,1]) # both male
Xdyad[,,4,]<-tcrossprod(Xnode[,2,1]) # both smokers
Xdyad[,,5,]<-outer( dutchcollege$X[,3],dutchcollege$X[,3],"==") # same program
dimnames(Xdyad)[[3]]<-c("Ylag","tYlag","bothmale","bothsmoke","sameprog")
## ----dc_ar1, cache=TRUE,fig.keep='none',results='hide'-------------------
fit_ar1<-ame_rep(Y,Xdyad,Xnode,Xnode,model="bin",plot=FALSE)
## ------------------------------------------------------------------------
summary(fit_ar1)
## ------------------------------------------------------------------------
Wnode<-Xnode
Wnode[,,1:3]<-0
XWnode<-array( dim=dim(Xnode)+c(0,2,0))
XWnode[,1:2,]<-Xnode ; XWnode[,3:4,]<-Wnode
dimnames(XWnode)[[2]]<-c(dimnames(Xnode)[[2]],paste0(dimnames(Xnode)[[2]],".w"))
Wdyad<-Xdyad
Wdyad[,,,1:3]<-0
XWdyad<-array( dim=dim(Xdyad)+c(0,0,5,0) )
XWdyad[,,1:5,]<-Xdyad ; XWdyad[,,6:10,]<-Wdyad
dimnames(XWdyad)[[3]]<-c(dimnames(Xdyad)[[3]],paste0(dimnames(Xdyad)[[3]],".w"))
## ----dc_ar1_vb, cache=TRUE,fig.keep='none',results='hide'----------------
fit_ar1_vb<-ame_rep(Y,XWdyad,XWnode,XWnode,model="bin")
## ------------------------------------------------------------------------
summary(fit_ar1_vb)
## ------------------------------------------------------------------------
data(coldwar)
# response
Y<-sign( apply(coldwar$cc,c(1,2), mean ) )
# nodal covariates
Xn<-cbind( apply( log(coldwar$gdp),1,mean ) , # log gdp
sign(apply(coldwar$polity ,1,mean ) ) ) # sign of polity
Xn[,1]<-Xn[,1]-mean(Xn[,1])
dimnames(Xn)[[2]]<-c("lgdp","polity")
# dyadic covariates
Xd<-array(dim=c(nrow(Y),nrow(Y),3))
Xd[,,1]<- tcrossprod(Xn[,1]) # gdp interaction
Xd[,,2]<- tcrossprod(Xn[,2]) # polity interaction
Xd[,,3]<-log(coldwar$distance) # log distance
dimnames(Xd)[[3]]<-c("igdp","ipol","ldist")
## ----cwfitR1,cache=TRUE,fig.keep='none',results='hide'-------------------
fit_cw_R1<-ame(Y,Xd,Xn,R=1,model="ord",symmetric=TRUE,burn=1000,nscan=100000,odens=100)
## ----echo=FALSE----------------------------------------------------------
plot(fit_cw_R1)
#plot(fit_cw_R2)
## ------------------------------------------------------------------------
summary(fit_cw_R1)
fit_cw_R1$L # eigenvalues
## ----fig.height=5,fig.width=7,echo=FALSE---------------------------------
u<-fit_cw_R1$U
cnames<-rownames(u)
urank<-rank(u)
plot(urank,u,type="n",xlab="rank order of u")
abline(h=0,col="gray")
addlines(Y>0 , cbind(urank,u),col="green")
addlines(Y<0 , cbind(urank,u),col="pink")
text(urank,u,cnames,srt=-45,cex=1.0)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.