#  G1.R for Generalizability Theory analyses


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


# Specify the data to be analyzed. 
# The data matrix must be named "scores".


#***************  START OF TRIAL-RUN DATA  ************************************


# one-facet crossed design, as in p x i
# Type 1, based on Brennan, 2001, page 431, Table A.1
# data from Brennan, 2001, p. 28, Table 2.2, & p. 36, Table 2.4
# 10 persons, 12 items
scores <- t(matrix(c(
1,  0,  1,  0,  0,  0,  0,  0,  0,  0,  0,  0,
1,  1,  1,  0,  0,  1,  0,  0,  0,  0,  0,  0,
1,  1,  1,  1,  1,  0,  0,  0,  0,  0,  0,  0,
1,  1,  0,  1,  1,  0,  0,  1,  0,  0,  0,  0,
1,  1,  1,  1,  1,  0,  1,  0,  0,  0,  0,  0,
1,  1,  1,  0,  1,  1,  1,  0,  0,  0,  0,  0,
1,  1,  1,  1,  1,  1,  1,  0,  0,  0,  0,  0,
1,  1,  1,  1,  0,  1,  1,  1,  1,  1,  0,  0,
1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  0,
1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1   ), 12, 10) )


# one-facet nested design, as in i : p 
# Type 2, based on Brennan, 2001, page 431, Table A.2
# data from Brennan, 2001, p. 43, Table 2.8, & p. 44
# 10 persons, 8 items per person, but none of the items are the same, so there are 80 items.
scores <- t(matrix(c(
   2,   6,   7,   5,   2,   5,   5,   5,
   4,   5,   6,   7,   6,   7,   5,   7,
   5,   5,   4,   6,   5,   4,   5,   5,
   5,   9,   8,   6,   5,   7,   7,   6,
   4,   3,   5,   6,   4,   5,   6,   4,
   4,   4,   4,   7,   6,   4,   7,   8,
   2,   6,   6,   5,   2,   7,   7,   5,
   3,   4,   4,   5,   6,   6,   6,   4,
   0,   5,   4,   5,   5,   5,   5,   3,
   6,   8,   7,   6,   6,   8,   8,   6   ), 8, 10) )




# two-facet crossed design, as in  p x i x o
# Type 3, based on Brennan, 2001, page 432, Table A.3
# data from Brennan, 2001, p. 72, Table 3.1, & p. 111
# 10 persons, 4 items, 2 occasions.
scores <- t(matrix(c(
2, 6, 7, 5,  2, 5, 5, 5,
4, 5, 6, 7,  6, 7, 5, 7,
5, 5, 4, 6,  5, 4, 5, 5,
5, 9, 8, 6,  5, 7, 7, 6,
4, 3, 5, 6,  4, 5, 6, 4,
4, 4, 4, 7,  6, 4, 7, 8,
2, 6, 6, 5,  2, 7, 7, 5,
3, 4, 4, 5,  6, 6, 6, 4,
0, 5, 4, 5,  5, 5, 5, 3,
6, 8, 7, 6,  6, 8, 8, 6   ), 8, 10) )





# two-facet nested design, as in  p x (r : t)
# Type 4, based on Brennan, 2001, page 432, Table A.4
# data from Brennan, 2001, p. 73, Table 3.2, & p. 116
# 10 persons, 4 raters, 3 tasks.
scores <- t(matrix(c(
5, 6, 5, 5,  5, 3, 4, 5,  6, 7, 3, 3,
9, 3, 7, 7,  7, 5, 5, 5,  7, 7, 5, 2,
3, 4, 3, 3,  5, 3, 3, 5,  6, 5, 1, 6,
7, 5, 5, 3,  3, 1, 4, 3,  5, 3, 3, 5,
9, 2, 9, 7,  7, 7, 3, 7,  2, 7, 5, 3,
3, 4, 3, 5,  3, 3, 6, 3,  4, 5, 1, 2,
7, 3, 7, 7,  7, 5, 5, 7,  5, 5, 5, 4,
5, 8, 5, 7,  7, 5, 5, 4,  3, 2, 1, 1,
9, 9, 8, 8,  6, 6, 6, 5,  5, 8, 1, 1,
4, 4, 4, 3,  3, 5, 6, 5,  5, 7, 1, 1   ), 12, 10) )



