######  Contributed by James Ng   ######

####Example 7.7#####

### Initial Part Sets up Data #

## Reads in the fMRI data in brainm.dat (a 128 by 234 matrix)
## and arranges them into a list with 54 components
## data$L1S1, data$L1S2, ... data$L1S6, ... data$L9S1, data$L9S2, ... data$L9S6.
## L is for 'location', S is for 'stimulus'
## Each component is a matrix with one column (a time series of length 128)
## for each subject.  Numbers of subjects varies:
##
## S1: "Awake-Brush"  5 subjects
## S2: "Awake-Heat"   4 subjects
## S3: "Awake-Shock"  5 subjects
## S4: "Low-Brush"    3 subjects
## S5: "Low-Heat"     5 subjects
## S6: "Low-Shock"    4 subjects

## Then Figure 7.1 in the text (mean responses per stimulus at location 1) is duplicated
## Similar figures are produced for the other 8 locations

## Read in the data - a 128 by 234 matrix
z.all=ts(read.table("http://www.stat.ualberta.ca/~wiens/stat479/S&Sdatasets/brainm.dat"))
# or:
z.all = read.table("brainm.dat")

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

## Arrange data in 9 groups of 26 columns; each group is a "location":
L = NULL
for (j in 1:9) {
        L = cbind(L, z.all[, j + 9*(0:25)])
}
colnames(L) = NULL
# L has all data at location 1 (first 26 columns),
# then all data at location 2 (next 26 columns), ...
# then all data at location 9 (final 26 columns)

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

## Set the indicators of "location"
where.loc = vector(length = 234)
for (l in 1:9) where.loc[(1:26) + (l-1)*26] = l
#where.loc   # L[,where.loc == l] consists of the data from location "l"


Location.names =c("Cortex 1", "Cortex 2", "Cortex 3", "Cortex 4", "Caudate",
        "Thalamus 1", "Thalamus 2", "Cerebellum 1", "Cerebellum 2")
cat("Here are the labels of the locations:","\n")
as.matrix(Location.names)

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

## Set the indicators of "stimuli"
stim = c(rep(1,5), rep(2,4), rep(3,5), rep(4,3), rep(5,5), rep(6,4))
where.stim = rep(stim,9)
#where.stim   # L[,where.stim == s] consists of the data from stimulus "s"


Stimulus.names = c("Awake-Brush", "Awake-Heat", "Awake-Shock",
       "Low-Brush", "Low-Heat", "Low-Shock")
cat("Here are the labels of the stimuli:","\n")
as.matrix(Stimulus.names)

##############
## Create a list with components L1S1, L1S2, ... L1S6, L2S1, ..., L9S6
## Each component is a matrix with one column (a series of length 128) per subject

data = vector("list", 54)
for (l in 1:9) {for(s in 1:6) {
        data[[(l-1)*6+s]] = L[, (where.loc == l & where.stim == s)]
        }}

datanames = vector(length=54)
for(l in 1:9) datanames[(l-1)*6+1:6]=paste("L",l,"S",1:6,sep="")
names(data) = datanames
#data

# A check:
#L5S3test = L[, (where.loc == 5 & where.stim == 3)]
#L5S3test - data$L5S3   # Should be all zeros

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

## Figure 7.1 in text; average reponses to stimuli at location 1:
#par(mfcol = c(3,2))
#for(s in 1:6) plot(ts(rowMeans(data[[s]])),
#     ylim = c(-1,1), ylab = Stimulus.names[s])
#title(main = paste("\nMean responses at", Location.names[1], "to various stimuli"), outer = T)

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

## Similar plots at other locations:
#for(l in 2:9) {
#win.graph()
#par(mfcol = c(3,2))
#for(s in 1:6) plot(ts(rowMeans(data[[(l-1)*6+s]])),
#     ylim = c(-1,1), ylab = Stimulus.names[s])
#title(main = paste("\nMean responses at", Location.names[l], "to various stimuli"), outer = T)
#}


