#  G2.R for Generalizability Theory analyses



# Follow the instructions until the "# End of user specifications" line,
# then simply Select All, and then Execute


# Enter a one-row data matrix of variance components with a comma 
# between each component. The data matrix must be named "vcs".

# data from Brennan, 2001, Data Set 3,  p 72 & p 111
vcs =  c(0.5528,  0.4417, 0.0074,  0.5750,  0.1009, 0.1565, 0.9352 )


# How many facets are in the design? Enter "1" or "2"
nfacets = 2

# Look at the vcs values and then indicate, in the "persons" matrix below, whether
# persons were involved in the computations for each term. Enter "1" if persons
# were involved, and enter "0" if persons were not involved.
# e.g., if there were 7 variance components for 
# p i o p*i p*o i*o p*i*o, then enter
# 1,0,0, 1,  1,  0,   1
persons = c(  1,  0,  0,  1,  1,  0,  1 )


# Look at the vcs values and then indicate, in the "vfacet1" matrix below, 
# whether Facet 1 was involved in the computations for each term. 
# Enter the number of levels of Facet 1 if Facet 1 was involved, and 
# enter "0" if Facet 1 was not involved in the computations for each term.
# e.g., if there were 7 variance components, and Facet 1 had 4 levels 
# p i o p*i p*o i*o p*i*o, then enter
# 0,4,0, 4,  0,  4,   4
vfacet1  = c( 0, 4, 0, 4, 0, 4, 4 )


# Look at the vcs values and then indicate, in the "vfacet2" matrix below, 
# whether Facet 2 was involved in the computations for each term. Enter the number of
# levels of Facet 2 if Facet 2 was involved, and enter "0" if Facet 2 was not 
# involved in the computations for each term.
# e.g., if there were 7 variance components, and Facet 2 had 2 levels 
# p i o p*i p*o i*o p*i*o, then enter
# 0,0,2, 0,  2,  2,   2;  
# You can ignore this step for single-facet designs
vfacet2 = c( 0,  0,  2,  0,  2,  2,  2  )


# Enter D-study values for Facet 1; enter the values inside curly
# brackets, and place a comma between the values
dfacet1 = c(4,  8,  16)


# Enter D-study values for Facet 2; enter the values inside curly
# brackets, and place a comma between the values.  You can ignore this
# step for single-facet designs
dfacet2 = c( 2,  4)


# The program produces a plot of the results for the D-study values that you specified above.  
# Specify the data that you would like to plot by entering the appropriate number on the 
# graphdat statement below:
  # enter "1" for relative error variances;
  # enter "2" for absolute error variances;
  # enter "3" for G-coefficients;
  # enter "4" for phi coefficients
graphdat = 1	
	
	
# End of user specifications

# Select All, and then Execute




##########################################################################################

