orthpoly<-function(w,x=0:(length(w)-1),qstart=1)
{
# orthogonal polynomial values using weights w
#
w<-w/sum(w)
n<-length(w)
p<-matrix(0,nrow=n+1,ncol=n)
A<-matrix(0,nrow=n,ncol=1)
B<-matrix(0,nrow=n,ncol=1)
q<-matrix(0,nrow=n,ncol=n)
p[(-1)+2,]<-0 # p_{-1}=0
q[(0)+1,]<-qstart # q_{0}=1 (default)
A[(0)+1]<-sqrt(sum(w*q[(0)+1,]^2))
p[(0)+2,]<-q[(0)+1,]/A[(0)+1]
for (j in 1:(n-1))
{
B[(j)+1]<-sum(x*w*q[(j-1)+1,]^2)/sum(w*q[(j-1)+1,]^2)
q[(j)+1,]<-(x-B[(j)+1])*p[(j-1)+2,]-A[(j-1)+1]*p[(j-2)+2,]
A[(j)+1]<-sqrt(sum(w*q[(j)+1,]^2))
p[(j)+2,]<-q[(j)+1,]/A[(j)+1]
}
return(p[((0)+2):((n-1)+2),])
}
orthpoly(rep(1/6,6))
#           [,1]       [,2]       [,3]       [,4]       [,5]      [,6]
#[1,]  1.0000000  1.0000000  1.0000000  1.0000000  1.0000000 1.0000000
#[2,] -1.4638501 -0.8783101 -0.2927700  0.2927700  0.8783101 1.4638501
#[3,]  1.3363062 -0.2672612 -1.0690450 -1.0690450 -0.2672612 1.3363062
#[4,] -0.9128709  1.2780193  0.7302967 -0.7302967 -1.2780193 0.9128709
#[5,]  0.4629100 -1.3887301  0.9258201  0.9258201 -1.3887301 0.4629100
#[6,] -0.1543033  0.7715167 -1.5430335  1.5430335 -0.7715167 0.1543033



doit<-function(Nij,mult,SSeffects.max=4,SSeffects.min=1,rankscores=T,residualSS.pval=0.15,SS.pval=0.10)
{
t<-nrow(Nij)
k<-ncol(Nij)
ni<-apply(Nij,1,sum)
nj<-apply(Nij,2,sum)
n<-sum(ni)
pj<-nj/n
if (rankscores) xj<-cumsum(c(0,nj)[-(k+1)])+(nj+1)/2
else xj<-1:k
G<-orthpoly(pj,xj)
#
tCri<-Nij%*%t(G[-1,])/sqrt(ni)*sqrt(mult)
colnames(tCri)<-1:ncol(tCri)
for (i in 1:ncol(tCri)) colnames(tCri)[i]<-paste(i,"-th order",sep="")
colnames(tCri)[1:min(4,ncol(tCri))]<-c("Location","Dispersion","Skewness","Kurtosis")[1:min(4,ncol(tCri))]
#
SS<-apply(tCri^2,2,sum)
df<-rep(t-1,k-1)
#choose i so that no individual pvals<SS.pval are hidden in the residual, and residual pval>residualSS.pval
for (i in 1:(k-1)) if (1-pchisq(sum(SS[i:(k-1)]),sum(df[i:(k-1)]))>residualSS.pval) break
i<-max(i,max(which(1-pchisq(SS,df)<SS.pval))+1)
i<-min(max(i,SSeffects.min+1),SSeffects.max+1,k-1)
table<-cbind(c(df[1:(i-1)],sum(df[i:(k-1)])),c(SS[1:(i-1)],sum(SS[i:(k-1)])),c(1-pchisq(SS[1:(i-1)],df[1:(i-1)]),1-pchisq(sum(SS[i:(k-1)]),sum(df[i:(k-1)]))))
table<-rbind(table,c(sum(df),sum(SS),1-pchisq(sum(SS),sum(df))))
rownames(table)[nrow(table)]<-"Total"
rownames(table)[nrow(table)-1]<-"Residual"
colnames(table)<-c("df","SS","pvalue")
#
out<-list(tCri,table)
names(out)<-list("Cri","partition")
return(out)
}

CRD<-function(Nij,SSeffects.max=4,SSeffects.min=1,rankscores=T,residualSS.pval=0.15,SS.pval=0.10)
{
SSeffects.min<-min(SSeffects.min,SSeffects.max)
SSeffects.max<-max(SSeffects.min,SSeffects.max)
n<-sum(sum(Nij))
#
return(doit(Nij,(n-1)/n,SSeffects.max,SSeffects.min,rankscores,residualSS.pval,SS.pval))
}