################ FIGURE 7.9 #######################

ndata=128;
n=128;
nvar=9;

n2=n/2+1;

#number of vector series for each cell
N = c(5,4,5,3,5,4)
ns=sum(N);
nc=1;
nq=6;


#Bandwidth
L=3;
m=(L-1)/2;

#Define Design Matrix
# For the 1-way anova (Figure 7.9); 6 treatment combinations S1,...S6
z1 = outer(X=rep(1,N[1]), Y = c(1,0,0,0,0,0)) # Mean = mu
z2 = outer(X=rep(1,N[2]), Y = c(1,1,0,0,0,0)) # Mean = mu + alpha1
z3 = outer(X=rep(1,N[3]), Y = c(1,0,1,0,0,0)) # Mean = mu + alpha2
z4 = outer(X=rep(1,N[4]), Y = c(1,0,0,1,0,0)) # Mean = mu + alpha3
z5 = outer(X=rep(1,N[5]), Y = c(1,0,0,0,1,0)) # Mean = mu + alpha4
z6 = outer(X=rep(1,N[6]), Y = c(1,0,0,0,0,1)) # Mean = mu + alpha5

z = rbind(z1, z2, z3, z4, z5, z6) #Design matrix (26 by 6)
s = t(z)%*%z;

####Preparing to plot###
par(mfrow=c(3,3))

########Calculations########

for(q in (1:9)) {

    p=6*q

    Y = cbind(data[[p-5]],data[[p-4]],data[[p-3]],data[[p-2]],data[[p-1]],data[[p]]) #Define Group

    # Y is 128 by 26
    # Yft = fft(Y)/sqrt(n)  # Also 128 by 26; Yft[k,] = the 26 dft's at frequency omega_k = Y(omega_k)

    #Taper

    r=seq(1,ndata);
    wt=0.5*(1+cos((2*pi*(t(r)-ndata/2)/(ndata+1))));
    w=t(wt);

    Ytap = array(NA,dim=dim(Y))
        for (j in (1:ns)) {
            Ytap[,j]=w*Y[,j]; #row by row multiplication
    }


    ###Compute DFT on data###

    Yft=matrix(NA,nrow=ndata,ncol=ns);
    for (i in 1:ns) {
            Yft[,i] = fft(Ytap[,i])/sqrt(ndata);      
	}

    Yf = t(Yft); #Yf is 26 by 128


    ###Calculation of Error Spectrum###

    sy = array(NA,dim=n)
    ssreg = array(NA,dim=n)
    ssef = array(NA,dim=n)
    sser = array(NA,dim=n)

    ser = array(NA,dim=(n/2+1))
    sef = array(NA,dim=(n/2+1))

    #Full Model Error Power

    for (k in (1:n)) {

        sy[k]=Re(Conj(t(Yf[,k]))%*%Yf[,k])
          ssreg[k]= Re(Conj(t(Yf[,k]))%*%z%*%solve(s)%*%t(z)%*%Yf[,k])
          ssef[k]=sy[k]-ssreg[k];

    #Reduced Model Error Power

        ssreg[k]= Re(Conj(t(Yf[,k]))%*%z[,1]%*%solve(s[1,1])%*%t(z[,1])%*%Yf[,k])
        sser[k]=sy[k]-ssreg[k];
    }


 

    ##Smooth Error Power Components##

	ser = filter(sser, rep(1/L, L), circular = T)
	sef = filter(ssef, rep(1/L, L), circular = T)


    ##F-Statistic##
    nq1 = 1;
    num_df = 2*L*(nq-nq1);
    den_df = 2*L*(ns-nq)
    con=den_df/num_df;

    F = array(NA,dim=(n/2+1)) #F Statistics

    for (k in (1:(n/2+1))) {
            F[k]=con*(ser[k]-sef[k])/sef[k];
    }


    ##Plotting F-stats with critical frequency##

    Fr = array(NA,dim=(n/2+1)) #frequencies
    for (k in (1:(n/2+1))) {
            Fr[k]=0.5*(k-1)/n;
    }

    plot(Fr,F,xlab="cycles/sec",ylab="F-Statistic",main=Location.names[q],type="l",xlim=c(0,0.25),ylim=c(0,8))
    abline(h=qf(.999, num_df, den_df),lty=2,col="red")
    }


