inst/doc/amen.R

In winnga/amen: Additive and Multiplicative Effects Models for Networks and Relational Data

## ------------------------------------------------------------------------
#### ---- 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)