RBD<-function(Nij,SSeffects.max=4,SSeffects.min=1,rankscores=T,residualSS.pval=0.15,SS.pval=0.10)
{
SSeffects.min<-min(SSeffects.min,SSeffects.max)
SSeffects.max<-max(SSeffects.min,SSeffects.max)
t<-nrow(Nij)
#
return(doit(Nij,(t-1)/t,SSeffects.max,SSeffects.min,rankscores,residualSS.pval,SS.pval))
}


ones<-function(r)
{
matrix(rep(1,r^2),ncol=r)
}

xtotab<-function(x)
{
 # Function to transform a (n-by-p) consumer-by-product data matrix
 # of ordered categorical data using category labels 1..max(x) into 
 # a (p-by-b) product-by-rank matrix counting "fractionals" 
n<-length(x[,1]) # Number of blocks #
p<-length(x[1,])  # Number of treatments # 
b<-length(x[1,][!is.na(x[1,])]) # Block size #
rx<-t(apply(x,1,rank))
y<-x
s<-max(x[!is.na(x)])
for (i in 1:p) for (j in 1:n) for (k in 1:s) 
y[j,][x[j,]==k]<-1/length((x[j,][!is.na(x[j,])]==k)[(x[j,][!is.na(x[j,])]==k)==T])
res<-matrix(rep(0,p*b),ncol=b)
for (j in 1:n) for (i in 1:p) if (!is.na(y[j,i]))
{
for (k in 1:b){
if (y[j,i]==1/k) { 
res[i,rx[j,i]+(1:k)-mean(1:k)]<-res[i,rx[j,i]+(1:k)-mean(1:k)]+1/k
}
}
}	
rownames(res)<-colnames(x)
colnames(res)<-1:ncol(res)
return(res)
}


U<-function(x)
{
n<-length(x[,1]) # Number of blocks #
p<-length(x[1,])  # Number of treatments # 
b<-length(x[1,][!is.na(x[1,])]) # Block size #
rx<-t(apply(x,1,rank))
y<-x
s<-max(x[!is.na(x)])
for (i in 1:p) for (j in 1:n) for (k in 1:s) 
y[j,][x[j,]==k]<-1/length((x[j,][!is.na(x[j,])]==k)[(x[j,][!is.na(x[j,])]==k)==T])
yy<-matrix(rep(0,n*b),ncol=b)
for (j in 1:n) for (i in 1:b)
yy[j,i]<-y[j,][order(rx[j,])[i]]
yyy<-t(apply(rx,1,sort))[,1:b]
Vs<-diag(apply(yy,2,sum))
for (i in 1:(b-1)) for (ii in (i+1):b){
ind<-yyy[,i]==yyy[,ii]
Vs[i,ii]<-sum(yy[ind,i])
Vs[ii,i]<-Vs[i,ii]
}
rownames(Vs)<-1:nrow(Vs)
colnames(Vs)<-1:ncol(Vs)
return(Vs)
}

ortho<-function(p,x)
{
k<-length(x)
g<-matrix(rep(0,k*(k-1)),ncol=k)
mu<-sum(x*p)
mu2<-sum(p*(x-mu)**2)
mu3<-sum(p*(x-mu)**3)
mu4<-sum(p*(x-mu)**4)
a<-(mu4-mu3*mu3/mu2-mu2**2)**(-0.5)
g[1,]<-(x-mu)/sqrt(mu2)
g[2,]<-a*((x-mu)**2-mu3*(x-mu)/mu2-mu2)
if (k>3) for (j in 3:(k-1)){
B<-sum(p*x*g[j-1,]**2)
C<-sum(p*x*g[j-1,]*g[j-2,])
A<-1/sqrt(sum(p*(x**2)*g[j-1,]**2)-B**2-C**2)
g[j,]<-A*((x-B)*g[j-1,]-C*g[j-2,])
}
g
}

or2U<-function(U)
{
r<-length(U[1,])
n<-sum(U[1,])
or1<-ortho(rep(1/r,r),1:r)[1,]
or1<-or1*sqrt(r*n/(or1%*%U%*%t(t(or1))))
w<-or1%*%U
B<--sum(((1:r)**2)*w)/sum((1:r)*w)
A<--B*(r+1)/2-sum((1:r)**2)/r
cc<-sqrt(r*n/(A*A*r*n+A*B*r*n*(r+1)+2*A*n*sum((1:r)**2)+B*B*(1:r)%*%U%*%t(t(1:r))+2*B*(1:r)%*%U%*%t(t(1:r)**2)+((1:r)**2)%*%U%*%t(t(1:r)**2)))
a<-A*cc
b<-B*cc
or2<-a+b*(1:r)+cc*(1:r)**2
matrix(c(or1,or2),ncol=r,byrow=T)
}

