#######################################################################################
# ADDITIONAL R FUNCTIONS USED IN DYADA PACKAGE #
#######################################################################################
# Function to detect square matrices #
squaremat <- function (X){
sqmat <- (nrow(X) == ncol(X));
return(sqmat)
}
# Function to prepare matrices for rank orders without self-addressed cells #
prepmat <- function(X){
dyad <- ncol(X)-1;
X1 <- matrix(t(X)[row(X)!=col(X)],nrow=nrow(X),ncol=dyad,byrow=T);
return(X1);
}
#######################################################################################
# R FUNCTIONS FOR ANALYZING LINEAR HIERARCHY IN GROUPS #
#######################################################################################
# Function to transform original sociomatrix into a matrix of abilities #
freq2abil <- function (X,method=c('landau','matman','schmid')){
method <- match.arg(method)
if (method=='landau'){
M1 <- matrix(1*as.numeric(X>t(X)),nrow=nrow(X),byrow=F)
M2 <- matrix(-1*as.numeric(X<t(X)),nrow=nrow(X),byrow=F)
M3 <- matrix(0.*as.numeric(X==t(X)),nrow=nrow(X),byrow=T)
V <- M1 + M2 + M3
}
if (method=='matman'){
M1 <- matrix(1*as.numeric(X>t(X)),nrow=nrow(X),byrow=F)
M2 <- matrix(1./2.*as.numeric(X==t(X)),nrow=nrow(X),byrow=T)
V <- M1 +M2
}
if (method=='schmid'){
M1 <- matrix(1*as.numeric(X>t(X)),nrow=nrow(X),byrow=F)
M2 <- matrix(-1*as.numeric(X<t(X)),nrow=nrow(X),byrow=F)
M3 <- matrix(0.*as.numeric(X==t(X) & X+t(X)==0),nrow=nrow(X),byrow=T)
V <- M1 + M2 + M3
}
diag(V) <- 0.
dimnames(V) <- dimnames(X)
return(V)
}
# Function to obtain the number of tied relationships #
gettied <- function(X,relfreq=TRUE){
dyad <- nrow(X)*(nrow(X)-1)/2;
tied <- sum(X==t(X) & X!=0.)/2
if (relfreq == TRUE) list(tied.dyads=tied, perc.tied.dyads=(tied/dyad*100))
else return(tied)
}
# Function to obtain the number of unknown dyadic relationships #
getunknown <- function(X,relfreq=TRUE){
dyad <- nrow(X)*(nrow(X)-1)/2;
unknown <- (sum(X+t(X)==0.)-nrow(X))/2;
if (relfreq == TRUE) list(unknown.dyads=unknown,
perc.unknown.dyads=(unknown/dyad*100))
else return(unknown)
}
# Function to obtain the number of dyads with one-way relationships #
getonewayrel <- function (X, relfreq=TRUE){
dyad <- nrow(X)*(nrow(X)-1)/2;
c <- X + t(X);
onewayrel <- sum(X-c==0. & c!=0)
if (relfreq == TRUE) list(oneway.dyads=onewayrel,
perc.oneway.dyads=(onewayrel/dyad*100))
else return(onewayrel)
}
# Function to obtain the number of dyads with two-way relationships #
gettwowayrel <- function (X, relfreq=TRUE){
dyad <- nrow(X)*(nrow(X)-1)/2;
c <- X + t(X);
twowayrel <- sum((X-c!= 0.) * (t(X) - c != 0.))/2
if (relfreq == TRUE) list(twoway.dyads=twowayrel, perc.twoway.dyads=
(twowayrel/dyad*100))
else return(twowayrel)
}
getsignificant <- function(X,alpha=0.05, relfreq=TRUE){
dyad <- nrow(X)*(nrow(X)-1)/2;
pvals <- mapply(x=X[row(X)!=col(X) & X+t(X)!=0 & !is.na(X+t(X))],
n=(X+t(X))[row(X)!=col(X) & X+t(X)!=0& !is.na(X+t(X))],
FUN=function(x,n) binom.test(x,n)$p.value)
binommat <- matrix(NA,nrow=nrow(X),ncol(X))
binommat[row(X)!=col(X) & X+t(X)!=0 & !is.na(X+t(X))]<-round(pvals,3)
diag(binommat) <- 0
binommat[upper.tri(binommat)] <- (X+t(X))[upper.tri(X)]
ndyads <- sum(binommat[lower.tri(binommat)]<=alpha,na.rm=TRUE)
if (!is.null(dimnames(X))) dimnames(binommat) <- dimnames(X)
if (relfreq == TRUE) res<-list(original.mat=X,binomtest.matrix=binommat,significant.dyads=ndyads,
perc.significant.dyads=(ndyads/dyad*100))
else res <- list(original.mat=X,binomtest.matrix=binommat,significant.dyads=ndyads)
class(res) <- "sigdyads"
res
}
plot.sigdyads <- function(x){
#library(cowplot)
#library(reshape2)
#library(ggplot2)
X <- x$original.mat
pvals <- mapply(x=X[row(X)!=col(X) & X+t(X)!=0 & !is.na(X+t(X))],
n=(X+t(X))[row(X)!=col(X) & X+t(X)!=0 & !is.na(X+t(X))],
FUN=function(x,n) binom.test(x,n)$p.value)
temp <- matrix(NA,nrow=nrow(X),ncol(X))
temp[row(X)!=col(X) & X+t(X)!=0 & !is.na(X+t(X))]<-pvals
myfun<-function(x){if (!is.na(x) & x >=0.01){res <-round(x,2)} else if(!is.na(x) & x <0.01)
{res<-'<0.01'} else {res<-'NA'}; res}
pvals2<- sapply(t(temp)[upper.tri(temp)],function(x) myfun(x))
binommat <- matrix(NA,nrow=nrow(X),ncol(X))
binommat[row(X)!=col(X) & X+t(X)!=0 & !is.na(X+t(X))]<-pvals
binommat[row(binommat)>col(binommat)] <- (X+t(X))[lower.tri(X)]
dyadn <- paste(t(X),'/',X+t(X),sep='')[lower.tri(X)]
p <- ggplot(subset(melt(binommat)),
aes(x=Var1,y=Var2))
mi.ids <- subset(melt(binommat), Var1 == Var2)
mi.low <- subset(melt(binommat), Var1 > Var2)
mi.upp <- subset(melt(binommat), Var1 < Var2)
p1 <- (p + geom_tile(data=mi.ids, colour='black',alpha=0.2) +
geom_text(data=mi.ids,aes(label=if (is.null(rownames(X))){paste("Ind.",Var1)} else{
rownames(X)}),colour='red',fontface=2)+
geom_tile(data=mi.low,colour='black',alpha=0.1)+
geom_text(data=mi.low,aes(label=dyadn))+
geom_tile(data=mi.upp,aes(fill = value),colour = "black") +
geom_text(data=mi.upp,aes(label = pvals2)) +
scale_fill_gradient2(name="P-value",midpoint=0.5,low = "green",
mid='orange', high = "red",na.value = 'grey') +
xlab(NULL)+ylab(NULL)) + scale_y_reverse()
p1 + theme(axis.ticks = element_blank(), axis.text.y = element_blank()) +
theme(axis.ticks = element_blank(), axis.text.x = element_blank())
}
# Function to obtain several win and loss measures at individual level #
getwinloss <- function(X,names=NULL,method=c("Dij","Pij")){
if (!squaremat(X))
return("Error: Matrix X is not square and cannot be analyzed");
if ( is.na(X) || !is.numeric(X))
return("Error: Sociomatrix must be numeric");
method <- match.arg(method)
dyadc <- X + t(X);
if (method == "Dij"){
Dij <- X/dyadc-(((X/dyadc)-0.5)/(dyadc+1));
Dij[is.nan(Dij)] <- 0.;
w1 <- rowSums(Dij);
w2 <- Dij%*%w1;
l1 <-colSums(Dij);
l2 <- t(l1)%*%Dij;
}
if (method == "Pij"){
Pij <- array(dim=c(nrow(X),ncol(X)),0.);
Pij <- X/dyadc;
Pij[is.nan(Pij)] <- 0.
w1 <- rowSums(Pij);
w2 <- Pij%*%w1;
l1 <-colSums(Pij);
l2 <- t(l1)%*%Pij;
}
wl<- data.frame(w=array(w1,dim=c(nrow(X),1)),weighted.w=w2,l=array(l1,dim=c(nrow(X),1)),
weighted.l=t(l2))
if (is.null(names)) names <- paste('Ind.',1:nrow(X))
rownames(wl) <- names
return(wl)
}
# Function to obtain the David's scores (DS) based on dyadic dominance indices #
getDS <- function(X,names=NULL,method=c("Dij","Pij")){
if (!squaremat(X))
return("Error: Matrix X is not square and cannot be analyzed")
if ( is.na(X) || !is.numeric(X))
return("Error: Sociomatrix must be numeric")
method <- match.arg(method)
dyadc <- X + t(X)
if (method == "Dij"){
Dij <- array(dim=c(nrow(X),ncol(X)),0.)
Dij <- X/dyadc-(((X/dyadc)-0.5)/(dyadc+1))
Dij[is.nan(Dij)] <- 0.
w1 <- rowSums(Dij)
w2 <- Dij%*%w1
l1 <-colSums(Dij)
l2 <- t(l1)%*%Dij
}
if (method == "Pij"){
Pij <- array(dim=c(nrow(X),ncol(X)),0.)
Pij <- X/dyadc
Pij[is.nan(Pij)] <- 0.
w1 <- rowSums(Pij)
w2 <- Pij%*%w1
l1 <-colSums(Pij)
l2 <- t(l1)%*%Pij
}
DS <- w1 + w2 - l1 - t(l2);
if (is.null(names)) names <- paste('Ind.',1:nrow(X))
DS <- array(DS,dim=c(nrow(X),1),dimnames=c(list(names),"David's Scores"))
return(DS)
}
# Function to obtain Landau's linear hierarchy index #
getlandau <- function(X){
V <- freq2abil(X,'matman');
Vvector <- rowSums(V)
sumabil <- sum((Vvector - (length(Vvector) - 1)/2)^2)
landauh <- (12/(nrow(X)^3-nrow(X)))*sumabil;
return(landauh)
}
# Function to obtain improved Landau's linear hierarchy index #
getimplandau <- function(X){
landauh <- getlandau(X);
unknown <- getunknown(X)$unknown.dyads;
improved.landauh <- landauh + (unknown*6/(nrow(X)^3-nrow(X)));
return(improved.landauh)
}
# Function to obtain expected value of linear hierarchy index #
getexpech <- function(X){
expectedh <- 3/(nrow(X)+1);
return(expectedh)
}
# Function to obtain variance for linear hierarchy index #
getvarh <- function(X){
varianceh <- 18/nrow(X)^3;
return(varianceh)
}
# Function to compute the number of circular triads #
getcirc <- function (X){
S <- freq2abil(X,'matman');
Svector <- rowSums(S)
sumabil <- sum(Svector^2)
circular <- length(Svector)*(length(Svector)-1)*(2*length(Svector)-1)/
12-(1/2*sumabil);
return(circular)
}
# Function to compute the expected number of circular triads #
getexpeccirc <- function(X){
expectedcirc <- nrow(X)*(nrow(X)-1)*(nrow(X)-2)/24;
return(expectedcirc)
}
# Function to obtain the maximum number of circular triads #
getmaxcirc <- function(X){
if (nrow(X) %% 2 == 0.) maxcirc <- 1/24*(nrow(X)^3-4*nrow(X))
else maxcirc <- 1/24*(nrow(X)^3-nrow(X));
return(maxcirc)
}
# Function to obtain Kendall's linearity index #
getkendall <- function(X){
circular <- getcirc(X);
maxcirc <- getmaxcirc(X);
K <- 1 - circular/maxcirc;
return(K)
}
# Function to obtain statistical significance for Kendall's linearity index #
getpvkendall <- function(X){
circular <- getcirc(X);
if (nrow(X) >= 10){
degfr <- nrow(X)*(nrow(X)-1)*(nrow(X)-2)/(nrow(X)-4)^2;
chisq <- 8/(nrow(X)-4)*(nrow(X)*(nrow(X)-1)*(nrow(X)-2)/24-circular+1/2)+
degfr;
pval <- 1-pchisq(chisq,degfr)
list(Chisq=chisq,df=degfr,pvalue=pval)
}
else cat("WARNING: ","\n","Matrix size < 10: ","\n",
"Chi square test has not been carried out since it would
overestimate statistical significance","\n")
}
# Function to carry out the two-step randomization test for #
# testing linear hierarchy #
linear.hierarchy.test <- function (X,rep=9999){
# Is the matrix X square? #
if (!squaremat(X))
return("Error: Matrix X is not square and cannot be analyzed")
# Limit the number of replications #
if ((rep < 1) | (rep > 1000000))
return("Error: Number of replications must be between 1 and 1000000")
# Matrix have to be transformed into vectors in order to pass to a .C routine #
vecX <- c(t(X));
# R wrapper of the .C routine #
out <- .C("linearh",
as.double(vecX),
as.integer(nrow(X)),
as.integer(rep),
pvhright=double(1),
pvhleft=double(1),
PACKAGE="DyaDA")
h <- getlandau(X);
himp <- getimplandau(X);
pvright <- out$pvhright;
pvleft <- out$pvhleft;
exph <- getexpech(X);
varh <- getvarh(X);
nind<-nrow(X);
dyads <- nind*(nind-1)/2
total<-sum(X);
circ <- getcirc(X);
expcirc <- getexpeccirc(X);
maxcirc <- getmaxcirc(X);
rho <- getkendall(X);
rho.test <- if (nrow(X)>=10){getpvkendall(X)} else {NULL};
desctable <- rbind(unlist(getunknown(X)),unlist(gettied(X)),
unlist(getonewayrel(X)),unlist(gettwowayrel(X)),
unlist(getsignificant(X)[3:4]))
colnames(desctable) <- c('Frequency','Percentage')
rownames(desctable) <- c('Unknown','Tied','One-way', 'Two-way', 'Significant')
res<-list(nind=nind,dyads=dyads,total=total,desctable=desctable,
randomizations=rep,landauh=h,improved.landauh=himp,
expected.h=exph,variance.h=varh,
right.pvalue=pvright,left.pvalue=pvleft,circ=circ,
expcirc=expcirc,maxcirc=maxcirc,rho=rho,rho.test=rho.test)
class(res) <- "hlineartest"
res
}
print.hlineartest <- function(x,digits=max(4,getOption("digits")-4),...)
{
cat("\n")
cat("Linearity test: Landau's h \n")
cat("========================== \n\n")
cat("\n Number of individuals:", x$nind," (",x$dyads,"dyads)\n")
cat(" Total number of interactions:", x$total,"\n\n")
cat(" Descriptive table of dyadic relationships: \n")
print(x$desctable,digits)
cat("\n")
cat(" Landau's h index:", round(x$landauh,digits),"\n")
cat(" E(h):", round(x$expected.h,digits),"\n")
cat(" Var(h):", round(x$variance.h,digits),"\n")
cat(" Improved Landau's h:", round(x$improved.landauh,digits),"\n")
cat(" Number of randomizations:", x$randomizations,"\n")
cat(" P-value (right):", round(x$right.pvalue,digits),"\n")
cat(" P-value (left):", round(x$left.pvalue,digits),"\n")
cat("\n\n")
cat("Linearity test: Kendall's rho \n")
cat("============================= \n\n")
cat(" Number of circular dyads:", round(x$circ,digits),"\n")
cat(" Expected number of circular dyads:", round(x$expcirc,digits),"\n")
cat(" Maximum number of circular dyads:", round(x$maxcirc,digits),"\n")
cat(" Kendall's rho index:", round(x$rho,digits),"\n")
if (is.null(x$rho.test)){cat("Cannot computed since n<10.")} else{
cat(" Chi-Squared:", round(x$rho.test[[1]],digits),"\n")
cat(" df:", round(x$rho.test[[2]],digits),"\n")
cat(" P-value:", round(x$rho.test[[3]],digits),"\n")
}
invisible(x)
}
#######################################################################################
# I&SI METHOD IN R #
#######################################################################################
DminS.comp <- function(X) {
DminS <- rowSums(sign(X-t(X)))
return(DminS)
}
ordering.matrix <- function(X,order) {
N <- nrow(X);
matord <- X[order,order];
names1 <- dimnames(X)[[1]];
dimnames(matord) <- list(c(names1[order]),c(names1[order]));
return(matord);
}
ISI.comp <- function(X) {
str.inc <- sum((col(X)*(X<t(X))-row(X)*(X<t(X)))[row(X)<col(X)]);
numb.inc <- sum((X<t(X))[row(X)<col(X)]);
list(number.inconsistencies=numb.inc,strength.inconsistencies=str.inc)
}
ISI.info <- function (X,names=NULL){
N <- nrow(X);
I <- sum((X<t(X))[row(X)<col(X)]);
SI <- sum((col(X)*(X<t(X))-row(X)*(X<t(X)))[row(X)<col(X)]);
if (is.null(names)) names <- paste('Ind.',c(1:N))
else names <- dimnames(X)[[1]];
info <- matrix(0,nrow=as.numeric(I),ncol=3,dimnames=list(c(1:as.numeric(I)),
c("Individual i","Individual j","Strength Inconsistency")));
inc.dyads <- matrix(0,nrow=as.numeric(I),ncol=2);
inc.dyads[,1]<-(row(X)*(X<t(X)))[row(X)<col(X)][which((row(X)*(X<t(X)))[row(X)<col(X)]!=0)];
inc.dyads[,2]<-(col(X)*(X<t(X)))[row(X)<col(X)][which((col(X)*(X<t(X)))[row(X)<col(X)]!=0)];
info[,1]<-names[inc.dyads[,1]];
info[,2]<-names[inc.dyads[,2]];
info[,3]<-((col(X)*(X<t(X))-row(X)*(X<t(X)))[row(X)<col(X)])[which((col(X)*(X<t(X))-row(X)*
(X<t(X)))[row(X)<col(X)]!=0)];
info <- data.frame(info)[order(inc.dyads[,1]),];
rownames(info) <- 1:nrow(inc.dyads);
inc.dyads <- inc.dyads[order(inc.dyads[,1]),];
list(inconsistent.dyads=inc.dyads,ISI.frame=info)
}
ISI.method <- function(X,names=NULL,tries=100){
N <- nrow(X)
if (is.null(names)) names <- paste('Ind.',c(1:N)) # IF vector names is NULL rename it #
# Matrices need to be transformed into vectors in order to be passed to a .C routine #
DminS <- rowSums(sign(X-t(X)))
ord <- rank(-DminS, ties.method='random')
Xc <- X[ord,ord]
namesord<- names[ord]
vecX<-c(t(Xc))
out <- .C("ISImethod",as.double(vecX),
as.integer(N),
as.character(namesord),
as.integer(tries),
matord=double(N*N),
namord=character(N),
PACKAGE="DyaDA")
names.ordered <- out$namord
BestXthusfar <- matrix(out$matord,nrow=nrow(X),ncol=ncol(X),byrow=TRUE,
dimnames=list(names.ordered,names.ordered))
res <- list (call=match.call(),dataIni=X,dataFinal=BestXthusfar,initialI=ISI.comp(X)[[1]],
initialSI=ISI.comp(X)[[2]], iniInfo=ISI.info(X)[[2]],finalI=ISI.comp(BestXthusfar)[[1]],
finalSI=ISI.comp(BestXthusfar)[[2]],finalInfo=ISI.info(BestXthusfar)[[2]])
class(res) <- "isi"
res
}
# A print method for isi object #
print.isi <- function(x,digits=max(4,getOption("digits")-4),...)
{
cat("\n")
cat("Ordering a dominance matrix: I&SI Method")
cat("\n\n Call: \n")
cat("",deparse(x$call), "\n\n")
cat(" Initial order: ", rownames(x$dataIni),"\n\n")
print.table(x$dataIni, zero.print = '.')
cat("\n")
cat(" Initial number of inconsistencies: ", x$initialI,"\n")
cat(" Initial strength of incosistencies: ", x$initialI,"\n\n")
cat(" Inconsistent dyads and strength of inconsistencies: ","\n")
print(x$iniInfo,digits=digits)
cat("\n")
cat("**********************************************************","\n")
cat(" Final order: ", rownames(x$dataFinal),"\n\n")
print.table(x$dataFinal, zero.print = '.')
cat("\n")
cat(" Final number of inconsistencies: ", x$finalI,"\n")
cat(" Final strength of inconsistencies: ", x$finalSI,"\n\n")
cat(" Inconsistent dyads and strength of inconsistencies: ","\n")
print(x$finalInfo, digits = digits)
cat("\n\n")
invisible(x)
}
################################################################################
# R FUNCTIONS FOR MATRIX CORRELATION METHODS #
################################################################################
# Function to obtain Mantel's Z statistic #
getZ <- function(X,Y){
# Original sociomatrices should have the same size #
if ((nrow(X) != nrow(Y)) | (ncol(X) != ncol(Y)))
return("Error: Sizes of original sociomatrices are expected to be equal");
XY <- X*Y;
Z <- sum(XY);
return(Z);
}
# Function to compute expected value for Mantel's Z statistic #
getexpecZ <- function(X,Y){
# Original sociomatrices should have the same size #
if ((nrow(X) != nrow(Y)) | (ncol(X) != ncol(Y)))
return("Error: Sizes of original sociomatrices are expected to be equal");
expecZ <- sum(X)*sum(Y)/(nrow(X)*(nrow(X)-1))
return(expecZ);
}
# Function to compute standar error for Mantel's Z statistic #
getSEZ <- function(X,Y){
# Original sociomatrices should have the same size #
if ((nrow(X) != nrow(Y)) | (ncol(X) != ncol(Y)))
return("Error: Sizes of original sociomatrices are expected to be equal");
if (squaremat(X) & squaremat(Y)){
Bx <- sum(X*X);
By <- sum(Y*Y);
Cx <- sum(X*t(X));
Cy <- sum(Y*t(Y));
Dx <- sum(rowSums(X)^2);
Dy <- sum(rowSums(Y)^2);
Ex <- sum(colSums(X)^2);
Ey <- sum(colSums(Y)^2);
Fx <- c(rowSums(X)%*%colSums(X));
Fy <- c(rowSums(Y)%*%colSums(Y));
Gx <- sum(X)**2;
Gy <- sum(Y)**2;
Hx <- Dx - Bx;
Hy <- Dy - By;
Ix <- Ex - Bx;
Iy <- Ey - By;
Jx <- Fx - Cx;
Jy <- Fy - Cy;
Kx <- Gx - Bx - Cx - Hx - Ix - 2*Jx;
Ky <- Gy - By - Cy - Hy - Iy - 2*Jy;
VarZ <- ((Bx*By)+(Cx*Cy)+((Hx*Hy+Ix*Iy+2*Jx*Jy)/(nrow(X)-2))+(Kx*Ky/((nrow(X)-2)*(nrow(X)-3)))-
(Gx*Gy/(nrow(X)*(nrow(X)-1))))/(nrow(X)*(nrow(X)-1));
SEZ <- sqrt(VarZ);
return(SEZ)
}
else
return("Error: Matrices are not square and can not be analyzed")
}
# Function to compute Pearson's correlation for the two sociomatrices #
getpearson <- function(X,Y){
# Original sociomatrices should have the same size #
if ((nrow(X) != nrow(Y)) | (ncol(X) != ncol(Y)))
return("Error: Sizes of original sociomatrices are expected to be equal");
if (squaremat(X) & squaremat(Y)){
X1 <- prepmat(X);
Y1 <- prepmat(Y);
vecX <- c(t(X1));
vecY <- c(t(Y1));
pearson <- cor(vecX,vecY);
}
else {
vecX <- c(t(X));
vecY <- c(t(Y));
pearson <- cor(vecX,vecY);
}
return(pearson)
}
# Function to compute Dietz's R statistic #
getR <- function(X,Y){
# Original sociomatrices should have the same size #
if ((nrow(X) != nrow(Y)) | (ncol(X) != ncol(Y)))
return("Error: Sizes of original sociomatrices are expected to be equal");
if (squaremat(X) & squaremat(Y)){
X1 <- prepmat(X);
Y1 <- prepmat(Y);
Xrank <- rank(X1);
Yrank <- rank(Y1);
XYrank <- Xrank * Yrank;
R <- sum(XYrank);
return(R);}
else {
Xrank <- rank(X);
Yrank <- rank(Y);
XYrank <- Xrank * Yrank;
R <- sum(XYrank);
return(R);}
}
# Function to compute expected value for Dietz's R statistic #
getexpecR <- function(X,Y){
# Original sociomatrices should have the same size #
if ((nrow(X) != nrow(Y)) | (ncol(X) != ncol(Y)))
return("Error: Sizes of original sociomatrices are expected to be equal");
if (squaremat(X) & squaremat(Y)){
X1 <- prepmat(X);
Y1 <- prepmat(Y);
Xrank <- rank(X1);
Yrank <- rank(Y1);
expecR <- sum(Xrank)*sum(Yrank)/(length(Xrank)*(length(Xrank)-1))
return(expecR);}
else {
Xrank <- rank(X);
Yrank <- rank(Y);
expecR <- sum(Xrank)*sum(Yrank)/(length(Xrank)*(length(Xrank)-1))
return(expecR);}
}
# Function to compute standar error for Dietz's R statistic #
getSER <- function(X,Y){
# Original sociomatrices should have the same size #
if ((nrow(X) != nrow(Y)) | (ncol(X) != ncol(Y)))
return("Error: Sizes of original sociomatrices are expected to be equal");
if ((nrow(X) > 3) & (nrow(Y) > 3)){
if (squaremat(X) & squaremat(Y)){
X1 <- prepmat(X);
Y1 <- prepmat(Y);
Xrank <- rank(X1);
Yrank <- rank(Y1);
X[row(X)!=col(X)]<-Xrank;
X<-t(X);
Y[row(Y)!=col(Y)]<-Yrank;
Y<-t(Y);
Bx <- sum(X*X);
By <- sum(Y*Y);
Cx <- sum(X*t(X));
Cy <- sum(Y*t(Y));
Dx <- sum(rowSums(X)^2);
Dy <- sum(rowSums(Y)^2);
Ex <- sum(colSums(X)^2);
Ey <- sum(colSums(Y)^2);
Fx <- sum(rowSums(X)*colSums(X));
Fy <- sum(rowSums(Y)*colSums(Y));
Gx <- sum(X)**2;
Gy <- sum(Y)**2;
Hx <- Dx - Bx;
Hy <- Dy - By;
Ix <- Ex - Bx;
Iy <- Ey - By;
Jx <- Fx - Cx;
Jy <- Fy - Cy;
Kx <- Gx - Bx - Cx - Hx - Ix - 2*Jx;
Ky <- Gy - By - Cy - Hy - Iy - 2*Jy;
VarR <- ((Bx*By)+(Cx*Cy)+((Hx*Hy+Ix*Iy+2*Jx*Jy)/(nrow(X)-2))+(Kx*Ky/((nrow(X)-2)*(nrow(X)-3)))-
(Gx*Gy/(nrow(X)*(nrow(X)-1))))/(nrow(X)*(nrow(X)-1));
SER <- sqrt(VarR);
return(SER)}
else stop("Matrices should be square");
}
else stop("Matrices should have size > 3x3");
}
# Function to compute Spearman's correlation for the two sociomatrices #
getspearman <- function(X,Y){
# Original sociomatrices should have the same size #
if ((nrow(X) != nrow(Y)) | (ncol(X) != ncol(Y)))
return("Error: Sizes of original sociomatrices are expected to be equal");
if (squaremat(X) & squaremat(Y)){
X1 <- prepmat(X);
Y1 <- prepmat(Y);
vecX <- array(dim=c(nrow(X1)*ncol(X1),1));
m <- 0.;
for (i in 1:nrow(X1))
for (j in 1:ncol(X1)){
m <- m+1;
vecX[m] <- X1[i,j];
}
vecY <- array(dim=c(nrow(Y1)*ncol(Y1),1));
m <- 0.;
for (i in 1:nrow(Y1))
for (j in 1:ncol(Y1)){
m <- m+1;
vecY[m] <- Y1[i,j];
}
spearman <- cor(vecX,vecY,method="spearman");
return(spearman)
}
else {
vecX <- array(dim=c(nrow(X)*ncol(X),1));
m <- 0.;
for (i in 1:nrow(X))
for (j in 1:ncol(X)){
m <- m+1;
vecX[m] <- X[i,j];
}
vecY <- array(dim=c(nrow(Y)*ncol(Y),1));
m <- 0.;
for (i in 1:nrow(Y))
for (j in 1:ncol(Y)){
m <- m+1;
vecY[m] <- Y[i,j];
}
spearman <- cor(vecX,vecY,method="spearman");
return(spearman)
}
}
# Function to compute t statistic by means of Mantel's measures #
gett <- function(X,Y,method=c("mantel","dietz")){
# Original sociomatrices should have the same size #
if ((nrow(X) != nrow(Y)) | (ncol(X) != ncol(Y)))
return("Error: Sizes of original sociomatrices are expected to be equal");
method <- match.arg(method)
if (squaremat(X) & squaremat(Y)){
if (method == "mantel"){
Z <- getZ(X,Y);
expecZ <- getexpecZ(X,Y);
SEZ <- getSEZ(X,Y);
t <- (Z-expecZ)/SEZ;
return(t)}
if (method == "dietz"){
R <- getR(X,Y);
expecR <- getexpecR(X,Y);
SER <- getSER(X,Y);
t <- (R-expecR)/SER;
return(t)}
}
else stop("Matrices should be square");
}
# Function to compute rowwise Zr statistic #
getZr <- function(X,Y,names=NULL){
# Original sociomatrices should have the same size #
if ((nrow(X) != nrow(Y)) | (ncol(X) != ncol(Y)))
return("Error: Sizes of original sociomatrices are expected to be equal");
if (is.null(names)) {names <- paste('Ind.',LETTERS[c(1:nrow(X))])}
else {names <- names};
wi <- array(dim=c(nrow(X),1),dimnames=list(c(names),'wi'));
Xterm <- array(dim=c(nrow(X),1));
Yterm <- array(dim=c(nrow(X),1));
if (squaremat(X) & squaremat(Y)){
for (i in 1:nrow(X)){
Xterm[i] <- 0.;
for (j in 1:ncol(X))
for (k in 1:ncol(X)){
if (j > k){
if ((i != j) & (i != k)) Xterm[i] <- Xterm[i] + (X[i,j]-X[i,k])**2
}
}
Yterm[i] <- 0.;
for (j in 1:ncol(X))
for (k in 1:ncol(X)){
if (j < k){
if ((i != j) & (i != k)) Yterm[i] <- Yterm[i] + (Y[i,j]-Y[i,k])**2
}
}
wi[i] <- sqrt(Xterm[i]*Yterm[i]);
}
ri <- array(dim=c(nrow(X),1),dimnames=list(c(names),'ri'));
contrib <- array(dim=c(nrow(X),1),dimnames=list(c(names),'contrib'));
Zr <- 0.;
for (i in 1:nrow(X)){
ri[i] <- 0.;
for (j in 1:ncol(X))
for (k in 1:ncol(X)){
if (j > k){
if ((i != j) & (i != k)) ri[i] <- ri[i] + ((X[i,j]-X[i,k])*(Y[i,j]-Y[i,k])/wi[i])
}
}
Zr <- Zr + wi[i]*ri[i];
contrib[i] <- wi[i]*ri[i];
}
}
else {
for (i in 1:nrow(X)){
Xterm[i] <- 0.;
for (j in 1:ncol(X))
for (k in 1:ncol(X)){
if (j > k){
Xterm[i] <- Xterm[i] + (X[i,j]-X[i,k])**2}
}
Yterm[i] <- 0.;
for (j in 1:ncol(X))
for (k in 1:ncol(X)){
if (j < k){
Yterm[i] <- Yterm[i] + (Y[i,j]-Y[i,k])**2}
}
wi[i] <- sqrt(Xterm[i]*Yterm[i]);
}
ri <- array(dim=c(nrow(X),1),dimnames=list(c(names),'ri'));
contrib <- array(dim=c(nrow(X),1),dimnames=list(c(names),'contrib'));
Zr <- 0.;
for (i in 1:nrow(X)){
ri[i] <- 0.;
for (j in 1:ncol(X))
for (k in 1:ncol(X)){
if (j > k){
ri[i] <- ri[i] + ((X[i,j]-X[i,k])*(Y[i,j]-Y[i,k])/wi[i])}
}
Zr <- Zr + wi[i]*ri[i];
contrib[i] <- wi[i]*ri[i];
}
}
rrwav <- Zr/sum(wi);
rrw <- Zr/(sqrt(sum(Xterm)*sum(Yterm)));
list(Zr=Zr,rrwav=rrwav,rrw=rrw,rowpearson=ri,weighted_contributions=wi,contributions=contrib)
}
# Function to compute rowwise Rr statistic #
getRr <- function(X,Y,names=NULL){
# Original sociomatrices should have the same size #
if ((nrow(X) != nrow(Y)) | (ncol(X) != ncol(Y)))
return("Error: Sizes of original sociomatrices are expected to be equal");
if (is.null(names)) {names <- paste('Ind.',LETTERS[c(1:nrow(X))])}
else {names <- names};
if (squaremat(X) & squaremat(Y)){
X1 <- prepmat(X);
Y1 <- prepmat(Y);
for (i in 1:nrow(X)){
Xrank <- rank(X1[i,]);
Yrank <- rank(Y1[i,])
m <- 1.;
for (j in 1:ncol(X)){
if (i == j) {X[i,j] <- 0.;
Y[i,j] <- 0.}
else {X[i,j] <- Xrank[m];
Y[i,j] <- Yrank[m];
m <- m+1}}
}
wi <- array(dim=c(nrow(X),1),dimnames=list(c(names),'wi'));
Xterm <- array(dim=c(nrow(X),1));
Yterm <- array(dim=c(nrow(X),1));
for (i in 1:nrow(X)){
Xterm[i] <- 0.;
for (j in 1:ncol(X))
for (k in 1:ncol(X)){
if (j > k){
if ((i != j) & (i != k)) Xterm[i] <- Xterm[i] + (X[i,j]-X[i,k])**2
}
}
Yterm[i] <- 0.;
for (j in 1:ncol(X))
for (k in 1:ncol(X)){
if (j < k){
if ((i != j) & (i != k)) Yterm[i] <- Yterm[i] + (Y[i,j]-Y[i,k])**2
}
}
wi[i] <- sqrt(Xterm[i]*Yterm[i]);
}
ri <- array(dim=c(nrow(X),1),dimnames=list(c(names),'ri'));
contrib <- array(dim=c(nrow(X),1),dimnames=list(c(names),'contrib'));
Rr <- 0.;
for (i in 1:nrow(X)){
ri[i] <- 0.;
for (j in 1:ncol(X))
for (k in 1:ncol(X)){
if (j > k){
if ((i != j) & (i != k)) ri[i] <- ri[i] + ((X[i,j]-X[i,k])*(Y[i,j]-Y[i,k])/wi[i])
}
}
Rr <- Rr + wi[i]*ri[i];
contrib[i] <- wi[i]*ri[i];
}
}
else {
for (i in 1:nrow(X)){
Xrank <- rank(X[i,]);
Yrank <- rank(Y[i,])
m <- 1.;
for (j in 1:ncol(X)){
X[i,j] <- Xrank[m];
Y[i,j] <- Yrank[m];
m <- m+1}}
wi <- array(dim=c(nrow(X),1),dimnames=list(c(names),'wi'));
Xterm <- array(dim=c(nrow(X),1));
Yterm <- array(dim=c(nrow(X),1));
for (i in 1:nrow(X)){
Xterm[i] <- 0.;
for (j in 1:ncol(X))
for (k in 1:ncol(X)){
if (j > k){
Xterm[i] <- Xterm[i] + (X[i,j]-X[i,k])**2}
}
Yterm[i] <- 0.;
for (j in 1:ncol(X))
for (k in 1:ncol(X)){
if (j < k){
Yterm[i] <- Yterm[i] + (Y[i,j]-Y[i,k])**2}
}
wi[i] <- sqrt(Xterm[i]*Yterm[i]);
}
ri <- array(dim=c(nrow(X),1),dimnames=list(c(names),'ri'));
contrib <- array(dim=c(nrow(X),1),dimnames=list(c(names),'contrib'));
Rr <- 0.;
for (i in 1:nrow(X)){
ri[i] <- 0.;
for (j in 1:ncol(X))
for (k in 1:ncol(X)){
if (j > k){
ri[i] <- ri[i] + ((X[i,j]-X[i,k])*(Y[i,j]-Y[i,k])/wi[i])}
}
Rr <- Rr + wi[i]*ri[i];
contrib[i] <- wi[i]*ri[i];
}
}
rhorwav <- Rr/sum(wi);
rhorw <- Rr/(sqrt(sum(Xterm)*sum(Yterm)));
list(Rr=Rr,rhorwav=rhorwav,rhorw=rhorw,rowrho=ri,weighted_contributions=wi,contributions=contrib)
}
# Function to compute rowwise Kr statistic #
getKr <- function(X,Y,names=NULL){
# Original sociomatrices should have the same size #
if ((nrow(X) != nrow(Y)) | (ncol(X) != ncol(Y)))
return("Error: Sizes of original sociomatrices are expected to be equal");
if (is.null(names)) {names <- paste('Ind.',LETTERS[c(1:nrow(X))])}
else {names <- names};
wi <- array(dim=c(nrow(X),1),dimnames=list(c(names),'wi'));
Xterm <- array(dim=c(nrow(X),1));
Yterm <- array(dim=c(nrow(X),1));
if (squaremat(X) & squaremat(Y)){
for (i in 1:nrow(X)){
Xterm[i] <- 0.;
for (j in 1:ncol(X))
for (k in 1:ncol(X)){
if (j > k){
if ((i != j) & (i != k)) Xterm[i] <- Xterm[i] + sign(X[i,j]-X[i,k])**2
}
}
Yterm[i] <- 0.;
for (j in 1:ncol(X))
for (k in 1:ncol(X)){
if (j < k){
if ((i != j) & (i != k)) Yterm[i] <- Yterm[i] + sign(Y[i,j]-Y[i,k])**2
}
}
wi[i] <- sqrt(Xterm[i]*Yterm[i]);
}
taui <- array(dim=c(nrow(X),1),dimnames=list(c(names),'taui'));
contrib <- array(dim=c(nrow(X),1),dimnames=list(c(names),'contrib'));
Kr <- 0.;
for (i in 1:nrow(X)){
taui[i] <- 0.;
for (j in 1:ncol(X))
for (k in 1:ncol(X)){
if (j > k){
if ((i != j) & (i != k)) taui[i] <- taui[i] + (sign(X[i,j]-X[i,k])*sign(Y[i,j]-Y[i,k])/wi[i])
}
}
Kr <- Kr + wi[i]*taui[i];
contrib[i] <- wi[i]*taui[i];
}
}
else {
for (i in 1:nrow(X)){
Xterm[i] <- 0.;
for (j in 1:ncol(X))
for (k in 1:ncol(X)){
if (j > k){
Xterm[i] <- Xterm[i] + sign(X[i,j]-X[i,k])**2}
}
Yterm[i] <- 0.;
for (j in 1:ncol(X))
for (k in 1:ncol(X)){
if (j < k){
Yterm[i] <- Yterm[i] + sign(Y[i,j]-Y[i,k])**2}
}
wi[i] <- sqrt(Xterm[i]*Yterm[i]);
}
taui <- array(dim=c(nrow(X),1),dimnames=list(c(names),'taui'));
contrib <- array(dim=c(nrow(X),1),dimnames=list(c(names),'contrib'));
Kr <- 0.;
for (i in 1:nrow(X)){
taui[i] <- 0.;
for (j in 1:ncol(X))
for (k in 1:ncol(X)){
if (j > k){
taui[i] <- taui[i] + (sign(X[i,j]-X[i,k])*sign(Y[i,j]-Y[i,k])/wi[i])}
}
Kr <- Kr + wi[i]*taui[i];
contrib[i] <- wi[i]*taui[i];
}
}
taurwav <- Kr/sum(wi);
taurw <- Kr/(sqrt(sum(Xterm)*sum(Yterm)));
list(Kr=Kr,taurwav=taurwav,taurw=taurw,rowtau=taui,weighted_contributions=wi,contributions=contrib)
}
# Function to compute Kendall's statistics proposed by Dietz #
gettaukendall <- function(X,Y){
# Original sociomatrices should have the same size #
if ((nrow(X) != nrow(Y)) | (ncol(X) != ncol(Y)))
return("Error: Sizes of original sociomatrices are expected to be equal");
if (squaremat(X) & squaremat(Y)){
Kr <- 0.;
for (i in 1:nrow(X))
for (j in 1:ncol(X))
for (k in 1:ncol(X)){
if (j > k){
sig <- sign((X[i,j]-X[i,k])*(Y[i,j]-Y[i,k]));
if ((i != j) & (i != k)) Kr <- Kr + sig
else Kr <- Kr;}}
Kc <- 0.;
for (j in 1:ncol(X))
for (i in 1:nrow(X))
for (k in 1:nrow(X)){
if (i > k ){
sig <- sign((X[i,j]-X[k,j])*(Y[i,j]-Y[k,j]));
if ((i !=j) & (j !=k)) Kc <- Kc + sig
else Kc <- Kc;}}
Kc <- Kc + Kr;
Ku <- 0.;
for (i in 1:nrow(X))
for (j in 1:ncol(X))
for (k in 1:nrow(X))
for (l in 1:ncol(X)){
sig <- sign((X[i,j]-X[k,l])*(Y[i,j]-Y[k,l]));
if ((i < j) & (i < k) & (k < l) & (j != k) & (j != l)) Ku <- Ku + sig
else Ku <- Ku;}
}
else {
Kr <- 0.;
for (i in 1:nrow(X))
for (j in 1:ncol(X))
for (k in 1:ncol(X)){
if (j > k){
sig <- sign((X[i,j]-X[i,k])*(Y[i,j]-Y[i,k]));
Kr <- Kr + sig}
}
Kc <- 0.;
for (j in 1:ncol(X))
for (i in 1:nrow(X))
for (k in 1:nrow(X)){
if (i > k ){
sig <- sign((X[i,j]-X[k,j])*(Y[i,j]-Y[k,j]));
Kc <- Kc + sig}
}
Kc <- Kc + Kr;
Ku <- 0.;
for (i in 1:nrow(X))
for (j in 1:ncol(X))
for (k in 1:nrow(X))
for (l in 1:ncol(X)){
sig <- sign((X[i,j]-X[k,l])*(Y[i,j]-Y[k,l]));
if ((i < j) & (i < k) & (k < l)) Ku <- Ku + sig}
}
K <- Kc + Ku;
list(K=K,Kr=Kr,Kc=Kc,Ku=Ku);
}
# Function to compute partial rowwise Zr statistic #
getpartialZr <- function(X,Y,Z,names=NULL,perm=9999){
if (is.null(names)) {names <- paste('Ind.',LETTERS[c(1:nrow(X))])}
else {names <- names}
# Original sociomatrices should have the same size #
if ((nrow(X) != nrow(Y)) | (ncol(X) != ncol(Y)) | (nrow(X) != nrow(Z)) | (ncol(X) != ncol(Z)))
return("Error: Sizes of original sociomatrices are expected to be equal");
contrib <- array(dim=c(nrow(X),4),dimnames=list(c(names),c('rxy','rxz','ryz','Partial rxy.z')));
rrwxy <- getZr(X,Y)$rrw;
rrwxz <- getZr(X,Z)$rrw;
rrwyz <- getZr(Y,Z)$rrw;
rrwxy.z <- (rrwxy - rrwxz*rrwyz)/sqrt((1-rrwxz**2)*(1-rrwyz**2));
contrib[,1] <- getZr(X,Y)$rowpearson;
contrib[,2] <- getZr(X,Z)$rowpearson;
contrib[,3] <- getZr(Y,Z)$rowpearson;
term <- ((1-contrib[,2]**2)*(1-contrib[,3]**2))**0.5;
for (i in 1:nrow(X)){
if (term[i] == 0.) contrib[i,4] <- 0.
else contrib[i,4] <- (contrib[i,1] - (contrib[i,2]*contrib[i,3]))/term[i];}
# Matrices have to be transformed into vectors in order to pass to a .C routine #
vecX <- array(dim=c(nrow(X)*ncol(X),1));
vecY <- array(dim=c(nrow(Y)*ncol(Y),1));
vecZ <- array(dim=c(nrow(Z)*ncol(Z),1));
m <- 0.;
for (i in 1:nrow(X))
for (j in 1:ncol(X)){
m <- m+1;
vecX[m] <- X[i,j];
vecY[m] <- Y[i,j];
vecZ[m] <- Z[i,j];
}
if (squaremat(X) & squaremat(Y) & squaremat(Z)){
# R wrapper of the .C routine #
if (nrow(X) < 9){
out <- .C("partialexactsigZr",
as.double(vecX),
as.double(vecY),
as.double(vecZ),
as.integer(nrow(X)),
pvpartialZrright=double(1),
pvpartialZrleft=double(1),
PACKAGE="DyaDA")}
else {
out <- .C("partialstatsigZr",
as.double(vecX),
as.double(vecY),
as.double(vecZ),
as.integer(nrow(X)),
as.integer(perm),
pvpartialZrright=double(1),
pvpartialZrleft=double(1),
PACKAGE="DyaDA")}
}
else {
out <- .C("rectangularstatsigpartZr",
as.double(vecX),
as.double(vecY),
as.double(vecZ),
as.integer(nrow(X)),
as.integer(ncol(X)),
as.integer(perm),
pvpartialZrright=double(1),
pvpartialZrleft=double(1),
PACKAGE="DyaDA")
}
pvright <- out$pvpartialZrright;
pvleft <- out$pvpartialZrleft;
list(rrwxy.z=rrwxy.z,p_value_right=pvright,p_value_left=pvleft,rrwxy=rrwxy,rrwxz=rrwxz,rrwyz=rrwyz,partial_rowwise_Z=contrib);
}
# Function to compute partial rowwise Rr statistic #
getpartialRr <- function(X,Y,Z,names=NULL,perm=9999){
if (is.null(names)) {names <- paste('Ind.',LETTERS[c(1:nrow(X))])}
else {names <- names}
# Original sociomatrices should have the same size #
if ((nrow(X) != nrow(Y)) | (ncol(X) != ncol(Y)) | (nrow(X) != nrow(Z)) | (ncol(X) != ncol(Z)))
return("Error: Sizes of original sociomatrices are expected to be equal");
contrib <- array(dim=c(nrow(X),4),dimnames=list(c(names),c('rhoxy','rhoxz','rhoyz','Partial rhoxy.z')));
rhorwxy <- getRr(X,Y)$rhorw;
rhorwxz <- getRr(X,Z)$rhorw;
rhorwyz <- getRr(Y,Z)$rhorw;
rhorwxy.z <- (rhorwxy - rhorwxz*rhorwyz)/sqrt((1-rhorwxz**2)*(1-rhorwyz**2));
contrib[,1] <- getRr(X,Y)$rowrho;
contrib[,2] <- getRr(X,Z)$rowrho;
contrib[,3] <- getRr(Y,Z)$rowrho;
term <- ((1-contrib[,2]**2)*(1-contrib[,3]**2))**0.5;
for (i in 1:nrow(X)){
if (term[i] == 0.) contrib[i,4] <- 0.
else contrib[i,4] <- (contrib[i,1] - (contrib[i,2]*contrib[i,3]))/term[i];}
if (squaremat(X) & squaremat(Y) & squaremat(Z)){
# Matrices have to be transformed into vectors in order to pass to a .C routine #
X1 <- prepmat(X);
Y1 <- prepmat(Y);
Z1 <- prepmat(Z);
for (i in 1:nrow(X)){
Xrank <- rank(X1[i,]);
Yrank <- rank(Y1[i,]);
Zrank <- rank(Z1[i,]);
m <- 1.;
for (j in 1:ncol(X)){
if (i == j) {X[i,j] <- 0.;
Y[i,j] <- 0.;
Z[i,j] <- 0.;}
else {X[i,j] <- Xrank[m];
Y[i,j] <- Yrank[m];
Z[i,j] <- Zrank[m];
m <- m+1}}
}
vecX <- array(dim=c(nrow(X)*ncol(X),1));
vecY <- array(dim=c(nrow(Y)*ncol(Y),1));
vecZ <- array(dim=c(nrow(Z)*ncol(Z),1));
m <- 0.;
for (i in 1:nrow(X))
for (j in 1:ncol(X)){
m <- m+1;
vecX[m] <- X[i,j];
vecY[m] <- Y[i,j];
vecZ[m] <- Z[i,j];
}
# R wrapper of the .C routine #
if (nrow(X) < 9){
out <- .C("partialexactsigRr",
as.double(vecX),
as.double(vecY),
as.double(vecZ),
as.integer(nrow(X)),
pvpartialRrright=double(1),
pvpartialRrleft=double(1),
PACKAGE=DYaDA)}
else {
out <- .C("partialstatsigRr",
as.double(vecX),
as.double(vecY),
as.double(vecZ),
as.integer(nrow(X)),
as.integer(perm),
pvpartialRrright=double(1),
pvpartialRrleft=double(1),
PACKAGE=DYaDA)}
}
else {
for (i in 1:nrow(X)){
Xrank <- rank(X[i,]);
Yrank <- rank(Y[i,]);
Zrank <- rank(Z[i,]);
m <- 1.;
for (j in 1:ncol(X)){
X[i,j] <- Xrank[m];
Y[i,j] <- Yrank[m];
Z[i,j] <- Zrank[m];
m <- m+1}}
vecX <- array(dim=c(nrow(X)*ncol(X),1));
vecY <- array(dim=c(nrow(Y)*ncol(Y),1));
vecZ <- array(dim=c(nrow(Z)*ncol(Z),1));
m <- 0.;
for (i in 1:nrow(X))
for (j in 1:ncol(X)){
m <- m+1;
vecX[m] <- X[i,j];
vecY[m] <- Y[i,j];
vecZ[m] <- Z[i,j];
}
out <- .C("rectangularstatsigpartRr",as.double(vecX),
as.double(vecY),
as.double(vecZ),
as.integer(nrow(X)),
as.integer(ncol(X)),
as.integer(perm),
pvpartialRrright=double(1),
pvpartialRrleft=double(1),
PACKAGE=DyaDA)
}
pvright <- out$pvpartialRrright;
pvleft <- out$pvpartialRrleft;
list(rhorwxy.z=rhorwxy.z,p_value_right=pvright,p_value_left=pvleft,rhorwxy=rhorwxy,rhorwxz=rhorwxz,rhorwyz=rhorwyz,partial_rowwise_R=contrib);
}
# Function to compute partial rowwise Kr statistic #
getpartialKr <- function(X,Y,Z,names=NULL,perm=9999){
if (is.null(names)) {names <- paste('Ind.',LETTERS[c(1:nrow(X))])}
else {names <- names}
# Original sociomatrices should have the same size #
if ((nrow(X) != nrow(Y)) | (ncol(X) != ncol(Y)) | (nrow(X) != nrow(Z)) | (ncol(X) != ncol(Z)))
return("Error: Sizes of original sociomatrices are expected to be equal");
contrib <- array(dim=c(nrow(X),4),dimnames=list(c(names),c('tauxy','tauxz','tauyz','Partial tauxy.z')));
taurwxy <- getKr(X,Y)$taurw;
taurwxz <- getKr(X,Z)$taurw;
taurwyz <- getKr(Y,Z)$taurw;
taurwxy.z <- (taurwxy - taurwxz*taurwyz)/sqrt((1-taurwxz**2)*(1-taurwyz**2));
contrib[,1] <- getKr(X,Y)$rowtau;
contrib[,2] <- getKr(X,Z)$rowtau;
contrib[,3] <- getKr(Y,Z)$rowtau;
term <- ((1-contrib[,2]**2)*(1-contrib[,3]**2))**0.5;
for (i in 1:nrow(X)){
if (term[i] == 0.) contrib[i,4] <- 0.
else contrib[i,4] <- (contrib[i,1] - (contrib[i,2]*contrib[i,3]))/term[i];}
# Matrices have to be transformed into vectors in order to pass to a .C routine #
vecX <- array(dim=c(nrow(X)*ncol(X),1));
vecY <- array(dim=c(nrow(Y)*ncol(Y),1));
vecZ <- array(dim=c(nrow(Z)*ncol(Z),1));
m <- 0.;
for (i in 1:nrow(X))
for (j in 1:ncol(X)){
m <- m+1;
vecX[m] <- X[i,j];
vecY[m] <- Y[i,j];
vecZ[m] <- Z[i,j];
}
if (squaremat(X) & squaremat(Y) & squaremat(Z)){
# R wrapper of the .C routine #
if (nrow(X) < 9){
out <- .C("partialexactsiKr",
as.double(vecX),
as.double(vecY),
as.double(vecZ),
as.integer(nrow(X)),
pvpartialKrright=double(1),
pvpartialKrleft=double(1),
PACKAGE="DyaDA")}
else {
out <- .C("partialstatsigKr",
as.double(vecX),
as.double(vecY),
as.double(vecZ),
as.integer(nrow(X)),
as.integer(perm),
pvpartialKrright=double(1),
pvpartialKrleft=double(1),
PACKAGE="DyaDA")}
}
else {
out <- .C("rectangularstatsigpartKr",
as.double(vecX),
as.double(vecY),
as.double(vecZ),
as.integer(nrow(X)),
as.integer(ncol(X)),
as.integer(perm),
pvpartialKrright=double(1),
pvpartialKrleft=double(1),
PACKAGE="DyaDA")
}
pvright <- out$pvpartialKrright;
pvleft <- out$pvpartialKrleft;
list(taurwxy.z=taurwxy.z,p_value_right=pvright,p_value_left=pvleft,taurwxy=taurwxy,taurwxz=taurwxz,taurwyz=taurwyz,partial_rowwise_Kr=contrib);
}
# Function to carry out the permutation test to estimate statistical significance for Mantel's Z statistic #
mantelZtest <- function (X,Y,perm=9999){
# Limit the number of replications #
if ((perm < 1) | (perm > 1000000))
return("Error: Number of permutations must be between 1 and 1000000");
# Original sociomatrices should have the same size #
if ((nrow(X) != nrow(Y)) | (ncol(X) != ncol(Y)))
return("Error: Sizes of original sociomatrices are expected to be equal");
# Matrices have to be transformed into vectors in order to pass to a .C routine #
vecX <- array(dim=c(nrow(X)*ncol(X),1));
m <- 0.;
for (i in 1:nrow(X))
for (j in 1:ncol(X))
{
m <- m+1;
vecX[m] <- X[i,j];
}
vecY <- array(dim=c(nrow(Y)*ncol(Y),1));
m <- 0.;
for (i in 1:nrow(Y))
for (j in 1:ncol(Y))
{
m <- m+1;
vecY[m] <- Y[i,j];
}
if (squaremat(X) & squaremat(Y)){
# R wrapper of the .C routine #
if (nrow(X) < 9){
out <- .C("exactsigZ",as.double(vecX),
as.double(vecY),
as.integer(nrow(X)),
pvZright=double(1),
pvZleft=double(1),
PACKAGE="DyaDA")}
else {
out <- .C("statsigZ",as.double(vecX),
as.double(vecY),
as.integer(nrow(X)),
as.integer(perm),
pvZright=double(1),
pvZleft=double(1),
PACKAGE="DyaDA")}
}
else {
out <- .C("rectangularstatsigZ",as.double(vecX),
as.double(vecY),
as.integer(nrow(X)),
as.integer(ncol(X)),
as.integer(perm),
pvZright=double(1),
pvZleft=double(1),
PACKAGE="DyaDA")
}
z <- getZ(X,Y);
pvright <- out$pvZright;
pvleft <- out$pvZleft;
list(Mantel_Z=z,right_pvalue=pvright,left_pvalue=pvleft)
}
# Function to carry out the permutation test to estimate statistical significance for Dietz's R statistic #
dietzRtest <- function (X,Y,perm=9999){
# Limit the number of replications #
if ((perm < 1) | (perm > 1000000))
return("Error: Number of permutations must be between 1 and 1000000");
# Original sociomatrices should have the same size #
if ((nrow(X) != nrow(Y)) | (ncol(X) != ncol(Y)))
return("Error: Sizes of original sociomatrices are expected to be equal");
if (squaremat(X) & squaremat(Y)){
# Matrices have to be transformed into vectors in order to pass to a .C routine #
X1 <- prepmat(X);
Y1 <- prepmat(Y);
Xrank <- rank(X1);
Yrank <- rank(Y1);
m <- 1.;
for (i in 1:nrow(X))
for (j in 1:ncol(X)){
if (i == j) {X[i,j] <- 0.}
else {X[i,j] <- Xrank[m];
m <- m+1;}
}
m <- 1.;
for (i in 1:nrow(Y))
for (j in 1:ncol(Y)){
if (i == j) {Y[i,j] <- 0.}
else {Y[i,j] <- Yrank[m];
m <- m+1;}
}
vecX <- array(dim=c(nrow(X)*ncol(X),1));
m <- 0.;
for (i in 1:nrow(X))
for (j in 1:ncol(X)){
m <- m+1;
vecX[m] <- X[i,j];
}
vecY <- array(dim=c(nrow(Y)*ncol(Y),1));
m <- 0.;
for (i in 1:nrow(Y))
for (j in 1:ncol(Y)){
m <- m+1;
vecY[m] <- Y[i,j];
}
# R wrapper of the .C routine #
if (nrow(X) < 9){
out <- .C("exactsigR",as.double(vecX),
as.double(vecY),
as.integer(nrow(X)),
pvRright=double(1),
pvRleft=double(1),
PACKAGE="DyaDA")}
else {
out <- .C("statsigR",as.double(vecX),
as.double(vecY),
as.integer(nrow(X)),
as.integer(perm),
pvRright=double(1),
pvRleft=double(1),
PACKAGE="DyaDA")}
}
else {
Xrank <- rank(X);
Yrank <- rank(Y);
m <- 1.;
for (i in 1:nrow(X))
for (j in 1:ncol(X)){
X[i,j] <- Xrank[m];
m <- m+1;}
m <- 1.;
for (i in 1:nrow(Y))
for (j in 1:ncol(Y)){
Y[i,j] <- Yrank[m];
m <- m+1;}
vecX <- array(dim=c(nrow(X)*ncol(X),1));
m <- 0.;
for (i in 1:nrow(X))
for (j in 1:ncol(X)){
m <- m+1;
vecX[m] <- X[i,j];
}
vecY <- array(dim=c(nrow(Y)*ncol(Y),1));
m <- 0.;
for (i in 1:nrow(Y))
for (j in 1:ncol(Y)){
m <- m+1;
vecY[m] <- Y[i,j];
}
out <- .C("rectangularstatsigR",as.double(vecX),
as.double(vecY),
as.integer(nrow(X)),
as.integer(ncol(X)),
as.integer(perm),
pvRright=double(1),
pvRleft=double(1),
PACKAGE="DyaDA")
}
R <- getR(X,Y);
pvright <- out$pvRright;
pvleft <- out$pvRleft;
list(Dietz_R=R,right_pvalue=pvright,left_pvalue=pvleft)
}
# Function to carry out the permutation test to estimate statistical significance for Mantel's Zr statistic #
rowwiseZrtest <- function (X,Y,perm=9999){
# Limit the number of replications #
if ((perm < 1) | (perm > 1000000))
return("Error: Number of permutations must be between 1 and 1000000");
# Original sociomatrices should have the same size #
if ((nrow(X) != nrow(Y)) | (ncol(X) != ncol(Y)))
return("Error: Sizes of original sociomatrices are expected to be equal");
# Matrices have to be transformed into vectors in order to pass to a .C routine #
vecX <- array(dim=c(nrow(X)*ncol(X),1));
m <- 0.;
for (i in 1:nrow(X))
for (j in 1:ncol(X)){
m <- m+1;
vecX[m] <- X[i,j];
}
vecY <- array(dim=c(nrow(Y)*ncol(Y),1));
m <- 0.;
for (i in 1:nrow(Y))
for (j in 1:ncol(Y)){
m <- m+1;
vecY[m] <- Y[i,j];
}
if (squaremat(X) & squaremat(Y)){
# R wrapper of the .C routine #
if (nrow(X) < 9){
out <- .C("exactsigZr",as.double(vecX),
as.double(vecY),
as.integer(nrow(X)),
pvZrright=double(1),
pvZrleft=double(1),
PACKAGE="DyaDA")}
else {
out <- .C("statsigZr",as.double(vecX),
as.double(vecY),
as.integer(nrow(X)),
as.integer(perm),
pvZrright=double(1),
pvZrleft=double(1),
PACKAGE="DyaDA")}
}
else {
out <- .C("rectangularstatsigZr",as.double(vecX),
as.double(vecY),
as.integer(nrow(X)),
as.integer(ncol(X)),
as.integer(perm),
pvZrright=double(1),
pvZrleft=double(1),
PACKAGE="DyaDA")
}
Zr <- getZr(X,Y)$Zr;
pvright <- out$pvZrright;
pvleft <- out$pvZrleft;
list(Zr=Zr,right_pvalue=pvright,left_pvalue=pvleft)
}
# Function to carry out the permutation test to estimate statistical significance for Dietz's Rr statistic #
rowwiseRrtest <- function (X,Y,perm=9999){
# Limit the number of replications #
if ((perm < 1) | (perm > 1000000))
return("Error: Number of permutations must be between 1 and 1000000");
# Original sociomatrices should have the same size #
if ((nrow(X) != nrow(Y)) | (ncol(X) != ncol(Y)))
return("Error: Sizes of original sociomatrices are expected to be equal");
if (squaremat(X) & squaremat(Y)){
# Matrices have to be transformed into vectors in order to pass to a .C routine #
X1 <- prepmat(X);
Y1 <- prepmat(Y);
for (i in 1:nrow(X)){
Xrank <- rank(X1[i,]);
Yrank <- rank(Y1[i,]);
m <- 1.;
for (j in 1:ncol(X)){
if (i == j) {X[i,j] <- 0.;
Y[i,j] <- 0.}
else {X[i,j] <- Xrank[m];
Y[i,j] <- Yrank[m];
m <- m+1}}
}
vecX <- array(dim=c(nrow(X)*ncol(X),1));
m <- 0.;
for (i in 1:nrow(X))
for (j in 1:ncol(X))
{
m <- m+1;
vecX[m] <- X[i,j];
}
vecY <- array(dim=c(nrow(Y)*ncol(Y),1));
m <- 0.;
for (i in 1:nrow(Y))
for (j in 1:ncol(Y))
{
m <- m+1;
vecY[m] <- Y[i,j];
}
# R wrapper of the .C routine #
if (nrow(X) < 9){
out <- .C("exactsigRr",as.double(vecX),
as.double(vecY),
as.integer(nrow(X)),
pvRrright=double(1),
pvRrleft=double(1),
PACKAGE="DyaDA")}
else {
out <- .C("statsigRr",as.double(vecX),
as.double(vecY),
as.integer(nrow(X)),
as.integer(perm),
pvRrright=double(1),
pvRrleft=double(1),
PACKAGE="DyaDA")}
}
else {
for (i in 1:nrow(X)){
Xrank <- rank(X[i,]);
Yrank <- rank(Y[i,])
m <- 1.;
for (j in 1:ncol(X)){
X[i,j] <- Xrank[m];
Y[i,j] <- Yrank[m];
m <- m+1}}
vecX <- array(dim=c(nrow(X)*ncol(X),1));
m <- 0.;
for (i in 1:nrow(X))
for (j in 1:ncol(X)){
m <- m+1;
vecX[m] <- X[i,j];
}
vecY <- array(dim=c(nrow(Y)*ncol(Y),1));
m <- 0.;
for (i in 1:nrow(Y))
for (j in 1:ncol(Y)){
m <- m+1;
vecY[m] <- Y[i,j];
}
out <- .C("rectangularstatsigRr",
as.double(vecX),
as.double(vecY),
as.integer(nrow(X)),
as.integer(ncol(X)),
as.integer(perm),
pvRrright=double(1),
pvRrleft=double(1),
PACKAGE="DyaDA")
}
# R wrapper of the .C routine #
Rr <- getRr(X,Y)$Rr;
pvright <- out$pvRrright;
pvleft <- out$pvRrleft;
list(Rr=Rr,right_pvalue=pvright,left_pvalue=pvleft)
}
# Function to carry out the permutation test to estimate statistical significance for Kendall's Kr statistic #
rowwiseKrtest <- function (X,Y,perm=9999){
# Limit the number of replications #
if ((perm < 1) | (perm > 1000000))
return("Error: Number of permutations must be between 1 and 1000000");
# Original sociomatrices should have the same size #
if ((nrow(X) != nrow(Y)) | (ncol(X) != ncol(Y)))
return("Error: Sizes of original sociomatrices are expected to be equal");
# Matrices have to be transformed into vectors in order to pass to a .C routine #
vecX <- array(dim=c(nrow(X)*(ncol(X)-1)/2,1));
m <- 0.;
for (i in 1:nrow(X))
for (j in 1:ncol(X))
{
m <- m+1;
vecX[m] <- X[i,j];
}
vecY <- array(dim=c(nrow(Y)*(ncol(Y)-1)/2,1));
m <- 0.;
for (i in 1:nrow(Y))
for (j in 1:ncol(Y))
{
m <- m+1;
vecY[m] <- Y[i,j];
}
if (squaremat(X) & squaremat(Y)){
# R wrapper of the .C routine #
if (nrow(X) < 9){
out <- .C("exactsigKr",as.double(vecX),
as.double(vecY),
as.integer(nrow(X)),
pvKrright=double(1),
pvKrleft=double(1),
PACKAGE="DyaDA")}
else {
out <- .C("statsigKr",as.double(vecX),
as.double(vecY),
as.integer(nrow(X)),
as.integer(perm),
pvKrright=double(1),
pvKrleft=double(1),
PACKAGE="DyaDA")}
}
else {
out <- .C("rectangularstatsigKr",as.double(vecX),
as.double(vecY),
as.integer(nrow(X)),
as.integer(ncol(X)),
as.integer(perm),
pvKrright=double(1),
pvKrleft=double(1),
PACKAGE="DyaDA")
}
Kr <- getKr(X,Y)$Kr;
pvright <- out$pvKrright;
pvleft <- out$pvKrleft;
list(kendall_Kr=Kr,right_pvalue=pvright,left_pvalue=pvleft)
}
################################################################################
# R FUNCTIONS FOR DYADIC INTERDEPENDENCE #
################################################################################
# A function for estimating dyadic interdependence in standard dyadic designs with distinguishable dyad members and interval responses #
distinguishable.dyad <- function(dataset,conf.int=0.95){
N<-length(na.omit(dataset)[,2])
na1<- sum(is.na(dataset[,2]))
na2<- sum(is.na(dataset[,3]))
dataset<-na.omit(dataset)
summary.statistics <- array(c(N,na1,min(dataset[,2]),max(dataset[,2]),mean(dataset[,2]),
sd(dataset[,2]),quantile(dataset[,2],.25,names=F),quantile(dataset[,2],.50,names=F),
quantile(dataset[,2],.75,names=F),N,na2,min(dataset[,3]),max(dataset[,3]),
mean(dataset[,3]),sd(dataset[,3]),quantile(dataset[,3],.25,names=F),quantile(dataset[,3],.50,names=F),
quantile(dataset[,3],.75,names=F)),dim=c(9,2),dimnames=list(c("N:", "NAs:","Min:","Max:","Mean:","Sd:",
"25th Pctl:","50th Pctl:", "75th Pctl:"),c(colnames(dataset[,2:3]))))
r <- cor(dataset[,2],dataset[,3])
df <- N-2
alpha <- 1 - conf.int
t.statistic <- r*sqrt(N-2)/(sqrt(1-r**2))
p.value <- 1 - pt(t.statistic,df)
zr.value <- 0.5*log((1+r)/(1-r))
zr.low <- zr.value - qnorm((1-alpha/2),0,1)/sqrt(N-3)
zr.upper <- zr.value + qnorm((1-alpha/2),0,1)/sqrt(N-3)
r.low <- (exp(1)**(2*zr.low)-1)/(exp(1)**(2*zr.low)+1)
r.upper <- (exp(1)**(2*zr.upper)-1)/(exp(1)**(2*zr.upper)+1)
res <- list(call=match.call(),data=dataset,stats=summary.statistics,r=r,rCI=c(r.low,r.upper),
t=t.statistic,df=df,pval=2*p.value)
class(res) <- 'dyadr'
res
}
print.dyadr <- function(x,digits=max(4,getOption("digits")-4),...)
{
cat("\n")
cat("Distinguishable Members: Pearson's Correlation")
cat("\n\n Call: \n")
cat("",deparse(x$call), "\n\n")
cat(" Data", if (length(x$data[,1]) <= 5) ": " else " (5 first rows shown): ", "\n")
print( if (length(x$data[,1]) >= 5) x$data[1:5,] else x$data[1:length(x$data),])
cat("\n\n")
cat("Descriptive Statistics","\n")
print.table(x$stats, digits=digits)
cat("\n\n")
cat("t Test: ","\n")
results <- cbind(x$t,as.integer(x$df),x$pval)
colnames(results) <- c("t Value", "Df", "Pr(>|t|)")
rownames(results) <- ""
print.table(results,digits=digits)
cat("\n\n")
cat(paste((1-x$alpha)*100,"% Confidence Interval: ",sep=''),"\n")
results <- t(x$rCI)
colnames(results) <- c("Lower","Upper")
rownames(results) <- ''
print.table(results,digits=digits)
cat("\n\n")
cat("Sample estimate","\n")
results <- cbind(x$r)
colnames(results) <- "cor"
rownames(results) <- ''
print.table(results,digits=digits)
cat("\n\n")
invisible(x)
}
# A function for estimating dyadic interdependence in standard dyadic designs with indistinguishable dyad members and interval responses #
indistinguishable.dyad <- function(dataset,conf.int=0.95){
N<-length(na.omit(dataset)[,2])
na1<- sum(is.na(dataset[,2]))
na2<- sum(is.na(dataset[,3]))
dataset<-na.omit(dataset)
summary.statistics <- array(c(N,na1,min(dataset[,2]),max(dataset[,2]),mean(dataset[,2]),
sd(dataset[,2]),quantile(dataset[,2],.25,names=F),quantile(dataset[,2],.50,names=F),
quantile(dataset[,2],.75,names=F),N,na2,min(dataset[,3]),max(dataset[,3]),
mean(dataset[,3]),sd(dataset[,3]),quantile(dataset[,3],.25,names=F),quantile(dataset[,3],.50,names=F),
quantile(dataset[,3],.75,names=F)),dim=c(9,2),dimnames=list(c("N:", "NAs:","Min:","Max:","Mean:","Sd:",
"25th Pctl:","50th Pctl:", "75th Pctl:"),c(colnames(dataset[,2:3]))))
MSb <- 2*var(apply(dataset[,2:3],1,mean))
MSw <- sum((dataset[,2]-dataset[,3])**2)/(2*N)
ICC <- (MSb-MSw)/(MSb+MSw)
alpha <- 1 - conf.int
if ( MSb > MSw ) {
F.statistic <- MSb/MSw
df1 <- N-1
df2 <- N}
if (MSw > MSb ){
F.statistic <- MSw/MSb
df1 <- N
df2 <- N-1}
p.value <- pf(F.statistic,df1,df2,lower.tail=FALSE)
fl <- MSb/MSw/qf(1-alpha/2,df1,df2)
fu <- MSb/MSw*qf(1-alpha/2,df1,df2)
icc.low <- (fl-1)/(fl+1)
icc.upper <- (fu-1)/(fu+1)
res=list(call=match.call(),data=dataset,stats=summary.statistics,
intracor=ICC,MSb=MSb,MSw=MSw,Fstat=F.statistic,df1=df1,df2=df2,pval=2*p.value,
alpha=alpha,iccCI=c(icc.low,icc.upper))
class(res)="dyadicc"
res
}
print.dyadicc <- function(x,digits=max(4,getOption("digits")-4),...)
{
cat("\n")
cat("Indistinguishable Members: Intraclass Correlation")
cat("\n\n Call: \n")
cat("",deparse(x$call), "\n\n")
cat(" Data", if (length(x$data[,1]) <= 5) ": " else " (5 first rows shown): ", "\n")
print( if (length(x$data[,1]) >= 5) x$data[1:5,] else x$data[1:length(x$data),])
cat("\n\n")
cat("Descriptive Statistics","\n")
print.table(x$stats, digits=digits)
cat("\n\n")
cat("F Test: ","\n")
results <- cbind(rbind(x$MSb,x$MSw),rbind(x$df1,x$df2),rbind(x$Fstat,NA),rbind(x$pval,NA))
colnames(results) <- c("Mean Sq", "Df", "F value", "Pr(>|F|)")
rownames(results) <- c("Between-Dyads","Within-Dyads")
print.table(results,digits=digits)
cat("\n\n")
cat(paste("Intraclass Correlation and ",(1-x$alpha)*100,"% Confidence Interval: ",sep=''),"\n")
results <- cbind(x$intracor,t(x$iccCI))
colnames(results) <- c("ICC", "Lower","Upper")
rownames(results) <- ''
print.table(results,digits=digits)
cat("\n\n")
invisible(x)
}
# A function for estimating dyadic interdependence in standard dyadic designs with indistinguishable dyad members and interval responses by means of the pairwise correlation #
pairwise.correlation <-
function (dataset, conf.int = 0.95){
alpha <- 1 - conf.int
double.entry.dataset <- array(c(c(1:(2 * length(dataset[,1]))), c(dataset[, 2], dataset[, 3]),
c(dataset[, 3],dataset[, 2])), dim = c((2 * length(dataset[, 1])), length(dataset)))
colnames(double.entry.dataset) <- colnames(dataset)
dataset <- double.entry.dataset
N <- length(na.omit(dataset)[, 2])
na1 <- sum(is.na(dataset[, 2]))
na2 <- sum(is.na(dataset[, 3]))
dataset <- na.omit(dataset)
summary.statistics <- array(c(N, na1, min(dataset[, 2]),
max(dataset[, 2]), mean(dataset[, 2]), sd(dataset[, 2]),
quantile(dataset[, 2], 0.25, names = F), quantile(dataset[,2],0.5, names = F),
quantile(dataset[, 2], 0.75,names = F), N, na2, min(dataset[, 3]), max(dataset[,
3]), mean(dataset[, 3]), sd(dataset[, 3]), quantile(dataset[,
3], 0.25, names = F), quantile(dataset[, 3], 0.5,
names = F), quantile(dataset[, 3], 0.75, names = F)),
dim = c(9, 2), dimnames = list(c("N:", "NAs:", "Min:",
"Max:", "Mean:", "Sd:", "25th Pctl:", "50th Pctl:",
"75th Pctl:"), c(colnames(dataset[, 2:3]))))
rp <- cor(dataset[, 2], dataset[, 3])
standard.error <- 1/sqrt(N/2)
z <- rp/standard.error
p.value <- 1 - pnorm(abs(z), 0, 1)
rp.low <- rp - (qnorm((1 - alpha/2), 0, 1)/sqrt(N/2))
rp.upper <- rp + (qnorm((1 - alpha/2), 0, 1)/sqrt(N/2))
res <- list(call = match.call(), data = dataset, stats = summary.statistics,
rp = rp, alpha = alpha, rpCI = c(rp.low, rp.upper), zVal = z,
stdErr = standard.error, pval = 2 * p.value)
class(res) <- "rp"
res
}
print.rp <- function(x,digits=max(4,getOption("digits")-4),...)
{
cat("\n")
cat("Indistinguishable Members: Pairwise Correlation")
cat("\n\n Call: \n")
cat("",deparse(x$call), "\n\n")
cat(" Data", if (length(x$data[,1]) <= 5) ": " else " (5 first rows shown): ", "\n")
print( if (length(x$data[,1]) >= 5) x$data[1:5,] else x$data[1:length(x$data),])
cat("\n\n")
cat("Descriptive Statistics","\n")
print.table(x$stats, digits=digits)
cat("\n\n")
cat("z Test: ","\n")
results <- cbind(x$zVal,x$stdErr,x$pval)
colnames(results) <- c("Z Value", "Std Error", "Pr(>|Z|)")
rownames(results) <- ""
print.table(results,digits=digits)
cat("\n\n")
cat(paste((1-x$alpha)*100,"% Confidence Interval: ",sep=''),"\n")
results <- t(x$rpCI)
colnames(results) <- c("Lower","Upper")
rownames(results) <- ''
print.table(results,digits=digits)
cat("\n\n")
cat("Sample estimate","\n")
results <- cbind(x$rp)
colnames(results) <- "rp"
rownames(results) <- ''
print.table(results,digits=digits)
cat("\n\n")
invisible(x)
}
# A function for estimating dyadic interdependence in standard dyadic designs with distinguishable dyad members and categorical responses #
categorical.distinguishable.dyad <- function(dataset,conf.int=0.95){
tabulated.data <- table(dataset[,2],dataset[,3])
alpha <- 1 - conf.int
no <- sum(diag(tabulated.data))
ne <- sum(margin.table(tabulated.data,1)*margin.table(tabulated.data,2)/sum(tabulated.data))
kappa <- (no - ne)/(sum(tabulated.data)-ne)
table.props <- prop.table(tabulated.data)
prop.margins <-array(c(margin.table(table.props,1),margin.table(prop.table(tabulated.data),2)),dim=c(nrow(tabulated.data),2))
sumterm1 <- 0.
sumterm2 <- 0.
sumterm3 <- 0.
pe <- ne/sum(tabulated.data)
for (i in 1:nrow(table.props)){
sumterm1 <- sumterm1+(table.props[i,i]*(1-sum(prop.margins[i,])*(1-kappa))**2+(1-kappa)**2)
sumterm3 <- sumterm3+(prop.margins[i,1]*prop.margins[i,2]*sum(prop.margins[i,]))}
for (i in 1:nrow(table.props)){
for (j in 1:ncol(table.props)){
sumterm2 <- sumterm2+(table.props[i,j]*(prop.margins[i,2]+prop.margins[j,1])**2-(kappa-pe*(1-kappa))**2)}}
se.kappa2 <- sqrt((sumterm1*sumterm2)/(sum(tabulated.data)*(1-pe)**2))
se.kappa <- sqrt((pe+pe**2-sumterm3)/(sum(tabulated.data)*(1-pe)**2))
tabulated.data <- addmargins(tabulated.data)
z.value <- kappa/se.kappa
p.value <- 1-pnorm(abs(z.value),0,1)
kappa.low <- kappa-qnorm((1-alpha/2),0,1)*se.kappa
kappa.upper <- kappa+qnorm((1-alpha/2),0,1)*se.kappa
res <-list(call=match.call(),contingencyTab=tabulated.data,kappa=kappa,alpha=alpha,kappaCI=c(kappa.low,kappa.upper),
zval=z.value,stdError=se.kappa,pval=2*p.value)
class(res) <- "dyadkappa"
res
}
print.dyadkappa <- function(x,digits=max(4,getOption("digits")-4),...)
{
cat("\n")
cat("Distinguishable Members with categorical data: Kappa Index")
cat("\n\n Call: \n")
cat("",deparse(x$call), "\n\n")
cat(" Contingency table: \n")
print.table(x$contingencyTab)
cat("\n\n")
cat("Agreement Test: ","\n")
results <- cbind(x$kappa,x$stdError,x$zval,x$pval)
colnames(results) <- c("Kappa Value", "Asymp. Std Error", "z Value", "Approx. Sig.")
rownames(results) <- ""
print.table(results,digits=digits,scientific=FALSE)
cat("\n\n")
cat(paste((1-x$alpha)*100,"% Confidence Interval: ",sep=''),"\n")
results <- t(x$kappaCI)
colnames(results) <- c("Lower","Upper")
rownames(results) <- ''
print.table(results,digits=digits)
cat("\n\n")
invisible(x)
}
# A function for estimating dyadic interdependence in standard dyadic designs with indistinguishable dyad members and categorical responses #
categorical.indistinguishable.dyad <- function(dataset,conf.int=0.95){
tabulated.data <- table(dataset[,2],dataset[,3])
tabulated.data <- (tabulated.data+t(tabulated.data))/2
alpha <- 1 - conf.int
no <- sum(diag(tabulated.data))
ne <- sum(margin.table(tabulated.data,1)*margin.table(tabulated.data,2)/sum(tabulated.data))
kappa <- (no - ne)/(sum(tabulated.data)-ne)
table.props <- prop.table(tabulated.data)
prop.margins <-array(c(margin.table(table.props,1),margin.table(prop.table(tabulated.data),2)),dim=c(nrow(tabulated.data),2))
sumterm1 <- 0.
sumterm2 <- 0.
sumterm3 <- 0.
pe <- ne/sum(tabulated.data)
for (i in 1:nrow(table.props)){
sumterm1 <- sumterm1+(table.props[i,i]*(1-sum(prop.margins[i,])*(1-kappa))**2+(1-kappa)**2)
sumterm3 <- sumterm3+(prop.margins[i,1]*prop.margins[i,2]*sum(prop.margins[i,]))}
for (i in 1:nrow(table.props)){
for (j in 1:ncol(table.props)){
sumterm2 <- sumterm2+(table.props[i,j]*(prop.margins[i,2]+prop.margins[j,1])**2-(kappa-pe*(1-kappa))**2)}}
se.kappa2 <- sqrt((sumterm1*sumterm2)/(sum(tabulated.data)*(1-pe)**2))
se.kappa <- sqrt((pe+pe**2-sumterm3)/(sum(tabulated.data)*(1-pe)**2))
tabulated.data <- addmargins(tabulated.data)
z.value <- kappa/se.kappa
p.value <- 1-pnorm(abs(z.value),0,1)
kappa.low <- kappa-qnorm((1-alpha/2),0,1)*se.kappa
kappa.upper <- kappa+qnorm((1-alpha/2),0,1)*se.kappa
res <-list(call=match.call(),contingencyTab=tabulated.data,kappa=kappa,alpha=alpha,kappaCI=c(kappa.low,kappa.upper),
zval=z.value,stdError=se.kappa,pval=2*p.value)
class(res) <- "dyadkappa2"
res
}
print.dyadkappa2 <- function(x,digits=max(4,getOption("digits")-4),...)
{
cat("\n")
cat("Indistinguishable Members with categorical data: Kappa Index")
cat("\n\n Call: \n")
cat("",deparse(x$call), "\n\n")
cat(" Contingency table: \n")
print.table(x$contingencyTab)
cat("\n\n")
cat("Agreement Test: ","\n")
results <- cbind(x$kappa,x$stdError,x$zval,x$pval)
colnames(results) <- c("Kappa Value", "Asymp. Std Error", "z Value", "Approx. Sig.")
rownames(results) <- ""
print.table(results,digits=digits)
cat("\n\n")
cat(paste((1-x$alpha)*100,"% Confidence Interval: ",sep=''),"\n")
results <- t(x$kappaCI)
colnames(results) <- c("Lower","Upper")
rownames(results) <- ''
print.table(results,digits=digits)
cat("\n\n")
invisible(x)
}
################################################################################
# R FUNCTIONS FOR SRM ROUND ROBIN DESIGNS #
################################################################################
# Function for formatting data according a round robin design
prepare.SRM.matrix <- function(X,numb.individuals,numb.times,numb.groups,names=NULL,times=NULL,groups=NULL){
if (is.null(names)) names <- paste('Ind.',c(1:numb.individuals))
if (is.null(times)) times <- paste('Time',c(1:numb.times))
if (is.null(groups)) groups <- paste('Group',c(1:numb.groups))
X <- array(X,c(numb.individuals,numb.individuals,numb.times,numb.groups))
temp <- array(0,c(numb.individuals,numb.individuals,numb.times,numb.groups))
for (i in 1:(dim(X)[1])){
for (j in 1:(dim(X)[1])){
for (k in 1:(dim(X)[3])){
for (l in 1:(dim(X)[4])){
temp[i,j,k,l] <- X[j,i,k,l]}}}}
count <- 1;
count2 <- 1;
for (k in 1:(dim(X)[3])){
for (i in 1:(dim(X)[1])){
if (count <= (dim(X)[3])){
X[count2,,count,] <- temp[i,,k,]
count <- count + 1}
else {
count <- 1
count2 <- count2 + 1
X[count2,,count,] <- temp[i,,k,]
count <- count + 1}}}
res<-list(data=X,individuals=numb.individuals,ntimes=numb.times,ngroups=numb.groups,
names=names,times=times,groups=groups)
class(res) <- 'srmRRMatrix'
res
}
# Function for obtaining SRM effects #
SRM.effects <- function (X){
if (class(X) != "srmRRMatrix")
stop("Data input need to be a srmRRMatrix object")
N <- X$individuals
T <- X$ntimes
G <- X$ngroups
# Compute row means #
row.means <- array(dim=c(N,G),0.)
for (l in 1:G){
row.means[,l] <- rowSums(X$data[,,,l])/(G*(N-1));}
time.rowmeans <- array(dim=c(N,T,G),0.)
for (k in 1:T){
for (l in 1:G){
time.rowmeans[,k,l] <- rowSums(X$data[,,k,l])/(N-1)}}
# Compute column means #
column.means <- array(dim=c(N,G),0.)
if (T <= 1) {
for (l in 1:G){
column.means[,l] <- colSums(X$data[,,,l])/(T*(N-1))}}
if (T > 1) {
for (l in 1:G){
column.means[,l] <- rowSums(colSums(X$data[,,,l]))/(T*(N-1))}}
time.colmeans <- array(dim=c(N,T,G),0.)
for (k in 1:T){
for (l in 1:G){
time.colmeans[,k,l] <- colSums(X$data[,,k,l])/(N-1)}}
# Compute grand mean #
grand.mean <- array(dim=c(G,1),0.)
for (l in 1:G){
grand.mean[l] <- sum(X$data[,,,l])/(T*(N*(N-1)))}
time.grandmean <- array(dim=c(T,G),0.)
for (k in 1:T){
for (l in 1:G){
time.grandmean[k,l] <- sum(X$data[,,k,l])/((N*(N-1)))}}
# Estimate actor effects for all individuals #
actor <- array(dim=c(N,T,G),dimnames=list(X$names,X$times,X$groups),0.)
for (k in 1:T){
for (l in 1:G){
for (i in 1:N){
actor[i,k,l] <- (N-1)*((N-1)*time.rowmeans[i,k,l]+time.colmeans[i,k,l]-
N*time.grandmean[k,l])/(N*(N-2))}}}
# Estimate partner effects for all individuals #
partner <- array(dim=c(N,T,G),dimnames=list(X$names,X$times,X$groups),0.)
for (k in 1:T){
for (l in 1:G){
for (i in 1:N){
partner[i,k,l] <- (N-1)*((N-1)*time.colmeans[i,k,l]+time.rowmeans[i,k,l]-
N*time.grandmean[k,l])/(N*(N-2))}}}
# Estimate relationships effects for all individuals #
relationship <- array(dim=c(N,N,T,G),
dimnames=list(X$names,X$names,X$times,X$groups),0.)
for (i in 1:N){
for (j in 1:N){
for (k in 1:T){
for (l in 1:G){
if (i != j)
{relationship[i,j,k,l] <- X$data[i,j,k,l]-actor[i,k,l]-partner[j,k,l]-time.grandmean[k,l]}}}}}
res<-list(actor.effects=actor,partner.effects=partner,relationship.effects=relationship)
class(res)="srmRREffects"
res
}
print.srmRREffects <- function(x,digits=max(4,getOption("digits")-4),...)
{
cat("\n")
cat("SRM Effects: Round Robin Design")
cat("\n\n")
N <- dim(x$actor.effects)[1]
T <- dim(x$actor.effects)[2]
G <- dim(x$actor.effects)[3]
names <- unlist(dimnames(x$actor.effects)[1])
times <- unlist(dimnames(x$actor.effects)[2])
groups <- unlist(dimnames(x$actor.effects)[3])
group <- rep(groups,2*N*T)
time <- rep(times,2*N*G)
individual <- rep(names,2)
type <- c(rep("actor",N*T*G),rep("partner",N*T*G))
effect <- c (x$actor.effects,x$partner.effects)
results<-data.frame(group,time,individual,type,effect)
cat("Actor and Partner Effects: \n")
print(results)
cat("\n\n")
results <- vector()
for (j in 1:G)
{
for (i in 1:T)
{
mat <- x$relationship.effects[,,i,j]
actornames <- c(rownames(mat)[(row(mat)*(row(mat)<col(mat)))[row(mat)<col(mat)]],
rownames(mat)[(row(mat)*(row(mat)>col(mat)))[row(mat)>col(mat)]])
partnernames <- c(colnames(mat)[(col(mat)*(row(mat)<col(mat)))[row(mat)<col(mat)]],
colnames(mat)[(col(mat)*(row(mat)>col(mat)))[row(mat)>col(mat)]])
effects <- mat[cbind(c((row(mat)*(row(mat)<col(mat)))[row(mat)<col(mat)],
(row(mat)*(row(mat)>col(mat)))[row(mat)>col(mat)]),
c((col(mat)*(row(mat)<col(mat)))[row(mat)<col(mat)],
(col(mat)*(row(mat)>col(mat)))[row(mat)>col(mat)]))]
group <- rep(groups[j],N*(N-1)/2)
time <- rep(times[i],N*(N-1)/2)
res <- data.frame(group,time,actor=actornames,partner=partnernames,relationship=effects)
results<- rbind(results,res)
}
}
results<-results[order(results$actor,results$partner),]
rownames(results) <- NULL
cat("Relationship Effects: \n")
print(results)
cat("\n\n")
invisible(x)
}
# Function to estimate SRM variances #
SRM.variances <- function(X){
if (class(X) != "srmRRMatrix")
stop("Data input need to be a srmRRMatrix object")
N <- X$individuals
T <- X$ntimes
G <- X$ngroups
names <- X$names
times <- X$times
groups <- X$groups
effects <- SRM.effects(X)
actor <- effects$actor.effects
partner <- effects$partner.effects
relationship <- effects$relationship.effects
# Estimate mean squares for actor, partner and cross products actor-partner #
MSactor <- array(dim=c(T,G),0.)
for (k in 1:T){
for (l in 1:G){
MSactor[k,l] <- sum(actor[,k,l]**2)/(N-1)}}
MSpartner <- array(dim=c(T,G),0.)
for (k in 1:T){
for (l in 1:G){
MSpartner[k,l] <- sum(partner[,k,l]**2)/(N-1)}}
MCP <- array(dim=c(T,G),0.)
for (k in 1:T){
for (l in 1:G){
MCP[k,l] <- sum(actor[,k,l]*partner[,k,l])/(N-1)}}
# Estimate mean squares between and within dyads #
sum.relationship <- array(dim=c(N,N,T,G),0.)
for (k in 1:T){
for (l in 1:G){
sum.relationship[,,k,l] <- relationship[,,k,l] + t(relationship[,,k,l])}}
MSbetween <- array(dim=c(T,G),dimnames=list(times,groups),0.)
for (k in 1:T){
for (l in 1:G){
MSbetween[k,l] <- (sum(sum.relationship[,,k,l]**2)/2)/((N-1)*(N-2)-2)}}
subs.relationship <- array(dim=c(N,N,T,G),0.)
for (k in 1:T){
for (l in 1:G){
subs.relationship[,,k,l] <- relationship[,,k,l] - t(relationship[,,k,l]);}}
MSwithin <- array(dim=c(T,G),dimnames=list(times,groups),0.)
for (k in 1:T){
for (l in 1:G){
MSwithin[k,l] <- (sum(subs.relationship[,,k,l]**2)/2)/((N-1)*(N-2))}}
# Estimate variance for actor effect #
variance.actor <- array(dim=c(T,G),dimnames=list(times,groups),0.)
for (k in 1:T){
for (l in 1:G){
variance.actor[k,l] <- MSactor[k,l] - 0.5*((MSbetween[k,l]/(N-2))+(MSwithin[k,l]/N))
if (variance.actor[k,l]<0) variance.actor[k,l]<- NA}}
# Estimate variance for partner effect #
variance.partner <- array(dim=c(T,G),dimnames=list(times,groups),0.)
for (k in 1:T){
for (l in 1:G){
variance.partner[k,l] <- MSpartner[k,l] - 0.5*((MSbetween[k,l]/(N-2))+(MSwithin[k,l]/N))
if (variance.partner[k,l]<0) variance.partner[k,l]<- NA}}
# Estimate covariance for actor-partner effect #
covariance.actorpartner <- array(dim=c(T,G),dimnames=list(times,groups),0.)
for (k in 1:T){
for (l in 1:G){
covariance.actorpartner[k,l] <- MCP[k,l] - 0.5*((MSbetween[k,l]/(N-2))-(MSwithin[k,l]/N))}}
# Estimate variance for relationship effect #
variance.relationship <- array(dim=c(T,G),dimnames=list(times,groups),0.)
for (k in 1:T){
for (l in 1:G){
variance.relationship[k,l] <- 0.5*(MSbetween[k,l]+MSwithin[k,l])
if (variance.relationship[k,l]<0) variance.relationship[k,l]<- NA}}
# Estimate covariance for relationship effect #
covariance.relationship <- array(dim=c(T,G),dimnames=list(times,groups),0.)
for (k in 1:T){
for (l in 1:G){
covariance.relationship[k,l] <- 0.5*(MSbetween[k,l]-MSwithin[k,l])}}
res <- list(MSbetween=MSbetween,MSwithin=MSwithin,actor.variance=variance.actor,partner.variance=variance.partner,
relationship.variance=variance.relationship,actorpartner.covariance=covariance.actorpartner,
dyadic.covariance=covariance.relationship)
class(res) <- "srmRRVars"
res
}
print.srmRRVars <- function(x,digits=max(4,getOption("digits")-4),...)
{
cat("\n")
cat("SRM Variances: Round Robin Design")
cat("\n\n")
T <- dim(x$actor.variance)[1]
G <- dim(x$actor.variance)[2]
times <- unlist(dimnames(x$actor.variance)[1])
groups <- unlist(dimnames(x$actor.variance)[2])
group <- rep(groups,T)
time <- rep(times,G)
actor <- array(x$actor.variance)
partner <- array(x$partner.variance)
relationship <- array(x$relationship.variance)
actor.partner <- array(x$actorpartner.covariance)
dyadic <- array(x$dyadic.covariance)
results<-data.frame(group,time,actor,partner,relationship,actor.partner,dyadic)
print(results)
cat("\n\n")
invisible(x)
}
# Function to compute relative variance due to actor, partner and relationship effects #
SRM.relative.variances <- function(X){
if (class(X) != "srmRRMatrix")
stop("Data input need to be a srmRRMatrix object")
N <- X$individuals
T <- X$ntimes
G <- X$ngroups
names <- X$names
times <- X$times
groups <- X$groups
variances <- SRM.variances(X)
variance.actor <- variances[[3]]
variance.partner <- variances[[4]]
variance.relationship <- variances[[5]]
prop.variance.actor <- array(dim=c(T,G),dimnames=list(times,groups),0.)
prop.variance.partner <- array(dim=c(T,G),dimnames=list(times,groups),0.)
prop.variance.relationship <- array(dim=c(T,G),dimnames=list(times,groups),0.)
total.variance <- rowSums(data.frame(variance.actor,variance.partner,variance.relationship),na.rm=TRUE)
prop.variance.actor <- variance.actor/total.variance
prop.variance.partner <- variance.partner/total.variance
prop.variance.relationship <- variance.relationship/total.variance
res<-list(actor.variance=prop.variance.actor,partner.variance=prop.variance.partner,
relationship.variance=prop.variance.relationship)
class(res) <- "srmRRrelVars"
res
}
print.srmRRrelVars <- function(x,digits=max(4,getOption("digits")-4),...)
{
cat("\n")
cat("SRM Relative Variances: Round Robin Design")
cat("\n\n")
T <- dim(x$actor.variance)[1]
G <- dim(x$actor.variance)[2]
times <- unlist(dimnames(x$actor.variance)[1])
groups <- unlist(dimnames(x$actor.variance)[2])
group <- rep(groups,T)
time <- rep(times,G)
actor <- array(x$actor.variance)
partner <- array(x$partner.variance)
relationship <- array(x$relationship.variance)
results<-data.frame(group,time,actor,partner,relationship)
rownames(results) <- ''
print(results)
cat("\n\n")
invisible(x)
}
# Function to compute reliability for the estimates of the actor and partner variances #
SRM.reliability <- function(X){
if (class(X) != "srmRRMatrix")
stop("Data input need to be a srmRRMatrix object")
N <- X$individuals
T <- X$ntimes
G <- X$ngroups
names <- X$names
times <- X$times
groups <- X$groups
variances <- SRM.variances(X)
variance.actor <- variances[[3]]
variance.partner <- variances[[4]]
variance.relationship <- variances[[5]]
covariance.relationship <- variances[[7]]
actor.reliability <- array(dim=c(T,G),dimnames=list(times,groups),0.)
actor.reliability <- variance.actor/(variance.actor + variance.relationship/(N-1)-
covariance.relationship/(N-1)**2)
partner.reliability <- array(dim=c(T,G),dimnames=list(times,groups),0.)
partner.reliability <- variance.partner/(variance.partner + variance.relationship/(N-1)-
covariance.relationship/(N-1)**2)
list(actor.reliability=actor.reliability,partner.reliability=partner.reliability)
}
# Testing SRM effects by means of Jackknife tests #
SRM.jackknife <- function(X){
if (class(X) != "srmRRMatrix")
stop("Data input need to be a srmRRMatrix object")
N <- X$individuals
T <- X$ntimes
G <- X$ngroups
names <- X$names
times <- X$times
groups <- X$groups
# Is group(s) size greater than or equal to 5 individuals? #
if (N < 5)
return("Error: Group(s) size must be greater than or equal to 5 individuals for carrying out the jackknife.")
source <- c("Actor Var","Partner Var","Relationship Var","Actor-Partner Cov","Dyadic Cov")
jackknife <- array(dim=c(N,5,T,G),dimnames=list(names,source,times,groups),0.)
jackknife.t.statistic <- array(dim=c(T,dim(jackknife)[2],G),dimnames=list(times,source,groups))
jackknife.p.value <- array(dim=c(T,dim(jackknife)[2],G),dimnames=list(times,source,groups))
jackknife.mean <- array(dim=c(T,dim(jackknife)[2],G),0.)
jackknife.variance <- array(dim=c(T,dim(jackknife)[2],G),0.)
Xtemp <- array(dim=c(N-1,N-1,T,G),0.)
variances <- SRM.variances(X)
variance.actor <- variances[[3]]
variance.partner <- variances[[4]]
variance.relationship <- variances[[5]]
actorpartner.covariance <- variances[[6]]
dyadic.covariance <- variances[[7]]
for (i in 1:N){
Xtemp <- as.numeric(X$data[-i,-i,,])
Xtemp <- prepare.SRM.matrix(Xtemp,(N-1),T,G)
variancestemp <- SRM.variances(Xtemp)
for (k in 1:T){
for (l in 1:G){
jackknife[i,1,k,l] <- N*variance.actor[k,l]-variancestemp$actor.variance[k,l]*(N-1)
jackknife[i,2,k,l] <- N*variance.partner[k,l]-variancestemp$partner.variance[k,l]*(N-1)
jackknife[i,3,k,l] <- N*variance.relationship[k,l]-variancestemp$relationship.variance[k,l]*(N-1)
jackknife[i,4,k,l] <- N*actorpartner.covariance[k,l]-variancestemp$actorpartner.covariance[k,l]*(N-1)
jackknife[i,5,k,l] <- N*dyadic.covariance[k,l]-variancestemp$dyadic.covariance[k,l]*(N-1)}}}
for (k in 1:T){
for (l in 1:G){
jackknife.mean[k,,l] <- apply(jackknife[,,k,l],2,mean)
jackknife.variance[k,,l] <- apply(jackknife[,,k,l],2,var)/N
jackknife.t.statistic[k,,l] <- jackknife.mean[k,,l]/sqrt(jackknife.variance[k,,l])}}
df <- N-1
for (k in 1:T){
for (m in 1:(dim(jackknife)[2])){
for (l in 1:G){
if ( !is.na(jackknife.t.statistic[k,m,l])){
if (sign(jackknife.t.statistic[k,m,l])>0)
jackknife.p.value[k,m,l] <- pt(jackknife.t.statistic[k,m,l],df,lower.tail=FALSE)
else jackknife.p.value[k,m,l] <- pt(jackknife.t.statistic[k,m,l],df)}}}}
list(t.statistic=jackknife.t.statistic,df=df,two.tailed.p.value=2*jackknife.p.value)
}
# Function to estimate generalized reciprocity #
SRM.generalized.reciprocity <- function(X){
if (class(X) != "srmRRMatrix")
stop("Data input need to be a srmRRMatrix object")
N <- X$individuals
T <- X$ntimes
G <- X$ngroups
names <- X$names
times <- X$times
groups <- X$groups
effects <- SRM.effects(X)
actor <- effects$actor.effects
partner <- effects$partner.effects
variances <- SRM.variances(X)
variance.actor <- variances[[3]]
variance.partner <- variances[[4]]
MSbetween <- variances[[1]]
MSwithin <- variances[[2]]
generalized.reciprocity <- array(dim=c(T,G),dimnames=list(times,groups),0.)
for (k in 1:T){
for (l in 1:G){
if (is.na(variance.actor[k,l]) | is.na(variance.partner[k,l])) generalized.reciprocity[k,l] <- NA
else{
generalized.reciprocity[k,l] <- (sum(actor[,k,l]*partner[,k,l])/(N-1) - MSbetween[k,l]/(2*(N-2)) +
MSwithin[k,l]/(2*N))/sqrt(variance.actor[k,l]*variance.partner[k,l])
if ((generalized.reciprocity[k,l] > 1.0) | (generalized.reciprocity[k,l] < -1.0))
generalized.reciprocity[k,l] <- sign(generalized.reciprocity[k,l])*1.0}}}
return(generalized.reciprocity)
}
# Compute dyadic.reciprocity #
SRM.dyadic.reciprocity <- function(X){
if (class(X) != "srmRRMatrix")
stop("Data input need to be a srmRRMatrix object")
N <- X$individuals
T <- X$ntimes
G <- X$ngroups
names <- X$names
times <- X$times
groups <- X$groups
variances <- SRM.variances(X)
MSbetween <- variances[[1]]
MSwithin <- variances[[2]]
dyadic.reciprocity <- array(dim=c(T,G),dimnames=list(times,groups),0.)
for (k in 1:T){
for (l in 1:G){
dyadic.reciprocity[k,l] <- (MSbetween[k,l] - MSwithin[k,l])/(MSbetween[k,l] + MSwithin[k,l])}}
return(dyadic.reciprocity)
}
# Testing dyadic covariance by means of a F ratio #
SRM.dyadic.covariance.Ftest <- function(X){
if (class(X) != "srmRRMatrix")
stop("Data input need to be a srmRRMatrix object")
N <- X$individuals
T <- X$ntimes
G <- X$ngroups
names <- X$names
times <- X$times
groups <- X$groups
variances <- SRM.variances(X)
MSbetween <- variances[[1]]
MSwithin <- variances[[2]]
F.value <- array(dim=c(T,G),dimnames=list(times,groups))
for (k in 1:T){
for (l in 1:G){
F.value[k,l] <- MSbetween[k,l]/MSwithin[k,l]}}
df1 <- (N-1)*(N-2)/2-1
df2 <- (N-1)*(N-2)/2
F.p.value <- array(dim=c(T,G),dimnames=list(times,groups))
for (k in 1:T){
for (l in 1:G){
if (sign(F.value[k,l])>0)
F.p.value[k,l] <- pf(F.value[k,l],df1,df2,lower.tail=FALSE)
else F.p.value[k,l] <- pf(F.value[k,l],df1,df2)}}
res=list(MSb=MSbetween,MSw=MSwithin,F.value=F.value,df1=df1,df2=df2,p.value=F.p.value)
class(res) <- "srmCovFTest"
res
}
print.srmCovFTest <- function(x,digits=max(4,getOption("digits")-4),...)
{
cat("\n")
cat("SRM Dyadic Covariance F Test: Round Robin Design")
cat("\n\n")
T <- dim(x$MSb)[1]
G <- dim(x$MSb)[2]
times <- unlist(dimnames(x$MSb)[1])
groups <- unlist(dimnames(x$MSb)[2])
group <- rep(groups,T)
time <- rep(times,G)
MSb <- x[[1]]
MSw <- x[[2]]
df1 <- x[[4]]
df2 <- x[[5]]
F <- x[[3]]
pval <- x[[6]]
results<-data.frame(group,time,MSb,MSw,df1,df2,F, pval)
colnames(results)<-c("Group","Time","MSb","MSw","df1","df2","F value",if (sign(F)>0) "Pr(> F)" else "Pr(< F)")
rownames(results) <- ''
print(results)
cat("\n\n")
invisible(x)
}
# Testing SRM effects by means of a between-group t Test #
SRM.between.groups.tTest <- function(X){
if (class(X) != "srmRRMatrix")
stop("Data input need to be a srmRRMatrix object")
N <- X$individuals
T <- X$ntimes
G <- X$ngroups
names <- X$names
times <- X$times
# Are there more than 1 group? #
if (G < 2)
return("Error: Number of groups analyzed should be greater than 1.")
variances <- SRM.variances(X)
variance.actor <- variances[[3]]
variance.partner <- variances[[4]]
variance.relationship <- variances[[5]]
actorpartner.covariance <- variances[[6]]
dyadic.covariance <- variances[[7]]
source <- c("Actor Var","Partner Var","Relationship Var","Actor-Partner Cov","Dyadic Cov")
between.group.t.statistic <- array(dim=c(T,5),dimnames=list(times,source))
between.group.p.value <- array(dim=c(T,5),dimnames=list(times,source))
actor.variance.mean <- apply(variance.actor,1,mean)
partner.variance.mean <- apply(variance.partner,1,mean)
relationship.variance.mean <- apply(variance.relationship,1,mean)
actorpartner.covariance.mean <- apply(actorpartner.covariance,1,mean)
dyadic.covariance.mean <- apply(dyadic.covariance,1,mean)
actor.variance.var <- apply(variance.actor,1,var)
partner.variance.var <- apply(variance.partner,1,var)
relationship.variance.var <- apply(variance.relationship,1,var)
actorpartner.covariance.var <- apply(actorpartner.covariance,1,var)
dyadic.covariance.var <- apply(dyadic.covariance,1,var)
for (k in 1:T){
between.group.t.statistic[k,1] <- actor.variance.mean[k]/sqrt(actor.variance.var[k]/G)
between.group.t.statistic[k,2] <- partner.variance.mean[k]/sqrt(partner.variance.var[k]/G)
between.group.t.statistic[k,3] <- relationship.variance.mean[k]/sqrt(relationship.variance.var[k]/G)
between.group.t.statistic[k,4] <- dyadic.covariance.mean[k]/sqrt(dyadic.covariance.var[k]/G)
between.group.t.statistic[k,5] <- actorpartner.covariance.mean[k]/sqrt(actorpartner.covariance.var[k]/G)}
df <- G-1
for (k in 1:T){
for (m in 1:(dim(between.group.t.statistic)[2])){
if ( !is.na(between.group.t.statistic[k,m])){
if (sign(between.group.t.statistic[k,m]) > 0)
between.group.p.value[k,m] <- pt(between.group.t.statistic[k,m],df,lower.tail=FALSE)
else between.group.p.value[k,m] <- pt(between.group.t.statistic[k,m],df)}}}
list(t.statistic=between.group.t.statistic,df=df,two.tailed.p.value=2*between.group.p.value)
}
# Testing SRM effects by means of a within-group t Test #
SRM.within.groups.tTest <- function(X){
if (class(X) != "srmRRMatrix")
stop("Data input need to be a srmRRMatrix object")
N <- X$individuals
T <- X$ntimes
G <- X$ngroups
names <- X$names
times <- X$times
groups <- X$groups
# Computing some terms that will be used below #
h1 <- 2*(G**2*N**6-7*G**2*N**5+18*G**2*N**4+10*G*N**4-20*G**2*N**3-
46*G*N**3+8*G**2*N**2+70*G*N**2+24*N**2-36*G*N-48*N+32)
h2 <- 2*(G**2*N**5-5*G**2*N**4+2*G*N**4+8*G**2*N**3-6*G*N**3-
4*G**2*N**2+14*G*N**2+8*N**2-20*G*N-16*N+32)
h3 <- 4*(G**2*N**5-5*G**2*N**4+8*G**2*N**3+12*G*N**3-
4*G**2*N**2-38*G*N**2+28*G*N+32*N-32)
h4 <- (N**2)*((N-2)**2)*(G*N-G+2)*(G*N**2-3*G*N+4)*(G*N**2-3*G*N+
2*G+4)
h5 <- G**3*N**7-7*G**3*N**6+19*G**3*N**5+10*G**2*N**5-25*G**3*N**4-
54*G**2*N**4-4*G*N**4+16*G**3*N**3+120*G**2*N**3+44*G*N**3-
4*G**3*N**2-132*G**2*N**2-124*G*N**2-16*N**2+56*G**2*N+
168*G*N+32*N-64*G-64
h6 <- G**3*N**7-7*G**3*N**6-2*G**2*N**6+19*G**3*N**5+26*G**2*N**5-
25*G**3*N**4-100*G**2*N**4-20*G*N**4+16*G**3*N**3+176*G**2*N**3+
124*G*N**3-4*G**3*N**2-156*G**2*N**2-236*G*N**2-48*N**2+56*G**2*N+
200*G*N+96*N-64*G-64
# Estimating mean and variance of round robin parameters #
variances <- SRM.variances(X)
variance.actor <- unlist(variances[[3]])
variance.partner <- unlist(variances[[4]])
variance.relationship <- unlist(variances[[5]])
actorpartner.covariance <- unlist(variances[[6]])
apcor <- unlist(SRM.generalized.reciprocity(X))
dyadcor <- unlist(SRM.dyadic.reciprocity(X))
dyadic.covariance <- unlist(variances[[7]])
actor.variance.mean <- apply(variance.actor,1,function(x) mean(x,na.rm=TRUE))
partner.variance.mean <- apply(variance.partner,1,function(x) mean(x,na.rm=TRUE))
relationship.variance.mean <- apply(variance.relationship,1,function(x) mean(x,na.rm=TRUE))
actorpartner.covariance.mean <- apply(actorpartner.covariance,1,function(x) mean(x,na.rm=TRUE))
dyadic.covariance.mean <- apply(dyadic.covariance,1,function(x) mean(x,na.rm=TRUE))
source <- c("Actor Var","Partner Var","Relationship Var","Actor-Partner Cov","Dyadic Cov")
SRM.variance.parameters <- array(dim=c(T,5),dimnames=list(times,source),0.)
within.group.t.statistic <- array(dim=c(T,5),dimnames=list(times,source),0.)
within.group.p.value <- array(dim=c(T,5),dimnames=list(times,source),0.)
for (k in 1:T){
SRM.variance.parameters[k,1] <- ((2*(actor.variance.mean[k])**2)/(G*N-G+2))+(((4*actor.variance.mean[k])*
((N-1)*relationship.variance.mean[k]+dyadic.covariance.mean[k]))/(N*(N-2)*(G*N-G+2)))+
((h1*(relationship.variance.mean[k])**2+h2*(dyadic.covariance.mean[k])**2+h3*
(relationship.variance.mean[k]*dyadic.covariance.mean[k]))/h4)
SRM.variance.parameters[k,2] <- ((2*(partner.variance.mean[k])**2)/(G*N-G+2))+(((4*partner.variance.mean[k])*
((N-1)*relationship.variance.mean[k]+dyadic.covariance.mean[k]))/(N*(N-2)*(G*N-G+2)))+
((h1*(relationship.variance.mean[k])**2+h2*(dyadic.covariance.mean[k])**2+h3*
(relationship.variance.mean[k]*dyadic.covariance.mean[k]))/h4)
SRM.variance.parameters[k,3] <- SRM.variance.parameters[k,5] <- ((2*((G*N**2)-(3*G*N)+G+4)*
(relationship.variance.mean[k]**2+dyadic.covariance.mean[k]**2))/(((G*N**2)-(3*G*N)+4)*
((G*N**2)-(3*G*N)+2*G+4)))+((4*G*relationship.variance.mean[k]*dyadic.covariance.mean[k])/(((G*N**2)-(3*G*N)+4)*
((G*N**2)-(3*G*N)+2*G+4)))
SRM.variance.parameters[k,4] <- ((G*N-G-2)*actorpartner.covariance.mean[k]**2+G*(N-1)*
actor.variance.mean[k]*partner.variance.mean[k])/((G*N-G+2)*(G*N-G-1))+(h5*relationship.variance.mean[k]**2+
h6*dyadic.covariance.mean[k]**2+h3*relationship.variance.mean[k]*dyadic.covariance.mean[k])/
((G*N-G-1)*h4)+(G*(N-1)*(actor.variance.mean[k]+partner.variance.mean[k])*((N-1)*
relationship.variance.mean[k]+dyadic.covariance.mean[k]))/(N*(N-2)*(G*N-G+2)*(G*N-G-1))+
(2*(G*N-G-2)*actorpartner.covariance.mean[k]*(relationship.variance.mean[k]+(N-1)*dyadic.covariance.mean[k]))/
(N*(N-2)*(G*N-G+2)*(G*N-G-1))}
SRM.relative.parameters <- array(dim=c(T,5),dimnames=list(times,source),0.)
SRM.mean.parameters <- array(dim=c(T,5),dimnames=list(times,source),c(actor.variance.mean,
partner.variance.mean,relationship.variance.mean,actorpartner.covariance.mean,dyadic.covariance.mean))
SRM.relative.parameters[1:3] <- SRM.mean.parameters[1:3]/sum(unlist(SRM.mean.parameters[1:3]),na.rm=TRUE)
SRM.relative.parameters[4] <- apply(apcor,1,function(x) mean(x,na.rm=TRUE))
SRM.relative.parameters[5] <- apply(dyadcor,1,function(x) mean(x,na.rm=TRUE))
SRM.stderror.parameters <- sqrt(SRM.variance.parameters)
for (k in 1:T){
within.group.t.statistic[k,1] <- actor.variance.mean[k]/sqrt(SRM.variance.parameters[k,1])
within.group.t.statistic[k,2] <- partner.variance.mean[k]/sqrt(SRM.variance.parameters[k,2])
within.group.t.statistic[k,3] <- relationship.variance.mean[k]/sqrt(SRM.variance.parameters[k,3])
within.group.t.statistic[k,4] <- actorpartner.covariance.mean[k]/sqrt(SRM.variance.parameters[k,4])
within.group.t.statistic[k,5] <- dyadic.covariance.mean[k]/sqrt(SRM.variance.parameters[k,5])}
df <- G*(N-1)
for (k in 1:T){
for (m in 1:(dim(within.group.t.statistic)[2])){
if ( !is.na(within.group.t.statistic[k,m])){
if (sign(within.group.t.statistic[k,m]) > 0)
within.group.p.value[k,m] <- pt(within.group.t.statistic[k,m],df,lower.tail=FALSE)
else within.group.p.value[k,m] <- pt(within.group.t.statistic[k,m],df)}}}
res <- list(estimates=SRM.mean.parameters,standard.values=SRM.relative.parameters,standard.error=SRM.stderror.parameters,
t.value=within.group.t.statistic,df=df,p.value=within.group.p.value)
class(res) <-"srmWithintTest"
res
}
print.srmWithintTest <- function(x,digits=max(4,getOption("digits")-4),...)
{
cat("\n")
cat("Within group t-Test: Round Robin Design")
cat("\n\n")
estimates <- as.numeric(x$estimates)
std.value <- as.numeric(x$standard.values)
std.error <- as.numeric(x$standard.error)
t <- as.numeric(x$t.value)
df <- as.numeric(x$df)
p.value<- as.numeric(x$p.value)
results<-data.frame(estimates,std.value,std.error,t,df,p.value)
results[is.na(results)]<-NA
results[is.na(results[,2]),3:6]<-NA
rownames(results) <- c('actor var','partner var','relationship var','actor-partner cov','dyadic cov')
print(results)
cat("\n\n")
invisible(x)
}
################################################################################
# R FUNCTIONS FOR SRM BLOCK DESIGNS #
################################################################################
# Function for preparing data in order to be analized as a SRM block design #
prepare.block.matrix <- function (X,subgroup1,subgroup2,numb.times,numb.groups,
names1=NULL,names2=NULL,times=NULL,groups=NULL,
subgroups=NULL){
numb.individuals <- subgroup1 + subgroup2
if (is.null(names1)) names1 <- paste('Ind.',c(1:(subgroup1)))
else names1 <- names1
if (is.null(names2)) names2 <- paste('Ind.',c((numb.individuals-subgroup2+1):numb.individuals))
else names2 <- names2
if (is.null(times)) times <- paste('Time',c(1:numb.times))
else times <- times
if (is.null(groups)) groups <- paste('Group',c(1:numb.groups))
else groups <- groups
if (is.null(subgroups)) subgroups <- paste('Subgroup',c(1,2))
else subgroups <- subgroups
X <- array(X,c(numb.individuals,numb.individuals,numb.times,numb.groups),
dimnames=list(c(names1,names2),c(names1,names2),times,groups))
temp <- array(0,c(numb.individuals,numb.individuals,numb.times,numb.groups))
for (i in 1:(dim(X)[1])){
for (j in 1:(dim(X)[2])){
for (k in 1:(dim(X)[3])){
for (l in 1:(dim(X)[4])){
temp[i,j,k,l] <- X[j,i,k,l]}}}}
count <- 1;
count2 <- 1;
for (k in 1:(dim(X)[3])){
for (i in 1:(dim(X)[1])){
if (count <= (dim(X)[3])){
X[count2,,count,] <- temp[i,,k,]
count <- count + 1}
else {
count <- 1
count2 <- count2 + 1
X[count2,,count,] <- temp[i,,k,]
count <- count + 1}}}
res<-list(data=X,individuals1=subgroup1,individuals2=subgroup2,
individuals=numb.individuals,ntimes=numb.times,ngroups=numb.groups,
names1=names1,names2=names2,times=times,groups=groups,
subgroups=subgroups)
class(res) <- 'srmBlockMatrix'
res
}
# Function for obtaining SRM effects in Block designs #
block.SRM.effects <- function (X){
if (class(X) != "srmBlockMatrix")
stop("Data input need to be a srmBlockMatrix object")
N <- X$individuals
T <- X$ntimes
G <- X$ngroups
G1 <- X$individuals1
G2 <- X$individuals2
names1 <- X$names1
names2 <- X$names2
groups<- X$groups
times <- X$times
# Compute row means #
time.rowmeans <- array(dim=c(N,T,G),0.)
for (k in 1:T){
for (l in 1:G){
for (i in 1:N){
if (i <= G1) time.rowmeans[i,k,l] <- sum(X$data[i,,k,l])/G2
else time.rowmeans[i,k,l] <- sum(X$data[i,,k,l])/G1}}}
# Compute column means #
time.colmeans <- array(dim=c(N,T,G),0.)
for (k in 1:T){
for (l in 1:G){
for (j in 1:N){
if (j <= G1) time.colmeans[j,k,l] <- sum(X$data[,j,k,l])/G2
else time.colmeans[j,k,l] <- sum(X$data[,j,k,l])/G1}}}
# Compute grand mean #
time.grandmean <- array(dim=c(2,T,G),0.)
for (i in 1:N){
for (j in 1:N){
for (k in 1:T){
for (l in 1:G){
if ((i <= G1) && (j > G1)){
time.grandmean[1,k,l] <- time.grandmean[1,k,l] + X$data[i,j,k,l]}
if ( (i > G1) && (j <= G1)){
time.grandmean[2,k,l] <- time.grandmean[2,k,l] + X$data[i,j,k,l]}}}}}
time.grandmean <- time.grandmean/(G1*G2)
# Estimate actor effects for all individuals #
actor <- array(dim=c(N,T,G),dimnames=list(c(names1,names2),times,groups),0.)
for (k in 1:T){
for (l in 1:G){
for (i in 1:N){
if (i <= G1){
actor[i,k,l] <- time.rowmeans[i,k,l]-time.grandmean[1,k,l]}
else actor[i,k,l] <- time.rowmeans[i,k,l]-time.grandmean[2,k,l]}}}
# Estimate partner effects for all individuals #
partner <- array(dim=c(N,T,G),dimnames=list(c(names1,names2),times,groups),0.)
for (j in 1:N){
for (k in 1:T){
for (l in 1:G){
if (j > G1){
partner[j,k,l] <- time.colmeans[j,k,l]-time.grandmean[1,k,l]}
if (j <= G1){
partner[j,k,l] <- time.colmeans[j,k,l]-time.grandmean[2,k,l]}}}}
# Estimate relationships effects for all individuals #
relationship <- array(dim=c(N,N,T,G),
dimnames=list(c(names1,names2),c(names1,names2),times,groups),0.)
for (i in 1:N){
for (j in 1:N){
for (k in 1:T){
for (l in 1:G){
if ((i <= G1) && (j > G1))
{relationship[i,j,k,l] <- X$data[i,j,k,l]-actor[i,k,l]-partner[j,k,l]-time.grandmean[1,k,l]}
if ( (i > G1) && (j <= G1))
{relationship[i,j,k,l] <- X$data[i,j,k,l]-actor[i,k,l]-partner[j,k,l]-time.grandmean[2,k,l]}}}}}
list(actor.effects=actor,partner.effects=partner,relationship.effects=relationship)
}
# Function to estimate SRM variances in block designs #
block.SRM.variances <- function(X){
if (class(X) != "srmBlockMatrix")
stop("Data input need to be a srmBlockMatrix object")
N <- X$individuals
T <- X$ntimes
G <- X$ngroups
G1 <- X$individuals1
G2 <- X$individuals2
names1 <- X$names1
names2 <- X$names2
groups <- X$groups
subgroups <- X$subgroups
times <- X$times
effects <- block.SRM.effects(X)
actor <- effects$actor.effects
partner <- effects$partner.effects
relationship <- effects$relationship.effects
# Estimate variance for relationship effect #
variance.relationship <- array(dim=c(2,T,G),dimnames=list(subgroups,times,groups),0.)
for (i in 1:N){
for (j in 1:N){
for (k in 1:T){
for (l in 1:G){
if ((i <= G1) && (j > G1))
{variance.relationship[1,k,l] <- variance.relationship[1,k,l]+relationship[i,j,k,l]**2}
if ( (i > G1) && (j <= G1))
{variance.relationship[2,k,l] <- variance.relationship[2,k,l]+relationship[i,j,k,l]**2}}}}}
variance.relationship <- variance.relationship/((G1-1)*(G2-1))
# Estimate variance for actor effect #
variance.actor <- array(dim=c(2,T,G),dimnames=list(subgroups,times,groups),0.)
for (i in 1:N){
for (k in 1:T){
for (l in 1:G){
if ( i <= G1 )
{variance.actor[1,k,l] <- variance.actor[1,k,l]+actor[i,k,l]**2}
if ( i > G1 )
{variance.actor[2,k,l] <- variance.actor[2,k,l]+actor[i,k,l]**2}}}}
variance.actor <- (variance.actor/(G1-1))-(variance.relationship/G2)
for (k in 1:T){
for (l in 1:G){
for (m in 1:length(subgroups)){
if ( sign(variance.actor[m,k,l]) < 0. ) variance.actor[m,k,l] <- 0.}}}
# Estimate variance for partner effect #
variance.partner <- array(dim=c(2,T,G),dimnames=list(subgroups,times,groups),0.)
for (j in 1:N){
for (k in 1:T){
for (l in 1:G){
if ( j <= G2 )
{variance.partner[1,k,l] <- variance.partner[1,k,l]+partner[j,k,l]**2}
if ( j > G1 )
{variance.partner[2,k,l] <- variance.partner[2,k,l]+partner[j,k,l]**2}}}}
variance.partner[1,,] <- (variance.partner[1,,]/(G2-1))-(variance.relationship[2]/G2)
variance.partner[2,,] <- (variance.partner[2,,]/(G2-1))-(variance.relationship[1]/G2)
for (k in 1:T){
for (l in 1:G){
for (m in 1:length(subgroups)){
if ( sign(variance.partner[m,k,l]) < 0. ) variance.partner[m,k,l] <- 0.}}}
# Estimate covariance for actor-partner effect #
dyads1 <- array(dim=c(G1*G2,T,G),0.)
dyads2 <- array(dim=c(G1*G2,T,G),0.)
for (k in 1:T){
for (l in 1:G){
count1 <- 1
count2 <- 1
for (i in 1:N){
for (j in 1:N){
if ((i <= G1) && (j > G1))
{dyads1[count1,k,l] <- relationship[i,j,k,l]
count1 <- count1 +1}
if ( (i > G1) && (j <= G1))
{dyads2[count2,k,l] <- relationship[i,j,k,l]
count2 <- count2 +1}}}}}
covariance.actorpartner <- array(dim=c(T,G),dimnames=list(times,groups),0.)
for (k in 1:T){
for (l in 1:G){
covariance.actorpartner[k,l] <- cov(dyads1[,k,l],dyads2[,k,l])}}
# Estimate covariance for relationship effect #
covariance.relationship <- array(dim=c(2,T,G),dimnames=list(subgroups,times,groups),0.)
crosspod.relationship <- array(dim=c(N,N,T,G),0.)
for (k in 1:T){
for (l in 1:G){
crosspod.relationship[,,k,l] <- relationship[,,k,l] * t(relationship[,,k,l])}}
for (i in 1:N){
for (k in 1:T){
for (l in 1:G){
if ( i <= G1 )
{covariance.relationship[1,k,l] <- covariance.relationship[1,k,l]+actor[i,k,l]*partner[i,k,l]}
if ( i > (N-G2) )
{covariance.relationship[2,k,l] <- covariance.relationship[2,k,l]+actor[i,k,l]*partner[i,k,l]}}}}
for (k in 1:T){
for (l in 1:G){
covariance.relationship[1,k,l] <- (covariance.relationship[1,k,l]/(G2-1))-(sum(crosspod.relationship[,,k,l])/
2/((G2-1)*(G1-1)))/G2
covariance.relationship[2,k,l] <- (covariance.relationship[2,k,l]/(G1-1))-(sum(crosspod.relationship[,,k,l])/
2/((G2-1)*(G1-1)))/G2}}
list(actor.variance=variance.actor,partner.variance=variance.partner,relationship.variance=variance.relationship,
actorpartner.covariance=covariance.actorpartner,dyadic.covariance=covariance.relationship)
}
# Function to compute relative variance due to actor, partner and relationship effects i block designs #
block.SRM.relative.variances <- function(X){
if (class(X) != "srmBlockMatrix")
stop("Data input need to be a srmBlockMatrix object")
N <- X$individuals
T <- X$ntimes
G <- X$ngroups
G1 <- X$individuals1
G2 <- X$individuals2
names1 <- X$names1
names2 <- X$names2
groups <- X$groups
subgroups <- X$subgroups
times <- X$times
variances <- block.SRM.variances(X)
variance.actor <- variances$actor.variance
variance.partner <- variances$partner.variance
variance.relationship <- variances$relationship.variance
prop.variance.actor <- array(dim=c(2,T,G),dimnames=list(subgroups,times,groups),0.)
prop.variance.partner <- array(dim=c(2,T,G),dimnames=list(subgroups,times,groups),0.)
prop.variance.relationship <- array(dim=c(2,T,G),dimnames=list(subgroups,times,groups),0.)
total.variance <- variance.actor + variance.partner + variance.relationship
prop.variance.actor <- variance.actor/total.variance
prop.variance.partner <- variance.partner/total.variance
prop.variance.relationship <- variance.relationship/total.variance
list(relative.actor.variance=prop.variance.actor,relative.partner.variance=prop.variance.partner,
relative.relationship.variance=prop.variance.relationship)
}
# Function to compute reliability for the estimates of the actor and partner variances in block designs #
block.SRM.reliability <- function(X){
if (class(X) != "srmBlockMatrix")
stop("Data input need to be a srmBlockMatrix object")
N <- X$individuals
T <- X$ntimes
G <- X$ngroups
G1 <- X$individuals1
G2 <- X$individuals2
names1 <- X$names1
names2 <- X$names2
subgroups <- X$subgroups
groups <- X$sgroups
times <- X$times
variances <- block.SRM.variances(X)
variance.actor <- variances$actor.variance
variance.partner <- variances$partner.variance
variance.relationship <- variances$relationship.variance
actor.reliability <- array(dim=c(2,T,G),dimnames=list(subgroups,times,groups),0)
partner.reliability <- array(dim=c(2,T,G),dimnames=list(subgroups,times,groups),0.)
for (k in 1:T){
for (l in 1:G){
if ( (variance.actor[1,k,l] != 0) | (variance.relationship[1,k,l] != 0) )
actor.reliability[1,k,l] <- variance.actor[1,k,l]/(variance.actor[1,k,l] +
variance.relationship[1,k,l]/G1)
if ( (variance.actor[2,k,l] != 0) | (variance.relationship[2,k,l] != 0) )
actor.reliability[2,k,l] <- variance.actor[2,k,l]/(variance.actor[2,k,l] +
variance.relationship[2,k,l]/G2)
if ( (variance.partner[1,k,l] != 0) | (variance.relationship[1,k,l] != 0) )
partner.reliability[1,k,l] <- variance.partner[1,k,l]/(variance.partner[1,k,l] +
variance.relationship[1,k,l]/G1)
if ( (variance.partner[2,k,l] != 0) | (variance.relationship[2,k,l] != 0) )
partner.reliability[2,k,l] <- variance.partner[2,k,l]/(variance.partner[2,k,l] +
variance.relationship[1,k,l]/G2)}}
list(actor.reliability=actor.reliability,partner.reliability=partner.reliability)
}
# Function to estimate generalized reciprocity in block designs #
block.SRM.generalized.reciprocity <- function(X){
if (class(X) != "srmBlockMatrix")
stop("Data input need to be a srmBlockMatrix object")
N <- X$individuals
T <- X$ntimes
G <- X$ngroups
G1 <- X$individuals1
G2 <- X$individuals2
names1 <- X$names1
names2 <- X$names2
subgroups<- X$subgroups
groups <- X$groups
times <- X$times
variances <- block.SRM.variances(X)
variance.actor <- variances$actor.variance
variance.partner <- variances$partner.variance
covariance.relationship <- variances$dyadic.covariance
generalized.reciprocity <- array(dim=c(2,T,G),dimnames=list(subgroups,times,groups),0.)
for (m in 1:length(subgroups)){
for (k in 1:T){
for (l in 1:G){
if ((variance.actor[m,k,l] == 0.) | (variance.partner[m,k,l] == 0.)) generalized.reciprocity[m,k,l] <- 0.
else
generalized.reciprocity[m,k,l] <- covariance.relationship[m,k,l]/sqrt(variance.actor[m,k,l]*variance.partner[m,k,l])}}
if ((generalized.reciprocity[m,k,l] > 1.0) | (generalized.reciprocity[m,k,l] < -1.0))
generalized.reciprocity[m,k,l] <- sign(generalized.reciprocity[m,k,l])*1.0}
return(generalized.reciprocity)
}
# Compute dyadic.reciprocity in block designs #
block.SRM.dyadic.reciprocity <- function(X,symmetric=FALSE){
if (class(X) != "srmBlockMatrix")
stop("Data input need to be a srmBlockMatrix object")
N <- X$individuals
T <- X$ntimes
G <- X$ngroups
G1 <- X$individuals1
G2 <- X$individuals2
names1 <- X$names1
names2 <- X$names2
subgroups<- X$subgroups
groups <- X$groups
times <- X$times
effects <- block.SRM.effects(X)
relationship <- effects$relationship.effects
variances <- block.SRM.variances(X)
variance.relationship <- variances$relationship.variance
# Estimate covariance for actor-partner effect #
dyads1 <- array(dim=c(G1*G2,T,G),0.)
dyads2 <- array(dim=c(G1*G2,T,G),0.)
for (k in 1:T){
for (l in 1:G){
count1 <- 1
count2 <- 1
for (i in 1:N){
for (j in 1:N){
if ((i <= G1) && (j > G1))
{dyads1[count1,k,l] <- relationship[i,j,k,l]
count1 <- count1 +1}
if ( (i > G1) && (j <= G1))
{dyads2[count2,k,l] <- relationship[i,j,k,l]
count2 <- count2 +1}}}}}
dyadic.reciprocity <- array(dim=c(T,G),dimnames=list(times,groups),0.)
if (symmetric == FALSE) {
for (k in 1:T){
for (l in 1:G){
dyadic.reciprocity[k,l] <- cor(dyads1[,k,l],dyads2[,k,l])}}}
else {
for (k in 1:T){
for (l in 1:G){
dyadic.reciprocity[k,l] <- cor(dyads1[,k,l],dyads2[,k,l])*(sqrt(variance.relationship[1,k,l]*
variance.relationship[2,k,l])/((variance.relationship[1,k,l]+variance.relationship[2,k,l])/2))}}}
return(dyadic.reciprocity)
}
# Testing SRM effects in block designs by means of Jackknife tests #
block.SRM.jackknife <- function(X){
if (class(X) != "srmBlockMatrix")
stop("Data input need to be a srmBlockMatrix object")
N <- X$individuals
T <- X$ntimes
G <- X$ngroups
G1 <- X$individuals1
G2 <- X$individuals2
names1 <- X$names1
names2 <- X$names2
subgroups<- X$subgroups
groups <- X$groups
times <- X$times
# Is group(s) size greater than or equal to 6 individuals? #
if (N < 6)
return("Error: Group(s) size must be greater than or equal to 6 individuals for carrying out the jackknife for a block design.")
if (G1 < 3 | G2 < 3)
return("Error: Subgroups must be greater than or equal to 3 individuals for carrying out the jackknife for a block design.")
source <- c("Actor Var","Partner Var","Relationship Var","Actor-Partner Cov","Dyadic Cov")
jackknife <- array(dim=c(length(subgroups),N,5,T,G),dimnames=list(subgroups,c(names1,names2),
source,times,groups),0.)
jackknife.t.statistic <- array(dim=c(length(subgroups),T,dim(jackknife)[3],
G),dimnames=list(subgroups,times,source,groups),0.)
jackknife.p.value <- array(dim=c(length(subgroups),T,dim(jackknife)[3],G),
dimnames=list(subgroups,times,source,groups),NA)
jackknife.mean <- array(dim=c(length(subgroups),T,dim(jackknife)[3],G),0.)
jackknife.variance <- array(dim=c(length(subgroups),T,dim(jackknife)[3],G),0.)
Xtemp <- array(dim=c(N-1,N-1,T,G),0.)
variances <- block.SRM.variances(X)
variance.actor <- variances$actor.variance
variance.partner <- variances$partner.variance
variance.relationship <- variances$relationship.variance
actorpartner.covariance <- variances$actorpartner.covariance
dyadic.covariance <- variances$dyadic.covariance
for (m in 1:length(subgroups)){
count <- 1
for (i in 1:N){
if (count <= G1) {
ind.temp.1 <- G1-1
ind.temp.2 <- G2}
else {
ind.temp.1 <- G1
ind.temp.2 <- G2-1}
if (count <= G1){
term1 <- G1
term2 <- ind.temp.1}
else {term1 <- G2
term2 <- ind.temp.2}
count <- count + 1
Xtemp <- prepare.block.matrix(as.numeric(X$data[-i,-i,,]),N-1,N-1,T,G)
variancestemp <- block.SRM.variances(Xtemp)
for (k in 1:T){
for (l in 1:G){
jackknife[m,i,1,k,l] <- term1*variance.actor[m,k,l]-variancestemp$actor.variance[m,k,l]*term2
jackknife[m,i,2,k,l] <- term1*variance.partner[m,k,l]-variancestemp$partner.variance[m,k,l]*term2
jackknife[m,i,3,k,l] <- term1*variance.relationship[m,k,l]-variancestemp$relationship.variance[m,k,l]*term2
jackknife[m,i,4,k,l] <- term1*dyadic.covariance[m,k,l]-variancestemp$dyadic.covariance[m,k,l]*term2
jackknife[m,i,5,k,l] <- term1*actorpartner.covariance[k,l]-variancestemp$actorpartner.covariance[k,l]*term2}}}}
for (m in 1:length(subgroups)){
for (k in 1:T){
for (l in 1:G){
jackknife.mean[m,k,,l] <- apply(jackknife[m,,,k,l],2,mean)
jackknife.variance[m,k,,l] <- apply(jackknife[m,,,k,l],2,var)/(N)
jackknife.t.statistic[m,k,,l] <- jackknife.mean[m,k,,l]/sqrt(jackknife.variance[m,k,,l])}}}
df1 <- G1-1
df2 <- G2-2
df <- array(dim=c(2,1),dimnames=list(subgroups),c(df1,df2))
for (m in 1:length(subgroups)){
for (k in 1:T){
for (n in 1:(dim(jackknife)[3])){
for (l in 1:G){
if (sign(jackknife.t.statistic[m,k,n,l])>0){
if (m == 1)
jackknife.p.value[m,k,n,l] <- pt(jackknife.t.statistic[m,k,n,l],df1,lower.tail=FALSE)
if (m == 2)
jackknife.p.value[m,k,n,l] <- pt(jackknife.t.statistic[m,k,n,l],df2,lower.tail=FALSE)}
else{
if (m == 1)
jackknife.p.value[m,k,n,l] <- pt(jackknife.t.statistic[m,k,n,l],df1)
if (m == 2)
jackknife.p.value[m,k,n,l] <- pt(jackknife.t.statistic[m,k,n,l],df2)}}}}}
list(t.statistic=jackknife.t.statistic,df=df,two.tailed.p.value=2*jackknife.p.value)
}
# Testing SRM effects in block designs by means of a between-group t Test #
block.SRM.between.groups.tTest <-function(X){
if (class(X) != "srmBlockMatrix")
stop("Data input need to be a srmBlockMatrix object")
N <- X$individuals
T <- X$ntimes
G <- X$ngroups
G1 <- X$individuals1
G2 <- X$individuals2
names1 <- X$names1
names2 <- X$names2
subgroups<- X$subgroups
groups <- X$groups
times <- X$times
# Are there more than 1 group? #
if (G < 2)
return("Error: Number of groups analyzed should be greater than 1.")
variances <- block.SRM.variances(X)
variance.actor <- variances$actor.variance
variance.partner <- variances$partner.variance
variance.relationship <- variances$relationship.variance
actorpartner.covariance <- variances$actorpartner.covariance
dyadic.covariance <- variances$dyadic.covariance
source <- c("Actor Var","Partner Var","Relationship Var","Actor-Partner Cov","Dyadic Cov")
between.group.t.statistic <- array(dim=c(2,T,5),dimnames=list(subgroups,times,source))
between.group.p.value <- array(dim=c(2,T,5),dimnames=list(subgroups,times,source))
actor.variance.mean <- array(dim=c(2,T))
partner.variance.mean <- array(dim=c(2,T))
relationship.variance.mean <- array(dim=c(2,T))
actorpartner.covariance.mean <- array(dim=c(T))
dyadic.covariance.mean <- array(dim=c(2,T))
actor.variance.var <- array(dim=c(2,T))
partner.variance.var <- array(dim=c(2,T))
relationship.variance.var <- array(dim=c(2,T))
actorpartner.covariance.var <- array(dim=c(T))
dyadic.covariance.var <- array(dim=c(2,T))
for (k in 1:T){
actor.variance.mean[,k] <- apply(variance.actor[,k,],1,mean)
partner.variance.mean[,k] <- apply(variance.partner[,k,],1,mean)
relationship.variance.mean[,k] <- apply(variance.relationship[,k,],1,mean)
actorpartner.covariance.mean[k] <- mean(actorpartner.covariance[k,])
dyadic.covariance.mean[,k] <- apply(dyadic.covariance[,k,],1,mean)
actor.variance.var[,k] <- apply(variance.actor[,k,],1,var)
partner.variance.var[,k] <- apply(variance.partner[,k,],1,var)
relationship.variance.var[,k] <- apply(variance.relationship[,k,],1,var)
actorpartner.covariance.var[k] <- var(actorpartner.covariance[k,])
dyadic.covariance.var[,k] <- apply(dyadic.covariance[,k,],1,var)}
for (m in 1:length(subgroups)){
for (k in 1:T){
between.group.t.statistic[m,k,1] <- actor.variance.mean[m,k]/sqrt(actor.variance.var[m,k]/G)
between.group.t.statistic[m,k,2] <- partner.variance.mean[m,k]/sqrt(partner.variance.var[m,k]/G)
between.group.t.statistic[m,k,3] <- relationship.variance.mean[m,k]/sqrt(relationship.variance.var[m,k]/G)
between.group.t.statistic[m,k,4] <- dyadic.covariance.mean[m,k]/sqrt(dyadic.covariance.var[m,k]/G)
between.group.t.statistic[m,k,5] <- actorpartner.covariance.mean[k]/sqrt(actorpartner.covariance.var[k]/G)}}
df <- G-1
for (m in 1:(dim(between.group.t.statistic)[1])){
for (k in 1:T){
for (n in 1:(dim(between.group.t.statistic)[3])){
if ( !is.na(between.group.t.statistic[m,k,n]) ){
if ( sign(between.group.t.statistic[m,k,n]) > 0 )
between.group.p.value[m,k,n] <- pt(between.group.t.statistic[m,k,n],df,lower.tail=FALSE)
else between.group.p.value[m,k,n] <- pt(between.group.t.statistic[m,k,n],df)}}}}
list(t.statistic=between.group.t.statistic,df=df,two.tailed.p.value=2*between.group.p.value)
}
# Testing dyadic covariance in block designs by means of r Pearson correlation test #
block.SRM.dyadic.reciprocity.tTest <- function(X,symmetric=FALSE){
if (class(X) != "srmBlockMatrix")
stop("Data input need to be a srmBlockMatrix object")
N <- X$individuals
T <- X$ntimes
G <- X$ngroups
G1 <- X$individuals1
G2 <- X$individuals2
names1 <- X$names1
names2 <- X$names2
subgroups<- X$subgroups
groups <- X$groups
times <- X$times
dyadic.reciprocity <- block.SRM.dyadic.reciprocity(X,symmetric)
t.value <- array(dim=c(T,G),dimnames=list(times,groups),0.)
for (k in 1:T){
for (l in 1:G){
t.value[k,l] <- dyadic.reciprocity[k,l]*sqrt(G1*(G2-1)-1)/
sqrt(1-dyadic.reciprocity[k,l]**2)}}
df <- (G1-1)*(G2-1)-1
t.p.value <- array(dim=c(T,G),dimnames=list(times,groups),NA)
for (k in 1:T){
for (l in 1:G){
if (sign(t.value[k,l])>0)
t.p.value[k,l] <- pt(t.value[k,l],df,lower.tail=FALSE)
else t.p.value[k,l] <- pt(t.value[k,l],df)}}
list(t.value=t.value,two.tailed.p.value=2*t.p.value)
}
################################################################################
# R FUNCTIONS FOR MEASURING RECIPROCITY WITH DYADIC DISCREPANCIES #
################################################################################
# Function to obtain DC index by means of the observed sociomatrix #
getdc <- function(X){
N = sum(X);
dif_n = abs(X-t(X));
sumdif <- sum(dif_n)
dc = sumdif/(2*N)
return(dc)
}
# Function to obtain delta index #
getdelta <- function(X){
psi <- getpsi(X);
phi <- getphi(X);
delta = phi/psi;
return(delta)
}
# Function to compute generalized reciprocity index -epsilon- #
getepsilon <- function(X){
nu <- getnu(X);
nuj <- getnuj(X);
nui <- colSums(nu);
epsilon = 0.;
for (i in 1:nrow(X)){
epsilon = epsilon + abs(nuj[i] - nui[i]);}
if (sum(nuj) != 0)
{epsilon = 1 - (epsilon/(sum(nuj)))}
else
{epsilon=1.}
return(epsilon)
}
# Function to compute dyadic reciprocity index -kappa- #
getkappa <- function(X){
rationu <- getrationu(X);
dif_rationu = abs(rationu - t(rationu));
kappa =1 - ((sum(dif_rationu)/2)/nrow(X));
return(kappa)
}
# Function to compute matrix of unweighted contributions to symmetry -lambda- #
getlambda <- function(X){
S = (X + t(X))/2;
trX = sum(diag(t(X)%*%X));
Srow = S*S;
vectss = rowSums(Srow);
K = (X - t(X))/2;
Krow = K*K;
vectkk = rowSums(Krow);
lambda = array(dim=c(nrow(X),ncol(X)));
sumvect = vectss+vectkk;
for(i in 1:nrow(X))
for(j in 1:ncol(X))
{if (sumvect[i] == 0)
{lambda[i,j] = 0.}
else
{lambda[i,j] = (S[i,j]*S[i,j])/sumvect[i]}
}
lambda;
}
# Function to compute individual contributions to symmetry -lambdaj- #
getlambdaj <- function(X){
lambda <- getlambda(X);
lambdaj = rowSums(lambda);
lambdaj;
}
# Function to compute matrix of unweighted contributions to skew-symmetry -nu- #
getnu <- function(X){
S = (X + t(X))/2;
trX = sum(diag(t(X)%*%X));
Srow = S*S;
vectss = rowSums(Srow);
K = (X - t(X))/2;
Krow = K*K;
vectkk = rowSums(Krow);
nu = array(dim=c(nrow(X),ncol(X)));
sumvect = vectss+vectkk;
for(i in 1:nrow(X))
for(j in 1:ncol(X))
{if (sumvect[i] == 0)
{nu[i,j] = 0.}
else
{nu[i,j] = (K[i,j]*K[i,j])/sumvect[i]}
}
nu;
}
# Function to compute individual contributions to skew-symmetry -nuj- #
getnuj <- function(X){
nu <- getnu(X);
nuj = rowSums(nu);
nuj;
}
# Function to compute matrix of ratios skew-symmetry/symmetry -omega- #
getomega <- function(X){
nu <- getnu(X);
lambda <- getlambda(X);
omega = array(dim=c(nrow(X),ncol(X)));
for(i in 1:nrow(X))
for(j in 1:ncol(X))
{if (lambda[i,j] == 0.)
{omega[i,j]=0.}
else
{omega[i,j] = nu[i,j]/lambda[i,j]}
}
omega;
}
# Function to obtain phi index by means of the observed sociomatrix #
getphi <- function(X){
K = (X - t(X))/2;
trX = sum(diag(t(X)%*%X));
trK = sum(diag(t(K)%*%K));
phi = trK / trX
return(phi)
}
# Function to obtain psi index by means of the observed sociomatrix #
getpsi <- function(X){
S = (X + t(X))/2;
trX = sum(diag(t(X)%*%X));
trS = sum(diag(t(S)%*%S));
psi = trS / trX
return(psi)
}
# Function to compute dyadic decomposition of symmetry -ratiolambda- #
getratiolambda <- function(X){
lambda <- getlambda(X);
lambdaj <- getlambdaj(X);
ratiolambda = array(dim=c(nrow(X),ncol(X)));
for(i in 1:nrow(X))
for(j in 1:ncol(X))
{if (lambdaj[i] == 0)
{ratiolambda[i,j] = 0.}
else
{ratiolambda[j,i] = lambda[i,j]/lambdaj[i]}
}
ratiolambda;
}
# Function to compute dyadic decomposition of skew-symmetry -rationu- #
getrationu <- function(X){
nu <- getnu(X);
nuj <- getnuj(X);
rationu = array(dim=c(nrow(X),ncol(X)));
for(i in 1:nrow(X))
for(j in 1:ncol(X))
{if (nuj[i] == 0)
{rationu[i,j] = 0.}
else
{rationu[j,i] = nu[i,j]/nuj[i]}
}
rationu;
}
# Function to estimate p values for social reciprocity statistics at different levels of analysis under the specified null hypothesis #
reciptest1<-function(X,pi,rep=9999,names=NULL,label=FALSE){
if (!squaremat(X))
return("Error: Matrix X is not square and cannot be analyzed")
if (!squaremat(pi))
return("Error: Matrix Pi is not square and cannot be analyzed")
# Limit the number of replications #
if ((rep < 1) | (rep > 1000000))
return("Error: Number of replications must be between 1 and 1000000");
# Compute the social reciprocity statistics at different levels for the original matrix #
phi <- getphi(X);
psi <- getpsi(X);
dc <- getdc(X);
delta <- getdelta(X);
epsilon <- getepsilon(X);
kappa <- getkappa(X);
nuj <- getnuj(X);
lambdaj <- getlambdaj(X);
rationu <- getrationu(X);
ratiolambda <- getratiolambda(X);
omega <- getomega(X);
# Matrices have to be transformed into vectors in order to pass to a .C routine #
vecX=array(dim=c(nrow(X)*(ncol(X)-1)/2,1));
m=0.;
for (i in 1:nrow(X))
for (j in 1:ncol(X))
{
m =m+1;
vecX[m]=X[i,j];}
vecpi=array(dim=c(nrow(X)*(ncol(X)-1)/2,1));
m=0.;
for (i in 1:nrow(pi))
for (j in 1:ncol(pi))
{
m =m+1;
vecpi[m]=pi[i,j];}
# R wrapper of the .C routine #
out <- .C("recip0",
as.double(vecX),
as.double(vecpi),
as.integer(nrow(X)),
as.integer(rep),
pvphi=double(1),
pvpsi=double(1),
pvdc=double(1),
pvdelta=double(1),
pvepsilon=double(1),
pvkappa=double(1),
pvnuj=double(nrow(X)),
pvlambdaj=double(nrow(X)),
pvrationu=double(length(vecX)),
pvratiolambda=double(length(vecX)),
pvomega=double(length(vecX)),
PACKAGE="DyaDA")
pvphi = out$pvphi;
pvpsi = out$pvpsi;
pvdelta = out$pvdelta;
pvdc = out$pvdc;
pvepsilon = out$pvepsilon;
pvkappa = out$pvkappa;
pvrationu = array(dim=c(nrow(X),ncol(X)),0.);
m = 1;
for (i in 1:nrow(X))
for (j in 1:ncol(X))
{
pvrationu[i,j] = out$pvrationu[m];
m = m+1;
}
pvratiolambda = array(dim=c(nrow(X),ncol(X)),0.);
m = 1;
for (i in 1:nrow(X))
for (j in 1:ncol(X))
{
pvratiolambda[i,j] = out$pvratiolambda[m];
m = m+1;
}
pvomega = array(dim=c(nrow(X),ncol(X)),0.);
m = 1;
for (i in 1:nrow(X))
for (j in 1:ncol(X))
{
pvomega[i,j] = out$pvomega[m];
m = m+1;
}
pvnuj = out$pvnuj;
pvlambdaj = out$pvlambdaj;
# Computation of some summary statistics and print results #
if (label == TRUE){dimnames(rationu)<- list(c(names),c(names))}
if (label == TRUE){dimnames(pvrationu)<- list(c(names),c(names))}
if (label == TRUE){dimnames(ratiolambda)<- list(c(names),c(names))}
if (label == TRUE){dimnames(pvratiolambda)<- list(c(names),c(names))}
if (label == TRUE){dimnames(omega)<- list(c(names),c(names))}
if (label == TRUE){dimnames(pvomega)<- list(c(names),c(names))}
res <- list(label=label,names=names,phi=phi,pvphi=pvphi,psi=psi,pvpsi=pvpsi,delta=delta,pvdelta=pvdelta,
dc=dc,pvdc=pvdc,epsilon=epsilon,pvepsilon=pvepsilon,kappa=kappa,pvkappa=pvkappa,nuj=nuj,
pvnuj=pvnuj,lambdaj=lambdaj,pvlambdaj=pvlambdaj,ratiolambda=ratiolambda,
pvratiolambda=pvratiolambda,rationu=rationu,pvrationu=pvrationu,omega=omega,
pvomega=pvomega)
class(res) <- "reciptest1"
res
}
print.reciptest1 <- function(x,digits=max(4,getOption("digits")-4),...)
{
cat("SOCIAL RECIPROCITY ANALYSIS AT DIFFERENT LEVELS OF ANALYSIS","\n")
cat("===========================================================","\n")
cat(" ","\n")
cat("OVERALL MEASURES OF SOCIAL RECIPROCITY","\n")
cat("======================================","\n")
cat("Overall skew-symmetry index -phi-","\n")
cat("=================================","\n")
cat(round(x$phi,digits),"\n");
cat("p value","\n")
cat("=======","\n")
cat(round(x$pvphi,digits),"\n")
cat(" ","\n")
cat("Overall symmetry index -psi-","\n")
cat("============================","\n")
cat(round(x$psi,digits),"\n");
cat("p value","\n")
cat("=======","\n")
cat(round(x$pvpsi,digits),"\n")
cat(" ","\n")
cat("Overall skew-symmetry/symmetry ratio -delta-","\n")
cat("============================================","\n")
cat(round(x$delta,digits),"\n");
cat("p value","\n")
cat("=======","\n")
cat(round(x$pvdelta,digits),"\n")
cat(" ","\n")
cat("Directional consistency index -dc-","\n")
cat("==================================","\n")
cat(round(x$dc,digits),"\n");
cat("p value","\n")
cat("=======","\n")
cat(round(x$pvdc,digits),"\n")
cat(" ","\n")
cat("Generalized reciprocity index -epsilon-","\n")
cat("=======================================","\n")
cat(round(x$epsilon,digits),"\n");
cat("p value","\n")
cat("=======","\n")
cat(round(x$pvepsilon,digits),"\n")
cat(" ","\n")
cat("Dyadic reciprocity index -kappa-","\n")
cat("================================","\n")
cat(round(x$kappa,digits),"\n");
cat("p value","\n")
cat("=======","\n")
cat(round(x$pvkappa,digits),"\n")
cat(" ","\n")
cat(" ","\n")
cat("INDIVIDUAL MEASURES OF SOCIAL RECIPROCITY","\n")
cat("=========================================","\n")
cat(" ","\n")
nujmat = array(dim=c(nrow(X),2));
if (x$label == TRUE){dimnames(nujmat)<- list(c(x$ames),c("nuj","p value"))}
if (x$label == FALSE){colnames(nujmat)=c("nuj","p value")}
nujmat[,1] = round(x$nuj,digits);
nujmat[,2] = round(x$pvnuj,digits);
cat("Individual contribution to skew-symmetry","\n")
cat("========================================","\n")
print(nujmat);
cat(" ","\n")
lambdajmat = array(dim=c(nrow(X),2));
if (x$label == TRUE){dimnames(lambdajmat)<- list(c(x$names),c("lambdaj","p value"))};
if (x$label == FALSE){colnames(lambdajmat)=c("lambdaj","p value")};
lambdajmat[,1] = round(x$lambdaj,digits);
lambdajmat[,2] = round(x$pvlambdaj,digits);
cat("Individual contribution to symmetry","\n")
cat("===================================","\n")
print(lambdajmat);
cat(" ","\n")
cat(" ","\n")
cat("DYADIC MEASURES OF SOCIAL RECIPROCITY","\n")
cat("=====================================","\n")
cat(" ","\n")
cat("Dyadic contributions to skew-symmetry","\n")
cat("=====================================","\n")
print.table(x$rationu,digits);
cat("p value","\n")
cat("=======","\n")
print.table(x$pvrationu,digits);
cat(" ","\n")
cat("Dyadic contributions to symmetry","\n")
cat("================================","\n")
print.table(x$ratiolambda,digits);
cat("p value","\n")
cat("=======","\n")
print.table(x$pvratiolambda,digits);
cat(" ","\n")
cat("Matrix of dyadic balanced reciprocity","\n")
cat("=====================================","\n")
print.table(x$omega,digits);
cat("p value","\n")
cat("=======","\n")
cat(" ","\n")
print.table(x$pvomega,digits);
cat(" ","\n")
cat(" ","\n")
invisible(x)
}
# Function to compute the dyadic contribution to asymmetry - PHIij - #
getPHIij <- function (X)
{
phirmat <- getphirmat(X);
PHIij = phirmat + t(phirmat);
return(PHIij)
}
# Function to compute the overall asymmetry index - PHIr - #
getPHIr <- function (X)
{
phirmat <- getphirmat(X);
PHIr = sum(phirmat);
return(PHIr)
}
# Function to compute standard error of the PHIr statistic - PHIrSE - #
getPHIrSE <- function (X,pi)
{
dyadc = X + t(X);
z = 0.;
for (i in 1:nrow(X))
for (j in 1:ncol(X))
if (i < j)
{
if (dyadc[i,j] != 0)
z = z + 1;
}
oddterm = 0.;
for (i in 1:nrow(X))
for (j in 1:ncol(X))
if (i < j)
{
if ((dyadc[i,j]%%2) != 0)
{oddterm = oddterm + (1/(dyadc[i,j]**2))}
else {oddterm = oddterm}
}
oddterm = oddterm/2;
min = (z/2) + oddterm;
qqplus2s = 0.;
for (i in 1:nrow(X))
for (j in 1:ncol(X))
if (i < j)
{
qqplus2s = qqplus2s + ((4*dyadc[i,j]*((pi[i,j] - 7*(pi[i,j]**2) + 12*(pi[i,j]**3)) - 6*(pi[i,j]**4)) +
+ 8*(dyadc[i,j]**2)*(-pi[i,j] + 6*(pi[i,j]**2) - 10*(pi[i,j]**3) + 5*(pi[i,j]**4)) +
+ 4*(dyadc[i,j]**3)*(pi[i,j] - 5*(pi[i,j]**2) + 8*(pi[i,j]**3) - 4*(pi[i,j]**4)))/(dyadc[i,j]**4));
}
PHIrSE = 2*sqrt(qqplus2s)/(2*z-2*min);
return(PHIrSE)
}
# Function to compute mathematical expectancy of the PHIr statistic - PHIrexpec - #
getPHIrexpec <- function (X,pi)
{
dyadc = X + t(X);
z = 0.;
for (i in 1:nrow(X))
for (j in 1:ncol(X))
if (i < j)
{
if (dyadc[i,j] != 0)
z = z + 1;
}
oddterm = 0.;
for (i in 1:nrow(X))
for (j in 1:ncol(X))
if (i < j)
{
if ((dyadc[i,j]%%2) != 0)
{oddterm = oddterm + (1/(dyadc[i,j]**2))}
else {oddterm = oddterm}
}
oddterm = oddterm/2;
min = (z/2) + oddterm;
firstterm = 0.;
for (i in 1:nrow(X))
for (j in 1:ncol(X))
if (i < j)
{
firstterm = firstterm + (1/dyadc[i,j])
}
secondterm = 0.;
for (i in 1:nrow(X))
for (j in 1:ncol(X))
{
if (dyadc[i,j] == 0.)
{secondterm = secondterm + 0.}
if (dyadc[i,j] != 0.)
{secondterm = secondterm + (((pi[i,j]**2)*(dyadc[i,j]-1))/dyadc[i,j])}
}
PHIrexpec = (2*(firstterm + secondterm - min))/(2*z-2*min);
return(PHIrexpec)
}
# Function to compute the individual contribution to asymmetry as actor - phii - #
getphii <- function (X)
{
phirmat <- getphirmat(X);
phii = rowSums(phirmat);
return(phii)
}
# Function to compute the individual contribution to asymmetry as partner - phij - #
getphij <- function (X)
{
phirmat <- getphirmat(X);
phij = colSums(phirmat);
return(phij)
}
# Function to compute the matrix of directional dyadic asymmetry measures - phirmat - #
getphirmat <- function (X)
{
dyadc = X + t(X);
z = 0.;
for (i in 1:nrow(X))
for (j in i:ncol(X))
if (i < j)
{
if (dyadc[i,j] != 0)
z = z + 1;
}
max = z;
oddterm = 0.;
for (i in 1:nrow(X))
for (j in 1:ncol(X))
if (i < j)
{
if ((dyadc[i,j]%%2) != 0)
{oddterm = oddterm + (1/(dyadc[i,j]**2))}
else {oddterm = oddterm}
}
oddterm = oddterm/2;
min = (z/2) + oddterm;
phirmat = array(dim=c(nrow(X),ncol(X)),0.);
for (i in 1:nrow(X))
for (j in 1:ncol(X))
if (i != j)
{
if ( ((dyadc[i,j]%%2) > 0) & (dyadc[i,j] != 0) )
{phirmat[i,j] =(((4*(X[i,j]**2)/(dyadc[i,j]**2)-1))-(1/dyadc[i,j]**2))/(4*max-4*min)}
if ( ((dyadc[i,j]%%2) == 0) & (dyadc[i,j] != 0) )
{phirmat[i,j] =((4*(X[i,j]**2)/(dyadc[i,j]**2)-1))/(4*max-4*min)}
}
return(phirmat)
}
# Function to estimate p values for asymmetry statistics at different levels of analysis under the null hypothesis #
reciptest2<-function(X,pi,rep=9999,names=NULL,label=FALSE){
if (!squaremat(X))
return("Error: Matrix X is not square and cannot be analyzed");
if ( is.na(X) || !is.numeric(X))
return("Error: Sociomatrix must be numeric");
if (!squaremat(pi))
return("Error: Matrix Pi is not square and cannot be analyzed");
if ( is.na(pi) || !is.numeric(pi))
return("Error: Matrix of probabilities must be numeric");
# Limit the number of replications #
if ((rep < 1) | (rep > 1000000))
return("Error: Number of replications must be between 1 and 1000000");
# Compute the asymmetry statistics at different levels for the original matrix #
PHIr <- getPHIr(X);
expec <- getPHIrexpec(X,pi);
SE <- getPHIrSE(X,pi);
phirmat <- getphirmat(X);
PHIij <- getPHIij(X);
phii <- getphii(X);
phij <- getphij(X);
# Matrices have to be transformed into vectors in order to pass to a .C routine #
vecX=array(dim=c(nrow(X)*(ncol(X)-1)/2,1));
m=0.;
for (i in 1:nrow(X))
for (j in 1:ncol(X))
{
m =m+1;
vecX[m]=X[i,j];}
vecpi=array(dim=c(nrow(X)*(ncol(X)-1)/2,1));
m=0.;
for (i in 1:nrow(pi))
for (j in 1:ncol(pi))
{
m =m+1;
vecpi[m]=pi[i,j];}
# R wrapper of the .C routine #
out <- .C("recip",
as.double(vecX),
as.double(vecpi),
as.integer(nrow(X)),
as.integer(rep),
pvPHIrright=double(1),
pvPHIrleft=double(1),
pvphirmatright=double(length(vecX)),
pvphirmatleft=double(length(vecX)),
pvphiiright=double(nrow(X)),
pvphiileft=double(nrow(X)),
pvphijright=double(nrow(X)),
pvphijleft=double(nrow(X)),
pvPHIijright=double(length(vecX)),
pvPHIijleft=double(length(vecX)),
PACKAGE="DyaDA")
pvPHIrright = out$pvPHIrright;
pvPHIrleft = out$pvPHIrleft;
pvphirmatright = array(dim=c(nrow(X),ncol(X)),0.);
m = 1;
for (i in 1:nrow(X))
for (j in 1:ncol(X))
{
pvphirmatright[i,j] = out$pvphirmatright[m];
m = m+1;
}
pvphirmatleft = array(dim=c(nrow(X),ncol(X)),0.);
m = 1;
for (i in 1:nrow(X))
for (j in 1:ncol(X))
{
pvphirmatleft[i,j] = out$pvphirmatleft[m];
m = m+1;
}
pvPHIijright = array(dim=c(nrow(X),ncol(X)),0.);
m = 1;
for (i in 1:nrow(X))
for (j in 1:ncol(X))
{
pvPHIijright[i,j] = out$pvPHIijright[m];
m = m+1;
}
pvPHIijleft = array(dim=c(nrow(X),ncol(X)),0.);
m = 1;
for (i in 1:nrow(X))
for (j in 1:ncol(X))
{
pvPHIijleft[i,j] = out$pvPHIijleft[m];
m = m+1;
}
pvphiiright = out$pvphiiright;
pvphiileft = out$pvphiileft;
pvphijright = out$pvphijright;
pvphijleft = out$pvphijleft;
res <- list(label=label,names=names,PHIr=PHIr,pvPHIrright=pvPHIrright,
pvPHIrleft=pvPHIrleft,expec=expec,SE=SE,
phii=phii,pvphiiright=pvphiiright,pvphiileft=pvphiileft,
phij=phij,pvphijright=pvphijright,pvphijleft=pvphijleft,
phirmat=phirmat,pvphirmatright=pvphirmatright,pvphirmatleft=pvphirmatleft,
PHIij=PHIij,pvPHIijright=pvPHIijright,pvPHIijleft=pvPHIijleft)
class(res) <- "reciptest2"
res
}
print.reciptest2 <- function(x,digits=max(4,getOption("digits")-4),...)
{
# Computation of some summary statistics and print results #
if (x$label == TRUE){
dimnames(x$phirmat)<- list(c(x$names),c(x$names))
dimnames(x$pvphirmatright)<- list(c(x$names),c(x$names))
dimnames(x$pvphirmatleft)<- list(c(x$names),c(x$names))
dimnames(x$PHIij)<- list(c(x$names),c(x$names))
dimnames(x$pvPHIijright)<- list(c(x$names),c(x$names))
dimnames(x$pvPHIijleft)<- list(c(x$names),c(x$names))}
cat("SOCIAL RECIPROCITY ANALYSIS AT DIFFERENT LEVELS OF ANALYSIS","\n")
cat("===========================================================","\n")
cat(" ","\n")
cat("OVERALL MEASURES OF SOCIAL RECIPROCITY","\n")
cat("======================================","\n")
cat("PHIr","\n")
cat("====","\n")
cat(round(x$PHIr,digits),"\n");
cat("Right p value","\n")
cat("=============","\n")
cat(round(x$pvPHIrright,digits),"\n")
cat("Left p value","\n")
cat("============","\n")
cat(round(x$pvPHIrleft,digits),"\n")
cat(" ","\n")
cat("Mathematical expectancy of PHIr","\n")
cat("===============================","\n")
cat(round(x$expec,digits),"\n");
cat("Standard error of PHIr","\n")
cat("======================","\n")
cat(round(x$SE,digits),"\n")
cat(" ","\n")
cat(" ","\n")
cat("INDIVIDUAL MEASURES OF SOCIAL RECIPROCITY","\n")
cat("=========================================","\n")
phiimat = array(dim=c(nrow(X),3));
if (x$label == TRUE){dimnames(phiimat)<- list(c(x$names),c("phii","Right p value","Left p value"))}
if (x$label == FALSE){colnames(phiimat)=c("phii","Right p value","Left p value")}
phiimat[,1] = round(x$phii,digits);
phiimat[,2] = round(x$pvphiiright,digits);
phiimat[,3] = round(x$pvphiileft,digits);
cat("Individual contribution to asymmetry as actor","\n")
cat("=============================================","\n")
print.table(phiimat,digits);
phijmat = array(dim=c(nrow(X),3));
if (x$label == TRUE){dimnames(phijmat)<- list(c(x$names),c("phij","Right p value","Left p value"))};
if (x$label == FALSE){colnames(phijmat)=c("phij","Right p value","Left p value")};
phijmat[,1] = round(x$phij,digits);
phijmat[,2] = round(x$pvphijright,digits);
phijmat[,3] = round(x$pvphijleft,digits);
cat(" ","\n")
cat("Individual contribution to asymmetry as partner","\n")
cat("===============================================","\n")
print.table(phijmat);
cat(" ","\n")
cat(" ","\n")
cat("DYADIC MEASURES OF SOCIAL RECIPROCITY","\n")
cat("=====================================","\n")
cat("Dyadic directional asymmetry","\n")
cat("============================","\n")
print.table(x$phirmat,digits);
cat(" ","\n")
cat("Right p value","\n")
cat("=============","\n")
print.table(x$pvphirmatright,digits);
cat(" ","\n")
cat("Left p value","\n")
cat("============","\n")
print.table(x$pvphirmatleft,digits);
cat(" ","\n")
cat("Dyadic contribution to asymmetry","\n")
cat("================================","\n")
print.table(x$PHIij,digits)
cat(" ","\n");
cat("Right p value","\n")
cat("=============","\n")
print.table(x$pvPHIijright,digits);
cat(" ","\n")
cat("Left p value","\n")
cat("============","\n")
print.table(x$pvPHIijleft,digits);
invisible(x)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.