#***************  END OF TRIAL-RUN DATA  ************************************



# Enter the number of levels/conditions of Facet 1 (e.g., # of items)
nfacet1 = 4

# Enter the number of levels/conditions of Facet 2 (e.g., # of occasions)
# You can ignore this step for single-facet designs
nfacet2 = 3

# For two-facet designs, Facet 1 is the facet with the fastest-changing conditions
# in the columns of your data matrix.  For example, if the first 10 columns/variables
# contained the data for 10 different items measured on occasion 1, and if the 
# next 10 columns/variables contained the data for the same 10 items measured on
# occasion 2, then items would be the fastest-changing facet.  As you slide from
# one column to the next across the data matrix, it is the item levels
# that change most quickly.  You would therefore enter a value of "10" for
# nfacet1 and a value of "2" for nfacet2 on the above statements

# Enter the design of your data on the "type =" statement below:
#    enter "1" for a single-facet fully-crossed design, as in P * F1
#    enter "2" for a single-facet nested nested design, as in F1 : P
#    enter "3" for a two-facet fully-crossed design, as in P * F1 * F2
#    enter "4" for a two-facet nested design, as in P * (F1 : F2)
#    enter "5" for a two-facet nested design, as in (F1 : P) * F2
#    enter "6" for a two-facet nested design, as in F1 : (P * F2)
#    enter "7" for a two-facet nested design, as in (F1 * F2) : P
#    enter "8" for a two-facet nested design, as in F1 : F2 : P
type = 3

# Enter D-study values for Facet 1, with a comma between the values
dfacet1 = c(4,  8,  16)
dfacet1 = c(12,6,4,3,2,2)
dfacet1 = c(5,10,15,20)

# Enter D-study values for Facet 2, with a comma between the values.  
# You can ignore this step for single-facet designs
dfacet2 = c( 2,  4)
dfacet2 = c(1,2,3,4,5,6)

# At the bottom of this file are commands for plotting
# 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:
#   enter "1" for relative error variances
#   enter "2" for absolute error variances
#   enter "3" for G-coefficients
#   enter "4" for phi coefficients
graphdat = 3


# End of user specifications


# Select All, and then Execute




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

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



cat('\n\n\nGENERALIZABILITY THEORY ANALYSES:\n')

if (type == 1) cat('\n\nDesign Type 1: single-facet fully-crossed design, as in P * F1\n\n')
if (type == 2) cat('\n\nDesign Type 2: single-facet nested nested design, as in F1 : P\n\n')
if (type == 3) cat('\n\nDesign Type 3: two-facet fully-crossed design, as in P * F1 * F2\n\n')
if (type == 4) cat('\n\nDesign Type 4: two-facet nested design, as in P * (F1 : F2)\n\n')
if (type == 5) cat('\n\nDesign Type 5: two-facet nested design, as in (F1 : P) * F2\n\n')
if (type == 6) cat('\n\nDesign Type 6: two-facet nested design, as in F1 : (P * F2)\n\n')
if (type == 7) cat('\n\nDesign Type 7: two-facet nested design, as in (F1 * F2) : P\n\n')
if (type == 8) cat('\n\nDesign Type 8: two-facet nested design, as in F1 : F2 : P\n\n')

p = nrow(scores)
i = nfacet1
h = nfacet2

cat('\nNumber of persons/objects ("P") = ', p);
cat('\n\nNumber of levels for Facet 1 ("F1") = ', i)
if (type > 2) cat('\n\nNumber of levels for Facet 2 ("F2") = ', h)