ascov2U<-function(U,p)
{
n<-sum(U[,1])
b<-length(U[1,])
kronecker(diag(rep(1,p))-ones(p)/p,U-n*ones(b)/b)/(p-1)
}

BlockAnalysis<-function(x)
{
# x should be a blocks-by-treatment data matrix #
# of ordered categorical data using category labels 1..max(x) #
# and with missing values (NAs) for 
n<-length(x[,1]) # Number of blocks #
p<-length(x[1,])  # Number of treatments # 
b<-length(x[1,][!is.na(x[1,])]) # Block size #
tab<-xtotab(x)    # Makes the treatment-by-rank matrix #
cz<-ascov2U(U(x),p) # U(x) constructs the U-matrix, ascov2U = cc times #
es<-eigen(cz,symmetric=T)                           # the Sigma-matrix #
evls<-es$values
evecs<-es$vectors   # The H (divided by root-cc) from proof of Prop2 is constructed #
evls2<-evls[evls>0.0001]  # to make the general version of the statistic#
dftot<-length(evls2)       # Finds the total degrees of freedom #
H<-evecs[,1:dftot]%*%diag(1/sqrt(evls2))
forv<-matrix(rep(n/p,b*p),ncol=b)
X<-t(c(t(tab-forv)))
chi2tot<-sum((X%*%H)^2)
or<-or2U(U(x))   # The generalized polynomials are constructed #
Z<-t(c(t((tab-forv)/sqrt(n/p))))
G<-kronecker(diag(rep(1,p)),t(t(or[1,])))/sqrt(b)
V1<-sqrt((p-1)/p)*t(G)%*%t(Z)
G<-kronecker(diag(rep(1,p)),t(t(or[2,])))/sqrt(b)
V2<-sqrt((p-1)/p)*t(G)%*%t(Z)
chi2lin<-sum(V1^2)
chi2sq<-sum(V2^2)
chi2s<-c(chi2tot,chi2lin,chi2sq,chi2tot-chi2lin-chi2sq)
dfs<-c(dftot,p-1,p-1,dftot-2*p+2)
pvals<-rep(NA,4)
pvals[dfs>0]<-1-pchisq(chi2s[dfs>0],dfs[dfs>0])
#list(matrix(c(chi2s,dfs,pvals),ncol=3),matrix(c(V1,V2),ncol=2,byrow=F))
#last bit changed to align output labels
out1<-matrix(c(V1,V2),ncol=2,byrow=F)
out2<-matrix(c(chi2s,dfs,pvals),ncol=3)
rownames(out1)<-colnames(x)
colnames(out1)<-1:ncol(out1)
for (i in 1:ncol(out1)) colnames(out1)[i]<-paste(i,"-th order",sep="")
colnames(out1)[1:min(4,ncol(out1))]<-c("Location","Dispersion","Skewness","Kurtosis")[1:min(4,ncol(out1))]
out2<-out2[c(2:nrow(out2),1),]
out2<-out2[,c(2,1,3:ncol(out2))]
rownames(out2)<-1:nrow(out2)
rownames(out2)[1:(nrow(out2)-2)]<-c("Location","Dispersion","Skewness","Kurotsis")[1:(nrow(out2)-2)]
rownames(out2)[nrow(out2)]<-"Total"
rownames(out2)[nrow(out2)-1]<-"Residual"
colnames(out2)<-c("df","SS","pvalue")
out<-list(out1,out2)
names(out)<-list("Cri","partition")
return(out)
}

xtotab2<-function(x)
{
n<-nrow(x)
rx<-apply(t(x),1,rank)
xi<-sort(unique(rx[,2]))
yj<-sort(unique(rx[,1]))
res<-matrix(0,length(yj),length(xi))
for (i in 1:n) 
for (k1 in 1:length(yj))
for (k2 in 1:length(xi))
if ((rx[i,1]==yj[k1])&(rx[i,2]==xi[k2])) res[k1,k2]<-res[k1,k2]+1
rownames(res)<-yj
colnames(res)<-xi
res
}

