## File Name: SRM_ARBSRM_BOND_MALLOY.R
## File Version: 0.26
## This function contains the code of
## Bond & Malloy (2018). ARBSRM - The General Social Relations Model
## http://thomasemalloy.org/arbsrm-the-general-social-relations-model/
## copied and slightly edited from the file
## http://thomasemalloy.org/wp-content/uploads/2017/09/arbcodeR.pdf
SRM_ARBSRM_BOND_MALLOY <- function(data, serror=TRUE)
{
CALL <- match.call()
s1 <- Sys.time()
miss <- -999
data[ is.na(data) ] <- miss
nact <- nrow(data)
npar <- ncol(data)
ntot <- 0
for (i in 1:nact) {
for (j in 1:npar) {
if (data[i,j]!=miss){
ntot<-ntot+1
}
}
}
inde <- matrix(0,ntot/2,2)
dv <- rep(0,ntot)
k <- 0
l <- nact/2+1
bata <- matrix(0,nact,npar)
for (i in 1:nact) {
for (j in 1:npar) {
if (data[i,j]!=miss){
k<-k+1
dv[k] <- data[i,j]
bata[i,j] <- 1
if (i<l) {
inde[k,1]<-i #For each data point, first
inde[k,2]<-j
}
} #and second subscripts.
}
}
ntot <- k/2
nact <- nact/2
ac <- matrix(0,ntot,nact)
pa <- matrix(0,ntot,npar)
for (k in 1:ntot){
ac[k,inde[k,1]]<-1 # indicator matrices for
pa[k,inde[k,2]]<-1
} # actor and partner
apar <- colSums(ac)
pact <- colSums(pa)
mirr <- rep(0,nact)
recip <- matrix(0,nact,npar)
aa <- matrix(0,ntot,(6*ntot)) #matrix representation of SCP
vv <- matrix(0,ntot,(6*ntot)) #var-cov matrix of observed scores
uvv <- matrix(0,ntot,ntot)
j0 <- c(1:(6*ntot))
j0[(2*ntot+1):(3*ntot)]<-c((3*ntot+1):(4*ntot))
j0[(3*ntot+1):(4*ntot)]<-c((2*ntot+1):(3*ntot))
for (i in 1:ntot){
aa[i,i] <- 1
vv[i,i] <- 1
ii <- ntot+i
aa[i,ii] <- 1
vv[i,ii] <- 1
ii <- 4*ntot+i
aa[i,ii] <- 1
vv[i,ii] <- 1
}
for (i in 1:ntot){
for (j in (4*ntot+1):(5*ntot)){
aa[i,j]<-aa[i,j]-1/ntot
}
}
for (i in 1:ntot) {
aa[i,i]<-1-1/pact[inde[i,2]]
aa[i,(ntot+i)]<-1-1/apar[inde[i,1]]
for (j in 1:ntot) {
if ((inde[i,1]!=inde[j,1]) | (inde[i,2]!=inde[j,2])) {
if ((inde[i,1]==inde[j,2]) & (inde[i,2]==inde[j,1])){
recip[inde[i,1],inde[i,2]] <- 1
aa[i,(5*ntot+i)]<-1
aa[i,(5*ntot+j)]<--1
vv[i,(2*ntot+j)]<-1
vv[i,(3*ntot+j)]<-1
vv[i,(5*ntot+j)]<-1
}
jj<-j+ntot
if (inde[i,1]==inde[j,1]) aa[i,jj]<--1/apar[inde[i,1]]
if (inde[i,1]==inde[j,1]) vv[i,j]<-1
if (inde[i,2]==inde[j,2]) vv[i,jj]<-1
if (inde[i,2]==inde[j,2]) aa[i,j]<--1/pact[inde[i,2]]
if (inde[i,1]==inde[j,2]) vv[i,(2*ntot+j)]<-1
if (inde[i,2]==inde[j,1]) vv[i,(3*ntot+j)]<-1
}
}
}
mirr <- rowSums(recip)
for (i in 1:ntot){
for (j in 1:ntot){
if ((recip[inde[i,1],inde[i,2]]!=0) & (recip[inde[j,1],inde[j,2]]!=0)){
if ((inde[i,1]==inde[j,2]) & (inde[i,2]==inde[j,1])){
aa[i,(3*ntot+j)]<-(mirr[inde[i,1]]-1)/mirr[inde[i,1]]
aa[i,(2*ntot+j)]<-(mirr[inde[j,1]]-1)/mirr[inde[j,1]]
}
if ((inde[i,1]!=inde[j,2]) & (inde[i,2]==inde[j,1])) aa[i,(2*ntot+j)]<--1/mirr[inde[j,1]]
if ((inde[i,1]==inde[j,2]) & (inde[i,2]!=inde[j,1])) aa[i,(3*ntot+j)]<--1/mirr[inde[i,1]]
}
}
}
uvv <- vv[,(2*ntot+1):(3*ntot)]+vv[,(3*ntot+1):(4*ntot)]
ss <- matrix(0,6,3) # Sums of Squares and Cross-Products
c <- matrix(0,6,6) # Searle, Casella, & McCullough, eq 18, p. 173
w <- matrix(0,6,276) # Var-Covar matrix for SS and SCP
uw <- matrix(0,6,126)
zz <- matrix(0,ntot,ntot)
uu <- matrix(0,ntot,ntot)
for (i in 1:6){
dumi<-((i-1)*ntot+1):(i*ntot)
if ((i<3) | (i>4)){
for (j in 1:6){
dumj<-((j-1)*ntot+1):(j*ntot)
c[i,j]<-sum(diag(aa[,dumi]%*%vv[,dumj]))
}
}
if ((i==3) | (i==4)){
zz<-aa[,dumi]
for (j in 1:6) {
dumj<-((j-1)*ntot+1):(j*ntot)
uu<-vv[,dumj]
c[i,j]<-.5*(sum(diag(zz%*%t(uu))) + sum(diag((t(zz)%*%uu))) )
}
}
}
cin <- solve(c)
for (i in 1:6){
dumi<-((i-1)*ntot+1):(i*ntot)
ss[i,1] <- t(dv[1:ntot])%*%aa[,dumi]%*%dv[1:ntot]
ss[i,2] <- t(dv[(ntot+1):(2*ntot)])%*%aa[,dumi]%*%dv[(ntot+1):(2*ntot)]
ss[i,3]<-t(dv[1:ntot])%*%aa[,dumi]%*%dv[(ntot+1):(2*ntot)]
}
est <- matrix(0,6,3)
est <- cin %*% ss # Here are the parameter estimates
s2 <- Sys.time()
#--- standard errors
se <- NA*est
if (serror) {
cwu=matrix(0,6,6)
j=1
foo=c(0,6,12,18,24,30)
for (i in foo){
cwu[1,j]=i
cwu[j,1]=i
j=j+1
}
j=2
foo=c(36,42,48,54,60)
for (i in foo){
cwu[2,j]=i
cwu[j,2]=i
j=j+1
}
j=3
foo=c(66,72,78,84)
for (i in foo){
cwu[3,j]=i
cwu[j,3]=i
j=j+1
}
j=4
foo=c(90,96,102)
for (i in foo){
cwu[4,j]=i
cwu[j,4]=i
j=j+1
}
cwu[5,5]=108
cwu[6,5]=114 # an accounting variable,
cwu[5,6]=114 # cwu will point us to the
cwu[6,6]=120 # right parms, SS, and SCP
# Here comes the time-consuming computation
for (i1 in 1:6){
dumi1=(ntot*(i1-1)+1):(i1*ntot)
w1=0
for (j1 in 1:6){
startw=w1
dumj1=(ntot*(j1-1)+1):(j1*ntot)
garmat1<- aa[,dumi1] %*% vv[,j0[dumj1]]
garmat2<- aa[,j0[dumi1]] %*% vv[,dumj1]
garmat3<- aa[,dumi1] %*% vv[,dumj1]
for (i2 in 1:6){
dumi2=(ntot*(i2-1)+1):(i2*ntot)
w1<-startw+i2-6
barmat1<- garmat1 %*% aa[,dumi2]
barmat2<- garmat2 %*% aa[,j0[dumi2]]
barmat3<- garmat3 %*% aa[,dumi2]
for (j2 in j1:6){
w1=w1+6
dumj2=(ntot*(j2-1)+1):(j2*ntot)
cons=.5
if (j2!=j1) cons=1
w[i1,w1]=w[i1,w1]+cons*sum(diag(barmat1 %*% vv[,j0[dumj2]]))
w[i1,w1]=w[i1,w1]+cons*sum(diag(barmat2 %*% vv[,dumj2]))
}
for (j2 in 1:6){ # setting up var-cov
jj2=j2 # matrix for SS & SCP
if ((j2==3) | (j2==4)) jj2=7-j2
dumj2=(ntot*(jj2-1)+1):(jj2*ntot)
w2=cwu[j1,j2]+i2
uw[i1,w2]=uw[i1,w2]+2*sum(diag(barmat3 %*% vv[,dumj2]))
}
} # end i2 loop
} # end j1 loop
for (i2 in 1:6){
dumi2=(ntot*(i2-1)+1):(i2*ntot)
w1=i2+120
for (j1 in 1:5) {
jj1=j1
if (j1>3) jj1=j1+1
uu=vv[,(ntot*(jj1-1)+1):(jj1*ntot)]
if (j1==3) uu=uvv
for (j2 in 1:5) {
w1=w1+6
jj2=j2
if (jj2>3) jj2=j2+1
zz=vv[,(ntot*(jj2-1)+1):(jj2*ntot)]
if (j2==3) zz=uvv
w[i1,w1]=w[i1,w1]+sum(diag(aa[,j0[dumi1]]%*%uu%*%aa[,dumi2]%*%zz))
}
}
}
} # End i2 and i1 loops and end time-consuming computation
b=matrix(0,46,46)
var=matrix(0,6,6)
for (j1 in 1:46) {
dum=((j1-1)*6+1):(j1*6)
var=cin%*%w[,dum]%*%t(cin) # Searle, Casella, & McCullough
r1=0 # Equation 21, page 176
for (i1 in 1:6){
for (i2 in i1:6) {
r1=r1+1
b[r1,j1]=var[i1,i2]
}
}
if (j1<22){
r2=21
var=cin%*%uw[,dum]%*%t(cin)
for (i1 in 1:5) {
ii1=i1
if (i1>3) ii1=i1+1
for (i2 in 1:5){
ii2=i2
if (i2>3) ii2=i2+1 # Accommodating the
r2=r2+1 # asymmetrical SCPs
bare=var[ii1,ii2]
if (i1==3) bare=(var[3,ii2]+var[4,ii2])/2
if (i2==3) bare=(var[ii1,3]+var[ii1,4])/2
if ((i1==3) & (i2==3))
bare=(var[3,3]+var[3,4]+var[4,3]+var[4,4])/4
b[r2,j1]=bare
}
} # b gives exact vars and covs # among parameter estimates
}
}
bi=diag(46)
bi=b+bi
bi=solve(bi) # the biggest inversion
bb<-matrix(0,6,46)
dn<-c(1,7,12,16,19,21)
for (i in 1:6) {
bb[i,]<-b[dn[i],]
}
bb= bb %*% bi
# Having set everything up,
# we now begin using the data
oot<-matrix(0,6,3)
kk=3
estp<-rep(0,46)
#Computing products of estimates
# via the R operator for outer product: %o%
xvec<-c(est[1,1],est[2,1],est[3,1],est[5,1],est[6,1])
yvec<-c(est[1,2],est[2,2],est[3,2],est[5,2],est[6,2])
for (kk in 1:3){
xyvec<- c(est[1,kk],est[2,kk],est[3,kk],est[4,kk],est[5,kk],est[6,kk])
matr<-xyvec %o% xyvec
estp[1:21]<- c(matr[,1],matr[2:6,2],matr[3:6,3],matr[4:6,4],matr[5:6,5],matr[6,6])
matrx<- xvec %o% xvec
if (kk==2) matrx<- yvec %o% yvec
if (kk==3) matrx<- yvec %o% xvec
estp[22:46]<-c(matrx)
oot[,kk] = bb %*% estp
}
oot = oot**.5 # converting var to se
se = oot
} # oot gives estimated variances and covs
#--- output
s3 <- Sys.time()
time <- list(start=s1, end=s3, time_est=s2-s1, time_se=s3-s2)
res <- list(est=est, se=se, serror=serror, time=time)
return(res)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.