# single-facet fully-crossed design, as in P * F1
if (type == 1) {

# degrees of freedom
dfp = p-1
dfi = i-1
dfpi = dfp*dfi

# raw data computations needed for later sum of squares calculations
hABCD = scores
pmean = rowSums(hABCD)/i
p2mean = sum(pmean**2)
imean = colSums(hABCD)/p
i2mean = sum(imean**2)
pi2tot = sum(hABCD**2)
gmean = sum(imean)/i

Tpers = i*p2mean
Titem = p*i2mean
Tpi   = pi2tot
Tu    = p*i*(gmean*gmean)

# sum of squares
pss  = Tpers - Tu
iss  = Titem - Tu
piss = Tpi - Tpers - Titem + Tu

# mean squares
pms = pss/dfp
ims = iss/dfi
pims = piss/dfpi

# variance components: Brennan, 2001, p 79, 83, 85, 442
vcs = t(solve(matrix(c(i, 0, 1, 0, p, 1, 0, 0, 1),3,3,byrow=T)) %*% matrix(c(pms, ims, pims),3,1))
# setting negative variances to zero
if (any(vcs < 0)) {
vcs = ifelse(vcs < 0, 0, vcs)
cat('\n\nOne or more negative variance estimates have been set to zero') }

pvar  = vcs[1,1]
ivar  = vcs[1,2]
pivar = vcs[1,3]
totvar = sum(vcs)

# relative and absolute error variances
relv = pivar/i
absv = (ivar/i) + relv

# matrices for D-studies
persons  = c( 1, 0, 1 )
vfacet1  = c( 0, i, i )

cat('\n\n\nANOVA Table:\n\n')
anovatab = matrix(c(dfp,pss,pms,pvar,(pvar/totvar), dfi,iss,ims,ivar,(ivar/totvar), dfpi,piss,pims,pivar,(pivar/totvar)),3,5,byrow=TRUE)
colnames(anovatab) = c('df','SS','MS','Variance','%')
rownames(anovatab) = c('P','F1','P*F1')
print(round(anovatab,2))
}




# single-facet nested nested design, as in F1 : P
if (type == 2) {

# degrees of freedom
dfp = p-1
dfip = p*(i-1)

# raw data computations needed for later sum of squares calculations
hABCD = scores
pmean = rowSums(hABCD)/i
p2mean = sum(pmean**2)
pi2tot = sum(hABCD**2)
gmean = sum(pmean)/p

Tpers = i*p2mean
Tip   = pi2tot
Tu    = p*i*(gmean*gmean)

# sum of squares
pss = Tpers - Tu
ipss = Tip - Tpers

# mean squares
pms = pss/dfp
ipms = ipss/dfip

# variance components: Brennan, 2001, p 79, 83, 85, 442
vcs = t(solve(matrix(c(i, 1, 0, 1),2,2,byrow=T)) %*% matrix(c(pms, ipms),2,1))
# setting negative variances to zero
if (any(vcs < 0)) {
vcs = ifelse(vcs < 0, 0, vcs)
cat('\n\nOne or more negative variance estimates have been set to zero') }

pvar   = vcs[1,1]
ipvar  = vcs[1,2]
totvar = sum(vcs)

# relative and absolute error variances
relv = ipvar/i
absv = relv

# matrices for D-studies
persons  = c( 1, 1 )
vfacet1  = c( 0, i )

cat('\n\n\nANOVA Table:\n\n')
anovatab = matrix(c(dfp,pss,pms,pvar,(pvar/totvar),  dfip,ipss,ipms,ipvar,(ipvar/totvar) ),2,5,byrow=TRUE)
colnames(anovatab) = c('df','SS','MS','Variance','%')
rownames(anovatab) = c('P','F1:P')
print(round(anovatab,2))

}




