# inst/doc/amen.R In pdhoff/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[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

## ------------------------------------------------------------------------
#### ---- 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,family="nrm")

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

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,family="nrm")

## ------------------------------------------------------------------------
summary(fit_srrm)

## ----IR90s_fit_rm,cache=TRUE,fig.keep='none',results='hide'--------------
fit_rm<-ame(Y,Xd=Xd,Xr=Xn,Xc=Xn,family="nrm",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)
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,family="nrm")

## ----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,20) )
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,family="nrm",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,family="bin")

## ----cache=TRUE,fig.keep='none',results='hide'---------------------------
fit_SRG<-ame(Y,family="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]) ) )
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",
plot( colSums(Y,na.rm=TRUE), fit_SRM\$BPM,xlab="indegrees",

## ----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, family="bin")

## ------------------------------------------------------------------------
summary(fit_SRRM)

## ----cache=TRUE,fig.keep='none',results='hide'---------------------------
fit_AME<-ame(Y, Xd=Xd[,,2], R=3, family="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, family="ord")

## ------------------------------------------------------------------------
summary(fit)

## ------------------------------------------------------------------------
Y<-sampsonmonks[,,3]

apply(Y>0,1,sum,na.rm=T)

## ----cache=TRUE,fig.keep='none',results='hide'---------------------------
fit<-ame(Y,R=2,family="frn",odmax=4)

## ------------------------------------------------------------------------
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,family="bin")   # fit based on full data (population)

fit_ess<-ame(Ys,Xd,Xn,Xn,family="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,family="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")

Xdyad[,,1,]<-1*( dutchcollege\$Y >= 2 )[,,1:6]                     # lagged value

## ----dc_ar1, cache=TRUE,fig.keep='none',results='hide'-------------------

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

## ----dc_ar1_vb, cache=TRUE,fig.keep='none',results='hide'----------------

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

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,family="ord",symmetric=TRUE,burn=1000,nscan=100000,odens=50)

## ----echo=FALSE----------------------------------------------------------
plot(fit_cw_R1)

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