gtheory2 <- function(vcs, nfacets, persons, vfacet1, vfacet2, dfacet1, dfacet2, graphdat ) {



vcs = matrix(vcs,1,length(vcs))

# setting negative variance components to zero
for (lupe in 1:ncol(vcs)) { if (vcs[1,lupe] < 0) { vcs[1,lupe] = 0 }}

if (nfacets == 1) { dfacet2 = 1 }

# absolute & relative error variances
absv = 0
relv = 0
for (lupe in 2:length(vfacet1)) {
f2dum = 1
f1dum = 1
if ( vfacet1[lupe] != 0 ) { f1dum = vfacet1[lupe] }
if ( vfacet2[lupe] != 0 ) { f2dum = vfacet2[lupe] }
absv = absv + ( vcs[1,lupe] / (f1dum * f2dum ) )
if ( persons[lupe] != 0 ) { relv = relv + ( vcs[1,lupe] / (f1dum * f2dum) ) }
}

relG = vcs[1,1] / ( vcs[1,1] + relv )
absG = vcs[1,1] / ( vcs[1,1] + absv )

percents = vcs / sum(vcs)

cat('\n\n\nGeneralizability Theory Analyses\n')

cat("\n\nEntered Specifications, Variance Components and\n")
cat("their Proportions of the Total Variance:\n\n")
 
vcomps <- cbind( persons, vfacet1, vfacet2, t(vcs) , t(percents) )
dimnames(vcomps) <-list(rep("", dim(vcomps)[1]))
colnames(vcomps) <- c( "Persons", "Facet1", "Facet2", "VarComp", "Proport" )
print(round(vcomps,3))

cat('\n\nRelative Error Variance = ', round(relv,3))
cat('\n\nAbsolute Error Variance = ', round(absv,3))
cat('\n\nG coefficient = ', round(relG,3))
cat('\n\nPhi coefficient = ', round(absG,3))



# D study

cat("\n\n\nD-Study:\n")
cat("\nEntered D-Study values for Facet 1:", dfacet1)
cat("\nEntered D-Study values for Facet 2:", dfacet2)

dabsv  = -9999
drelv  = dabsv
drelG  = dabsv
dabsG  = dabsv
# dabsv2  = c( -9999, -9999, -9999 )
# drelv2  = dabsv2
# drelG2  = dabsv2
# dabsG2  = dabsv2

if (nfacets == 1) { dfacet2 = 1; vfacet2 = matrix(1,1,ncol(vcs)) }

for (f2 in 1:length(dfacet2)) {
resdum = matrix(-9999,length(dfacet1),4)
counter = 1

for (f1 in 1:length(dfacet1)) {
absvdum = 0
relvdum = 0

for (lupe in 2:ncol(vcs)) {
f2dum = 1
f1dum = 1
if ( vfacet1[lupe] > 1 ) { f1dum = dfacet1[f1] }
if ( vfacet2[lupe] > 1 ) { f2dum = dfacet2[f2] }
absvdum = absvdum + ( vcs[1,lupe] / (f1dum * f2dum ) )
if ( persons[lupe] != 0 ) {
relvdum = relvdum + ( vcs[1,lupe] / (f1dum * f2dum ) )
}
}
gdum = vcs[1,1] / ( vcs[1,1] + relvdum )
pdum = vcs[1,1] / ( vcs[1,1] + absvdum )

resdum[counter,] = c( absvdum, relvdum, gdum, pdum )

counter = counter + 1
}

dabsv = c( dabsv , resdum[,1] )
drelv = c( drelv , resdum[,2] )
drelG = c( drelG , resdum[,3] )
dabsG = c( dabsG , resdum[,4] )
}

dabsv = dabsv[-1]
drelv = drelv[-1]
drelG = drelG[-1]
dabsG = dabsG[-1]

dabsv = matrix(dabsv,length(dfacet1),length(dfacet2))
colnames(dabsv) = dfacet2; rownames(dabsv) = dfacet1
drelv = matrix(drelv,length(dfacet1),length(dfacet2))
colnames(drelv) = dfacet2; rownames(drelv) = dfacet1
drelG = matrix(drelG,length(dfacet1),length(dfacet2))
colnames(drelG) = dfacet2; rownames(drelG) = dfacet1
dabsG = matrix(dabsG,length(dfacet1),length(dfacet2))
colnames(dabsG) = dfacet2; rownames(dabsG) = dfacet1

cat("\n\nIn the D-study results below, the levels of Facet 1 appear in")
cat("\nthe first column, and the levels of Facet 2 appear in the first row.\n")

cat("\n\nD-Study Absolute Error Variances:\n\n")
print(round(dabsv,3))

cat("\n\nD-Study Relative Error Variances\n\n")
print(round(drelv,3))

cat("\n\nD-Study G Coefficients\n\n")
print(round(drelG,3))

cat("\n\nD-Study Phi Coefficients\n\n")
print(round(dabsG,3))


par(ann=F, font.main=1, font.lab=1, font.axis=1, cex=1, cex.main=1, cex.lab=1, cex.axis=1,
    lwd=2, las=1, mai=c(1.02,0.82,0.82,0.42), pty="s" ) 


if ( graphdat == 1) { 
matplot(dfacet1, drelv, type = c("b"), lty=1, xlim = c( min(dfacet1), max(dfacet1) ) )
title(main='D-Study: Relative Error Variances', ylab='Relative Error Variance', xlab='Facet 1')
legend("topright", legend = dfacet2, title='Facet 2:', lwd=2, col=dfacet2 )
}
if ( graphdat == 2) { 
matplot(dfacet1, dabsv, type = c("b"), lty=1, xlim = c( min(dfacet1), max(dfacet1) ) )
title(main='D-Study: Absolute Error Variances', ylab='Absolute Error Variance', xlab='Facet 1')
legend("topright", legend = dfacet2, title='Facet 2:', lwd=2, col=dfacet2 )
}
if ( graphdat == 3) { 
matplot(dfacet1, drelG, type = c("b"), lty=1, ylim = c( 0, 1 ), xlim = c( min(dfacet1), max(dfacet1)) )
title(main='D-Study: G Coefficients', ylab='G Coefficient', xlab='Facet 1')
legend("bottomright", legend = dfacet2, title='Facet 2:', lwd=2, col=dfacet2 )
}
if ( graphdat == 4) { 
matplot(dfacet1, dabsG, type = c("b"), lty=1, ylim = c( 0, 1 ), xlim = c( min(dfacet1), max(dfacet1)) )
title(main='D-Study: Phi Coefficients', ylab='Phi Coefficient', xlab='Facet 1')
legend("bottomright", legend = dfacet2, title='Facet 2:', lwd=2, col=dfacet2 )
}


}

gtheory2( vcs, nfacets, persons, vfacet1, vfacet2, dfacet1, dfacet2, graphdat )