# two-facet fully-crossed design, as in P * F1 * F2
if (type == 3) {

# degrees of freedom
dfp = nrow(scores) - 1
dfo = (h - 1)
dfi = (ncol(scores)/(h)) - 1
dfpo = (dfp*dfo)
dfpi = (dfp*dfi)
dfio = (dfi*dfo)
dfpoi = (dfp*dfo*dfi)

# pointers/counters
cntfirst = 1
cntlast = i

# raw data computations needed for later sum of squares calculations
pimean = matrix(0,p,i)
pmean = matrix(0,p,1)
omean = matrix(-9999,h,1)
pomean = matrix(-9999,p,h)
iomean = matrix(-9999,h,i)
for (cols in 1:h) {
hABCD = scores[,cntfirst:cntlast]
iomean[cols,] = colSums(hABCD)/p
pomean[,cols] = rowSums(hABCD)/i
pimean = pimean+hABCD
cntfirst = cntfirst + i
cntlast = cntlast + i
}

pmean = rowSums(pomean)/h
omean = colSums(pomean)/p
pimean = pimean/h
imean = colSums(pimean)/p
gmean = sum(pmean)/p

Tpers = (i*h)*sum(pmean**2)
Titem = (p*h)*sum((imean**2))
Tocca = (p*i)*sum(omean**2)
Tu  = p*(i)*(h)*((gmean)*(gmean))
Tpi = h*sum(sum(pimean**2))
Tpo = colSums((pomean**2))
Tpo = i*sum(Tpo)
Tio = p*sum(sum(iomean**2))
Tpoi = sum(sum(scores**2))

# sum of squares
pss = (Tpers-Tu)
iss = (Titem-Tu)
oss = (Tocca-Tu)
piss = (Tpi-Tpers-Titem+Tu)
poss = (Tpo-Tpers-Tocca+Tu)
ioss = (Tio-Titem-Tocca+Tu)
poiss = (Tpoi-Tpi-Tpo-Tio+Tpers+Titem+Tocca-Tu)

# mean squares
poims = poiss/dfpoi
ioms = ioss/dfio
poms = poss/dfpo
pims = piss/dfpi
oms = oss/dfo
ims = iss/dfi
pms = pss/dfp

# variance components: Brennan, 2001, p 79, 83, 85, 442
vcs = t(solve(matrix(c(i*h,0,0,h,i,0,1, 0,p*h,0,h,0,p,1, 0,0,p*i,0,i,p,1, 0,0,0,h,0,0,1,
          0,0,0,0,i,0,1,   0,0,0,0,0,p,1,   0,0,0,0,0,0,1),7,7,byrow=T))
           %*% matrix(c(pms, ims, oms, pims, poms, ioms, poims),7,1))
# setting negative variances to zero
if (any(vcs < 0)) {
vcs = ifelse(vcs < 0, 0, vcs)
cat('\n\nOne or more negative variance estimates have been set to zero') }

pvar  = vcs[1,1]
ivar  = vcs[1,2]
ovar  = vcs[1,3]
pivar = vcs[1,4]
povar = vcs[1,5]
iovar = vcs[1,6]
poivar = vcs[1,7]
totvar = sum(vcs)

# relative and absolute error variances
relv = (pivar/i)+(povar/h)+(poivar/(i*h))
absv = (ivar/i)+(ovar/h)+(iovar/(i*h))+relv

# matrices for D-studies
persons  = c( 1, 0, 0, 1, 1, 0, 1 )
vfacet1  = c( 0, i, 0, i, 0, i, i )
vfacet2  = c( 0, 0, h, 0, h, h, h )

cat('\n\n\nANOVA Table:\n\n')
anovatab = matrix(c(dfp,pss,pms,pvar,(pvar/totvar),   dfi,iss,ims,ivar,(ivar/totvar),   dfo,oss,oms,ovar,(ovar/totvar),
  dfpi,piss,pims,pivar,(pivar/totvar),  dfpo,poss,poms,povar,(povar/totvar),   dfio,ioss,ioms,iovar,(iovar/totvar), 
  dfpoi,poiss,poims,poivar,(poivar/totvar)  ),7,5,byrow=TRUE)
colnames(anovatab) = c('df','SS','MS','Variance','%')
rownames(anovatab) = c('P','F1', 'F2', 'P*F1', 'P*F2', 'F1*F2', 'P*F1*F2')
print(round(anovatab,2))

}