CorrelationAnalysis<-function(x)
{
n<-nrow(x)
rx<-apply(t(x),1,rank)
xi<-sort(unique(rx[,2]))
yj<-sort(unique(rx[,1]))
res<-xtotab2(x)
pxi<-apply(res,2,sum)/sum(res)
pyj<-apply(res,1,sum)/sum(res)
corrs<-matrix(0,length(yj)-1,length(xi)-1)
for (i in 1:(length(xi)-1))
for (j in 1:(length(yj)-1))
{
poly1<-t(as.matrix(ortho(pxi,xi)[i,]))
poly2<-t(as.matrix(ortho(pyj,yj)[j,]))
corrs[j,i]<-poly2%*%(res/sum(res))%*%t(poly1)
}
out<-list(corrs,2*pnorm(-abs(corrs*sqrt(sum(res)))))
names(out)<-list("correlations","pvalues")
return(out)
}


xtotab3<-function(x)
{
t<-ncol(x)
rx<-t(apply(x,1,rank))
yj<-sort(unique(c(rx)))
tab<-matrix(0,t,length(yj))
rownames(tab)<-colnames(x)
colnames(tab)<-yj
for (i in 1:length(yj))
for (j in 1:t)
{
tab[j,i]<-sum(yj[i]==rx[,j])
}
return(tab)
}

OrderedBlockAnalysis<-function(x)
{
n<-nrow(x)
t<-ncol(x)
rx<-t(apply(x,1,rank))
yj<-sort(unique(c(rx)))
tab<-matrix(0,length(yj),t)
colnames(tab)<-colnames(x)
rownames(tab)<-yj
for (i in 1:length(yj))
for (j in 1:t)
{
tab[i,j]<-sum(yj[i]==rx[,j])
}
xi<-1:t
pxi<-apply(tab,2,sum)/sum(tab)
pyj<-apply(tab,1,sum)/sum(tab)
corrs<-matrix(0,length(xi)-1,length(yj)-1)
for (i in 1:(length(xi)-1))
for (j in 1:(length(yj)-1))
{
poly1<-t(as.matrix(ortho(pxi,xi)[i,]))
poly2<-t(as.matrix(ortho(pyj,yj)[j,]))
corrs[i,j]<-poly2%*%(tab/sqrt(sum(tab)))%*%t(poly1)
}
Z<-corrs*sqrt((t-1)/t)
table<-cbind(c(1,1,(t-1)^2-2),c(Z[1,1]^2,Z[2,1]^2,sum(Z^2)-Z[1,1]^2-Z[2,1]^2))
table<-rbind(table,apply(table,2,sum))
table<-cbind(table,1-pchisq(table[,2],table[,1]))
rownames(table)<-c("Page","Umbrella","Residual","Total")
colnames(table)<-c("df","Value","pvalue")
out<-list(Z,table)
names(out)<-list("Z","partition")
return(out)
}

xtoCMH<-function(x)
{
n<-nrow(x)
m<-ncol(x)
scores<-sort(unique(c(x)))
out<-array(0,dim=c(m,length(scores),n),dimnames=list(colnames(x),scores,Consumer=1:n))
for (i in 1:n)
{
for (j in 1:m)
{
code<-(1:length(scores))[x[i,j]==scores]
out[j,code,i]<-out[j,code,i]+1
}
}
return(out)
}