################# END OF EXAMPLE 7.7 #####################


################# EXAMPLE 7.8 ######################

#number of vector series for each cell
N1 = c(5,4,5,3,5,4)
ns=sum(N);
nc=1;
nq=6;


#Define Design Matrix
# For the 3-way anova:
z1 = outer(X=rep(1,N[1]), Y = c(1,0,0,0,0,0)) 
z2 = outer(X=rep(1,N[2]), Y = c(1,0,1,0,0,0)) 
z3 = outer(X=rep(1,N[3]), Y = c(1,0,0,1,0,0)) 
z4 = outer(X=rep(1,N[4]), Y = c(1,1,0,0,0,0)) 
z5 = outer(X=rep(1,N[5]), Y = c(1,1,1,0,1,0)) 
z6 = outer(X=rep(1,N[6]), Y = c(1,1,0,1,0,1)) 

z = rbind(z1, z2, z3, z4, z5, z6) #Design matrix (26 by 6)
s = t(z)%*%z;

####Preparing to plot###

win.graph()
par(mfrow=c(5,3))

########Calculations########

for(q in (1:5)) {

    if (q<5) {
        p=6*q
        Y = cbind(data[[p-5]],data[[p-4]],data[[p-3]],data[[p-2]],data[[p-1]],data[[p]]) #Define Group
    }else{
        r=49
        Y = cbind(data[[r]],data[[r+1]],data[[r+2]],data[[r+3]],data[[r+4]],data[[r+5]]) #Define Cerebellum2
    }
        # Y is 128 by 26
        #Yft = fft(Y)/sqrt(n)  # Also 128 by 26; Yft[k,] = the 26 dft's at frequency omega_k = Y(omega_k)

    ##Taper

    r=seq(1,ndata);
    wt=0.5*(1+cos((2*pi*(t(r)-ndata/2)/(ndata+1))))
    w=t(wt);

    Ytap = array(NA,dim=dim(Y))
    for (j in (1:ns)) {
            Ytap[,j]=w*Y[,j]  } #row by row multiplication


    ###Compute DFT on data ###

    Yft=matrix(NA,nrow=ndata,ncol=ns)
    for (i in 1:ns) {
            Yft[,i] = fft(Ytap[,i])/sqrt(ndata) }

    Yf = t(Yft); #Yf is 26 by 128


    #####Calculation of Error Spectrum######



    #set up inner loop

    

        sy = array(NA,dim=n)
        ssreg = array(NA,dim=n)
        ssef = array(NA,dim=n)
        sser = array(NA,dim=n)
        ser = array(NA,dim=(n/2+1))
        sef = array(NA,dim=(n/2+1))

        ssreg.inter = array(NA,dim=n)
        ssef.inter = array(NA,dim=n)
        sser.inter = array(NA,dim=n)
        ser.inter = array(NA,dim=(n/2+1))
        sef.inter = array(NA,dim=(n/2+1))
	  nq2.inter = 2

        ssreg.stim = array(NA,dim=n)
        ssef.stim = array(NA,dim=n)
        sser.stim = array(NA,dim=n)
        ser.stim = array(NA,dim=(n/2+1))
        sef.stim = array(NA,dim=(n/2+1))
	  nq2.stim = 2

        ssreg.con = array(NA,dim=n)
        ssef.con  = array(NA,dim=n)
        sser.con  = array(NA,dim=n)
        ser.con  = array(NA,dim=(n/2+1))
        sef.con  = array(NA,dim=(n/2+1))
	  nq2.con = 1


        for (k in (1:n)) {

            #Full Model Power

            sy[k]=Re(Conj(t(Yf[,k]))%*%Yf[,k])
            ssreg[k]= Re(Conj(t(Yf[,k]))%*%z%*%solve(s)%*%t(z)%*%Yf[,k])
            ssef[k]=sy[k]-ssreg[k]

            #Reduced Models Error Power

		#Interaction
		ssreg.inter[k]= Re(Conj(t(Yf[,k]))%*%z[,-c(5:6)]%*%solve(s[-c(5:6),-c(5:6)])%*%t(z[,-c(5:6)])%*%Yf[,k])
 

		# Stimulus
            ssreg.stim[k]= Re(Conj(t(Yf[,k]))%*%z[,-c(3:4)]%*%solve(s[-c(3:4),-c(3:4)])%*%t(z[,-c(3:4)])%*%Yf[,k])

            # Consciousness
            ssreg.con[k]= Re(Conj(t(Yf[,k]))%*%z[,-2]%*%solve(s[-2,-2])%*%t(z[,-2])%*%Yf[,k])

            sser.inter[k]=sy[k]-ssreg.inter[k]
		sser.stim[k]=sy[k]-ssreg.stim[k]
		sser.con[k]=sy[k]-ssreg.con[k]
            }

        ##Smooth Error Power Components##


	sef = filter(ssef, rep(1/L, L), circular = T)
	ser.inter = filter(sser.inter, rep(1/L, L), circular = T)
	sef.inter = sef
	ser.stim = filter(sser.stim, rep(1/L, L), circular = T)
	sef.stim = sef
	ser.con = filter(sser.con, rep(1/L, L), circular = T)
	sef.con = sef

 

        ##F-Statistic##
        
        num_df.inter = 2*L*(nq2.inter);
        num_df.con = 2*L*(nq2.con);
        num_df.stim = 2*L*(nq2.stim);


        den_df = 2*L*(ns-nq)
        con.inter=den_df/num_df.inter;
        con.con=den_df/num_df.con;
        con.stim=den_df/num_df.stim;

        F.inter = array(NA,dim=(n/2+1)) #F Statistics

        for (k in (1:(n/2+1))) {
                F.inter[k]=con.inter*(ser.inter[k]-sef.inter[k])/sef.inter[k] }


	F.stim = array(NA,dim=(n/2+1)) #F Statistics

        for (k in (1:(n/2+1))) {
                F.stim[k]=con.stim*(ser.stim[k]-sef.stim[k])/sef.stim[k] }


	F.con = array(NA,dim=(n/2+1)) #F Statistics

        for (k in (1:(n/2+1))) {
                F.con[k]=con.con*(ser.con[k]-sef.con[k])/sef.con[k] }




        ##Plotting F-stats with critical frequency##

        Fr = array(NA,dim=(n/2+1)) #frequencies
        for (k in (1:(n/2+1))) {
                Fr[k]=0.5*(k-1)/n }
	 main.name = Location.names[q]
	 if(q==5) main.name = "Cerebellum"
       plot(Fr,F.stim,xlab="cycles/sec",ylab="F-Stimulus",main=main.name,type="l",xlim=c(0,0.25),ylim=c(0,8))
        	abline(h=qf(.999, num_df.stim, den_df),lty=2,col="red")
       plot(Fr,F.con,xlab="cycles/sec",ylab="F-Consciousness",main=main.name,type="l",xlim=c(0,0.25),ylim=c(0,8))
           abline(h=qf(.999, num_df.con, den_df),lty=2,col="red")
       plot(Fr,F.inter,xlab="cycles/sec",ylab="F-Interactions",main=main.name,type="l",xlim=c(0,0.25),ylim=c(0,8))
           abline(h=qf(.999, num_df.inter, den_df),lty=2,col="red")

 
    }