# two-facet nested design, as in P * (F1 : F2)
if (type == 4) {

# degrees of freedom
dfp = nrow(scores) - 1
dfh = (h-1)
dfih = (h)*(i-1)
dfph = dfp*dfh
dfpih = h*dfp*(i-1)
dftot = dfp+dfh+dfih+dfph+dfpih

# pointers/counters
cntfirst = 1
cntlast = i

# raw data computations needed for later sum of squares calculations
pmean = matrix(0,p,1)
hmean = matrix(-9999,h,1)
phmean = matrix(-9999,p,h)
ihmean = matrix(-9999,h,i)
for (cols in 1:h) {
hABCD = scores[, cntfirst:cntlast]
ihmean[cols,] = colSums(hABCD)/p
phmean[,cols] = rowSums(hABCD)/i
cntfirst = cntfirst + i
cntlast = cntlast + i
}

pmean = rowSums(phmean)/h
hmean = colSums(phmean)/p
gmean = sum(pmean)/p

Tp = (i*h)*(sum(pmean**2))
Th = (p*i)*sum((hmean**2))
Tih = p*sum(sum(ihmean**2))
Tph = i*sum(sum(phmean**2))
Tpih = sum(sum(scores**2))
Tu = p*i*h*gmean*gmean

# sum of squares
pss = Tp-Tu
hss = Th-Tu
ihss = Tih-Th
phss = Tph-Tp-Th+Tu
pihss = Tpih-Tph-Tih+Th
sstot = pss+hss+ihss+phss+pihss

# mean squares
pms = pss/dfp
hms = hss/dfh
ihms = ihss/dfih
phms = phss/dfph
pihms = pihss/dfpih

# variance components: Brennan, 2001, p 79, 83, 85, 442
vcs = t(solve(matrix(c( i*h,0,0,i,1, 0,p*i,p,i,1, 0,0,p,0,1, 0,0,0,i,1, 0,0,0,0,1),5,5,byrow=T)) 
    %*% matrix(c( pms, hms, ihms, phms, pihms),5,1))
# setting negative variances to zero
if (any(vcs < 0)) {
vcs = ifelse(vcs < 0, 0, vcs)
cat('\n\nOne or more negative variance estimates have been set to zero') }

pvar   = vcs[1,1]
hvar   = vcs[1,2]
ihvar  = vcs[1,3]
phvar  = vcs[1,4]
pihvar = vcs[1,5]
totvar = sum(vcs)

# relative and absolute error variances
relv = (pihvar/(i*h))+(phvar/h)
absv = (hvar/h)+(ihvar/(i*h)) + relv

# matrices for D-studies
persons  = c( 1, 0, 0, 1, 1 )
vfacet1  = c( 0, 0, i, 0, i )
vfacet2  = c( 0, h, h, h, h )

cat('\n\n\nANOVA Table:\n\n')
anovatab = matrix(c(dfp,pss,pms,pvar,(pvar/totvar),  dfh,hss,hms,hvar,(hvar/totvar),
   dfih,ihss,ihms,ihvar,(ihvar/totvar),  dfph,phss,phms,phvar,(phvar/totvar), 
   dfpih,pihss,pihms,pihvar,(pihvar/totvar)  ),5,5,byrow=TRUE)
colnames(anovatab) = c('df','SS','MS','Variance','%')
rownames(anovatab) = c('P','F2','F1:F2','P*F2','P*F1:F2')
print(round(anovatab,2))

}