CMHanalysis<-function(x,ttscore=F,tscore=F)
{
#perform CMH analysis
#Note my data has the last two dimensions swapped compared to John's n
#also my x is the transpose of Johns idat
#variables with same name as Johns mean exactly same thing (dimension order etc)
data<-xtoCMH(x)
nt<-dim(data)[1]
nj<-dim(data)[3]
ncat<-dim(data)[2]
nuw<-apply(data,c(1,2),sum)
nw<-apply(data,2,sum)
nvw<-apply(data,c(3,2),sum)
ntx<-nt^2
#Extended Stuart Test
s<-0
for (j in 1:nj)
{
s<-s-nvw[j,]%o%nvw[j,]/ntx
}
s<-s+diag(nw)/nt
m<-ncat-1
th<-(nuw-rep(1,nt)%o%nw/nt)[,-ncat] #don't use last column in th
y<-solve(s[-ncat,-ncat])
b<-0
a<-array(0,dim=c(m))
for (i in 1:nt)
{
for (j in 1:m)
{
a[j]<-0
for (k in 1:m) {a[j]<-a[j]+th[i,k]*y[k,j]}
}
for (l in 1:m)
{
b<-b+a[l]*th[i,l]
}
}
df<-nt-1
b<-b*df/nt
table<-rbind(c(b,df*m,1-pchisq(b,df*m)))
#Mean Score Differences Test
if (!ttscore) {ttscore<-1:ncat}
if (!tscore) {tscore<-1:nt}
w<-0
for (j in 1:nj)
{
v<-0
av<-0
for (i in 1:nt)
{
av<-av+ttscore[x[j,i]]
v<-v+ttscore[x[j,i]]^2
}
av<-av/nt
v<-(v/nt-av^2)*nt/(nt-1)
w<-w+v
}
gt<-0
ts<-array(0,dim=c(nt))
for (i in 1:nt)
{
ts[i]<-0
for (j in 1:nj)
{
gt<-gt+ttscore[x[j,i]]
ts[i]<-ts[i]+ttscore[x[j,i]]
}
}
ave<-0
for (i in 1:nt)
{
ave<-ave+tscore[i]
}
ave<-ave/nt
bl<-0
bll<-0
ww<-0
for (i in 1:nt)
{
ww<-ww+(tscore[i]-ave)^2
bll<-bll+(ts[i]-gt/nt)*tscore[i]
bl<-bl+(ts[i]-gt/nt)^2
}
bl<-bl/w
df<-nt-1
table<-rbind(table,c(bl,df,1-pchisq(bl,df)))
#Non-Mean Score Differences Test
q2<-b-bl
df<-(nt-1)*(m-1)
table<-rbind(table,c(q2,df,1-pchisq(q2,df)))
#Correlation Test
bll<-(bll^2)/(w*ww)
table<-rbind(table,c(bll,1,1-pchisq(bll,1)))
#
rownames(table)<-c("Extended Stuart","Mean Score Differences","Non-Mean Score Differences","Correlation")
colnames(table)<-c("Statistic","df","pvalue")
return(table)
}

ADdiscrete<-function(tab)
{
jf<-tab[,1]
qx<-tab[,2]/sum(tab[,2])
nc<-dim(tab)[1]
n<-sum(tab[,1])
ncc<-nc-1
#
hj<-0
zj<-0
a2<-0
for (i in 1:ncc)
{
zj<-zj+jf[i]-n*qx[i]
hj<-hj+qx[i]
a2<-a2+zj*zj*qx[i]/(hj*(1-hj))
}
a2<-a2/n
names(a2)<-NULL
return(a2)
}

GOF<-function(y,q=0,qcomps=(dim(y)[1]-q):(dim(y)[1]-1))
{
#input should be a table with first column the observed and second column the expected
#input columns should sum to exactly the same answer (don't round expected)
x<-y[,1]
probs<-y[,2]/sum(y[,2])
m<-length(x)
n<-sum(x)
if (!sum(y[,1])==sum(y[,2]))
{
stop("input columns should sum to exactly the same answer (don't round expected)")
}
#q is the number of params that were estimated to get the expected values
#qcomp "assigns" which components were "used up" in this process - defaults to highest order components
if (q>=m-1) stop("invalid q")
if (q>0)
{
qcomps<-qcomps[1:q]
qcomps<-(1:(m-1))[qcomps]
cat("Note ",q," parameters were estimated to calculate expected values\n",sep="")
cat("The particular component indices given zero df to compensate are:\n",sep="")
cat(qcomps,"\n")
if (!length(qcomps)==q) stop("invalid qcomps")
cat("This information is only used to calculate partition dfs and p-values\n\n")
}
#
Vr<-apply((rep(1,m)%o%x)*(orthpoly(probs,1:m))/sqrt(n),1,sum)[-1]
dfs<-rep(1,length(Vr))
if (q>0) {dfs[qcomps]<-0} #kill dfs of qcomps indices
table<-cbind(dfs,Vr^2)
if (m-1>4) {table<-rbind(table[1:4,],apply(table[5:(m-1),],2,sum))}
table<-rbind(table,apply(table,2,sum))
pvals<-array(NA,dim=length(table[,1]))
pvals[table[,1]>0]<-(1-pchisq(table[table[,1]>0,2],table[table[,1]>0,1]))
table<-cbind(table,pvals)
rownames(table)<-c("Location","Dispersion","Skewness","Kurtosis","Residual")[1:(dim(table)[1])]
rownames(table)[dim(table)[1]]<-"Total"
colnames(table)<-c("df","SS","pvalue")
#
out<-list(ADdiscrete(tab),Vr,table)
names(out)<-c("AndersonDarling","Vr","partition")
return(out)
}