# two-facet nested design, as in (F1 : P) * F2
if (type == 5) {

# degrees of freedom
dfp = nrow(scores) -1
dfh = h-1
dfip = p*(i-1)
dfph = (p-1)*(h-1)
dfhip = p*(i-1)*(h-1)
dftot = dfp+dfh+dfip+dfph+dfhip

# initializing pointers/counters
cntfirst = 1
cntlast = i

# raw data computations needed for later sum of squares calculations
pmean = matrix(0,p,1)
hmean = matrix(-9999,h,1)
phmean = matrix(-9999,p,h)
ipmean = matrix(0,p,i)
for (cols in 1:h) {
hABCD = scores[, cntfirst:cntlast]
phmean[,cols] = rowSums(hABCD)/i
ipmean = ipmean+ hABCD
cntfirst = cntfirst + i
cntlast = cntlast + i
}

pmean = rowSums(phmean)/h
hmean = colSums(phmean)/p
gmean = sum(pmean)/p
ipmean = ipmean/h

Tp = (i*h)*(sum(pmean**2))
Th = (p*i)*sum((hmean**2))
Tip = h*sum(sum(ipmean**2))
Tph = i*sum(sum(phmean**2))
Thip = sum(sum(scores**2))
Tu = p*i*h*gmean*gmean

# sum of squares
pss = Tp-Tu
hss = Th-Tu
ipss = Tip-Tp
phss = Tph-Tp-Th+Tu
hipss = Thip-Tph-Tip+Tp
sstot = pss+hss+ipss+phss+hipss

# mean squares
pms = pss/dfp
hms = hss/dfh
ipms = ipss/dfip
phms = phss/dfph
hipms = hipss/dfhip

# variance components: Brennan, 2001, p 79, 83, 85, 442
vcs = t(solve(matrix(c(i*h,0,0,i,1, 0,p*i,0,i,1, 0,0,h,0,1, 0,0,0,i,1, 0,0,0,0,1),5,5,byrow=T))
     %*% matrix(c(pms, hms, ipms, phms, hipms),5,1))
# setting negative variances to zero
if (any(vcs < 0)) {
vcs = ifelse(vcs < 0, 0, vcs)
cat('\n\nOne or more negative variance estimates have been set to zero') }

pvar   = vcs[1,1]
hvar   = vcs[1,2]
ipvar  = vcs[1,3]
phvar  = vcs[1,4]
hipvar = vcs[1,5]
totvar = sum(vcs)

# relative and absolute error variances
relv = (ipvar/i)+(phvar/h)+(hipvar/(i*h))
absv = (hvar/h)+relv

# matrices for D-studies
persons  = c( 1, 0, 1, 1, 1 )
vfacet1  = c( 0, 0, i, 0, i )
vfacet2  = c( 0, h, 0, h, h )

cat('\n\n\nANOVA Table:\n\n')
anovatab = matrix(c(dfp,pss,pms,pvar,(pvar/totvar),   dfh,hss,hms,hvar,(hvar/totvar), 
   dfip,ipss,ipms,ipvar,(ipvar/totvar),  dfph,phss,phms,phvar,(phvar/totvar),  
   dfhip,hipss,hipms,hipvar,(hipvar/totvar)  ),5,5,byrow=TRUE)
colnames(anovatab) = c('df','SS','MS','Variance','%')
rownames(anovatab) = c('P','"F2','F1:P','P*F2','F2*F1:P')
print(round(anovatab,2))

}




# two-facet nested design, as in F1 : (P * F2)
if (type == 6) {

# degrees of freedom
dfp = nrow(scores) -1
dfh = h-1
dfph = (p-1)*(h-1)
dfiph = p*h*(i-1)
dftot = dfp+dfh+dfph+dfiph

# pointers/counters
cntfirst = 1
cntlast = i

# raw data computations needed for later sum of squares calculations
pmean = matrix(0,p,1)
hmean = matrix(-9999,h,1)
phmean = matrix(-9999,p,h)
for (cols in 1:h) {
hABCD = scores[, cntfirst:cntlast]
phmean[,cols] = rowSums(hABCD)/i
cntfirst = cntfirst + i
cntlast = cntlast + i
}

pmean = rowSums(phmean)/h
hmean = colSums(phmean)/p
gmean = sum(pmean)/p

Tp = (i*h)*(sum(pmean**2))
Th = (p*i)*sum((hmean**2))
Tph = i*sum(sum(phmean**2))
Tiph = sum(sum(scores**2))
Tu = p*i*h*gmean*gmean

# sum of squares
pss = Tp-Tu
hss = Th-Tu
phss = Tph-Tp-Th+Tu
iphss = Tiph-Tph
sstot = pss+hss+phss+iphss

# mean squares
pms = pss/dfp
hms = hss/dfh
phms = phss/dfph
iphms = iphss/dfiph

# variance components: Brennan, 2001, p 79, 83, 85, 442
vcs = t(solve(matrix(c(i*h,0,i,1, 0,p*i,i,1, 0,0,i,1, 0,0,0,1),4,4,byrow=T))
     %*% matrix(c(pms, hms, phms, iphms),4,1))
# setting negative variances to zero
if (any(vcs < 0)) {
vcs = ifelse(vcs < 0, 0, vcs)
cat('\n\nOne or more negative variance estimates have been set to zero') }

pvar   = vcs[1,1]
hvar   = vcs[1,2]
phvar  = vcs[1,3]
iphvar = vcs[1,4]
totvar = sum(vcs)

# relative and absolute error variances & G & phi coefficients
relv = (phvar/h)+(iphvar/(i*h))
absv = (hvar/h)+relv

# matrices for D-studies
persons  = c( 1, 0, 1, 1 )
vfacet1  = c( 0, 0, 0, i )
vfacet2  = c( 0, h, h, h )

cat('\n\n\nANOVA Table:\n\n')
anovatab = matrix(c(dfp,pss,pms,pvar,(pvar/totvar),  dfh,hss,hms,hvar,(hvar/totvar),
   dfph,phss,phms,phvar,(phvar/totvar),   dfiph,iphss,iphms,iphvar,(iphvar/totvar) ),4,5,byrow=TRUE)
colnames(anovatab) = c('df','SS','MS','Variance','%')
rownames(anovatab) = c('P','F2','P*F2','F1:P*F2')
print(round(anovatab,2))

}




# two-facet nested design, as in (F1 * F2) : P
if (type == 7) {

# degrees of freedom
dfp = nrow(scores) -1
dfip = p*(i-1)
dfhp = p*(h-1)
dfihp = p*(i-1)*(h-1)
dftot = dfp+dfip+dfhp+dfihp

# pointers/counters
cntfirst = 1
cntlast = i

# raw data computations needed for later sum of squares calculations
pmean = matrix(0,p,1)
phmean = matrix(-9999,p,h)
ipmean = matrix(0,p,i)
hpmean = matrix(-9999,p,h)
for (cols in 1:h) {
hABCD = scores[, cntfirst:cntlast]
phmean[,cols] = rowSums(hABCD)/i
ipmean = ipmean+ hABCD
hpmean[,cols] = rowSums(hABCD)/i
cntfirst = cntfirst + i
cntlast = cntlast + i
}

pmean = rowSums(hpmean)/h
hmean = colSums(phmean)/p
gmean = sum(pmean)/p
ipmean = ipmean/h

Tp = (i*h)*(sum(pmean**2))
Th = (p*i)*sum(sum(hmean**2))
Tip = h*sum(sum(ipmean**2))
Thp = i*sum(sum(hpmean**2))
Tihp = sum(sum(scores**2))
Tu = p*i*h*gmean*gmean

# sum of squares
pss = Tp-Tu
ipss = Tip-Tp
hpss = Thp-Tp
ihpss = Tihp-Tip-Thp+Tp
sstot = pss+ipss+hpss+ihpss

# mean squares
pms = pss/dfp
ipms = ipss/dfip
hpms = hpss/dfhp
ihpms = ihpss/dfihp

# variance components: Brennan, 2001, p 79, 83, 85, 442
vcs = t(solve(matrix(c(i*h,h,i,1, 0,h,0,1, 0,0,i,1, 0,0,0,1),4,4,byrow=T)) 
      %*% matrix(c(pms, ipms, hpms, ihpms),4,1))
# setting negative variances to zero
if (any(vcs < 0)) {
vcs = ifelse(vcs < 0, 0, vcs)
cat('\n\nOne or more negative variance estimates have been set to zero') }

pvar   = vcs[1,1]
ipvar  = vcs[1,2]
hpvar  = vcs[1,3]
ihpvar = vcs[1,4]
totvar = sum(vcs)

# relative and absolute error variances
relv = (ipvar/i)+(hpvar/h)+(ihpvar/(i*h))
absv = relv

# matrices for D-studies
persons  = c( 1, 1, 1, 1 )
vfacet1  = c( 0, i, 0, i )
vfacet2  = c( 0, 0, h, h )

cat('\n\n\nANOVA Table:\n\n')
anovatab = matrix(c(dfp,pss,pms,pvar,(pvar/totvar),  dfip,ipss,ipms,ipvar,(ipvar/totvar),
     dfhp,hpss,hpms,hpvar,(hpvar/totvar),  dfihp,ihpss,ihpms,ihpvar,(ihpvar/totvar)),4,5,byrow=TRUE)
colnames(anovatab) = c('df','SS','MS','Variance','%')
rownames(anovatab) = c('P','F1*P','F2*P','F1:F2:P')
print(round(anovatab,2))

}




# two-facet nested design, as in F1 : F2 : P
if (type == 8) {

# degrees of freedom
dfp = nrow(scores) -1
dfhp = p*(h-1)
dfihp = p*h*(i-1)
dftot = dfp+dfhp+dfihp

# pointers/counters
cntfirst = 1
cntlast = i

# raw data computations needed for later sum of squares calculations
pmean = matrix(0,p,1)
phmean = matrix(-9999,p,h)
hpmean = matrix(-9999,p,h)
for (cols in 1:h) {
hABCD = scores[, cntfirst:cntlast]
phmean[,cols] = rowSums(hABCD)/i
hpmean[,cols] = rowSums(hABCD)/i
cntfirst = cntfirst + i
cntlast = cntlast + i
}

pmean = rowSums(hpmean)/h
hmean = colSums(phmean)/p
gmean = sum(pmean)/p

Tp = (i*h)*(sum(pmean**2))
Th = (p*i)*sum((hmean**2))
Thp = i*sum(sum(hpmean**2))
Tihp = sum(sum(scores**2))
Tu = p*i*h*gmean*gmean

# sum of squares
pss = Tp-Tu
hpss = Thp-Tp
ihpss = Tihp-Thp
sstot = pss+hpss+ihpss

# mean squares
pms = pss/dfp
hpms = hpss/dfhp
ihpms = ihpss/dfihp

# variance components: Brennan, 2001, p 79, 83, 85, 442
vcs = t(solve(matrix(c(i*h, i, 1, 0, i, 1, 0, 0, 1),3,3,byrow=T)) 
    %*% matrix(c(pms, hpms, ihpms),3,1))
# setting negative variances to zero
if (any(vcs < 0)) {
vcs = ifelse(vcs < 0, 0, vcs)
cat('\n\nOne or more negative variance estimates have been set to zero') }

pvar   = vcs[1,1]
hpvar  = vcs[1,2]
ihpvar = vcs[1,3]
totvar = sum(vcs)

# relative and absolute error variances
relv =  (hpvar/h)+(ihpvar/(i*h))
absv = relv

# matrices for D-studies
persons  = c( 1, 1, 1 )
vfacet1  = c( 0, 0, i )
vfacet2  = c( 0, h, h )

cat('\n\n\nANOVA Table:\n\n')
anovatab = matrix(c(dfp,pss,pms,pvar,(pvar/totvar),  
   dfhp,hpss,hpms,hpvar,(hpvar/totvar),  dfihp,ihpss,ihpms,ihpvar,(ihpvar/totvar)  ),3,5,byrow=TRUE)
colnames(anovatab) = c('df','SS','MS','Variance','%')
rownames(anovatab) = c('P','F2:P','F1:F2:P')
print(round(anovatab,2))

}



# G & phi coefficients
relG = pvar/(pvar+relv)
absG = pvar/(pvar+absv)

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)

nfacets = 2
if (type <= 2) nfacets = 1

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")
print(round(dabsv,3))

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

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

cat("\n\nD-Study Phi Coefficients\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 )
}

}

gtheory1( scores, nfacets, persons, vfacet1, vfacet2, dfacet1, dfacet2, graphdat )