######## R program by: Andrada E. Ivanescu

######## R program for Sec 4.3. Functional predictors sampled with moderate sparsity


library(SemiPar)
library(mgcv)
library(refund)


seed.start = 500

B=1						#number of simulated datasets


by = 0.1				# density grid on which functions are observed 
t=seq(0.1,6, by=by)		# grid on which response curves are observed
s=seq(0.1,5, by=by)		# grid on which X1(s) predictor curves are observed
r=seq(0.1,7, by=by)		# grid on which X2(r) predictor curves are observed
n=200					# number of curves
VarE=1
VarX1=0.2^2
VarX2=0.2^2



# keep track of the seed for reproducibility
seeds = matrix(seed.start + 0:(B-1), nrow=B, ncol=1)

###----------------------intercept ---------------------------------###
trueBeta0<-matrix(2*exp(-1*(t-2.5)^2),ncol=1); trueBeta0<-matrix(trueBeta0,nrow=1);			

trueBeta0_mat<-matrix(trueBeta0,nrow=n,ncol=length(t),byrow=T)


###---------------------bivariate beta param -----------------------###
trueBeta1<-matrix(0,nrow=length(s),ncol=length(t))
trueBeta2<-matrix(0,nrow=length(r),ncol=length(t))
for(s1 in 1:length(s)){
for(t1 in 1:length(t)){
	trueBeta1[s1,t1]<-sin(2*s[s1]*pi/10)*cos(2*t[t1]*pi/6)
	}}

for(t1 in 1:length(t)){
for(r1 in 1:length(r)){
	trueBeta2[r1,t1]<-sqrt(r[r1]*t[t1])/4.2
	}}
	

#############################################################
############## Initialize matrices to #######################
#######  store estimates for each simulated dataset  ########
#############################################################

SNR<-matrix(nrow=B,ncol=1)

timesX1<-matrix(nrow=B,ncol=1)
timesX2<-matrix(nrow=B,ncol=1)


m100_hatBeta0<-matrix(nrow=B, ncol=length(t))
m100_hatBeta1_sparse<-array(dim=c(B,nrow=length(s) , ncol=length(t)))
m100_hatBeta2_sparse<-array(dim=c(B,nrow=length(r) , ncol=length(t)))

betaHat0.mult<-matrix(nrow=B, ncol=length(t))
betaHat1.mult<-array(dim=c(B,nrow=length(s) , ncol=length(t)))
betaHat2.mult<-array(dim=c(B,nrow=length(r) , ncol=length(t)))


ISE_Beta0<-matrix(nrow=B, ncol=1)
ISE_Beta1<-matrix(nrow=B, ncol=1)
ISE_Beta2<-matrix(nrow=B, ncol=1)


ISE_Beta0_pfr<-matrix(nrow=B, ncol=1)
ISE_Beta1_pfr<-matrix(nrow=B, ncol=1)
ISE_Beta2_pfr<-matrix(nrow=B, ncol=1)


ICBeta0<-matrix(nrow=B, ncol=1)
ICBeta1<-matrix(nrow=B, ncol=1)
ICBeta2<-matrix(nrow=B, ncol=1)



ICBeta0_pfr<-matrix(nrow=B, ncol=1)
ICBeta1_pfr<-matrix(nrow=B, ncol=1)
ICBeta2_pfr<-matrix(nrow=B, ncol=1)


covBeta0<-matrix(0,nrow=B, ncol=length(t))
covBeta1<-array(0,dim=c(B,nrow=length(s) , ncol=length(t)))
covBeta2<-array(0,dim=c(B,nrow=length(r) , ncol=length(t)))


IWBeta0<-matrix(nrow=B, ncol=1)
IWBeta1<-matrix(nrow=B, ncol=1)
IWBeta2<-matrix(nrow=B, ncol=1)



IWBeta0_pfr<-matrix(nrow=B, ncol=1)
IWBeta1_pfr<-matrix(nrow=B, ncol=1)
IWBeta2_pfr<-matrix(nrow=B, ncol=1)

IR2<-matrix(nrow=B, ncol=1)
IR2_pfr<-matrix(nrow=B, ncol=1)


######################################################################
################# pfr function (Goldsmith et al. 2011)  ##############

pfr <-
function(Y, covariates=NULL, funcs, kz=30, kb=30, smooth.cov=FALSE, family="gaussian", ...) {
	require(mgcv)

	kb = min(kz, kb)
	n = length(Y)
	t = seq(0, 1, length = dim(funcs)[2])
	p = ifelse(is.null(covariates), 0, dim(covariates)[2])
	
	## get the eigen decomposition of the smoothed variance matrix
	if(smooth.cov){
		G2 <- var(funcs)
		M <- length(t)
		diag(G2)= rep(NA, M)
		g2 <- as.vector(G2)
		## define a N*N knots for bivariate smoothing
		N <- 10

		## bivariate smoothing using the gamm function
		t1 <- rep(t, each=M)
		t2 <- rep(t, M)
		newdata <- data.frame(t1 = t1, t2 = t2)
		K.0  <- matrix(predict(gamm(as.vector(g2)  ~ te(t1, t2, k = N))$gam, newdata), M, M) # smooth K.0
		K.0 <- (K.0 + t(K.0)) / 2    
		
		eigenDecomp <- eigen(K.0)
	} else {
		varFuncs <- var(funcs)
		eigenDecomp <- eigen(varFuncs)
	}
	
	psi = eigenDecomp$vectors[,1:kz]

	# set the basis to be used for beta(t)
	num=kb-2
	qtiles <- seq(0, 1, length = num + 2)[-c(1, num + 2)]
	knots <- quantile(t, qtiles)
	phi = cbind(1, t, sapply(knots, function(k) ((t - k > 0) * (t - k))))

	## set up the design and penalty matrix
	C = funcs %*% eigenDecomp$vectors[ ,1:kz ]
	J = t(eigenDecomp$vectors[,1:kz]) %*% phi
	CJ = C %*% J

	X = cbind(1, covariates, CJ)
	D = diag(c(rep(0, 1+p+2), rep(1, kb-2)))
	
	## fit the model
	gam = gam(Y~X-1, paraPen=list(X=list(D)), family=family, method="REML", ...)

	## get the coefficient and betaHat estimates
	coefs = gam$coef
	fitted.vals <- as.matrix(X[,1:length(coefs)]) %*% coefs
	
	beta.covariates = coefs[1:(p+1)]
	betaHat <- phi %*% coefs[-(1:(p+1))]

	## get the covariance matrix of the estimated functional coefficient
	varBeta=gam$Vp[-(1:(1+p)),-(1:(1+p))]
	varBetaHat=phi%*%varBeta%*%t(phi)
	
	## construct upper and lower bounds for betahat
	Bounds = cbind(betaHat + 1.96*(sqrt(diag(varBetaHat))), 
		betaHat - 1.96*(sqrt(diag(varBetaHat))))
	
	ret <- list(gam, fitted.vals, betaHat, beta.covariates, X, phi, psi, varBetaHat, Bounds)
	names(ret) <- c("gam", "fitted.vals", "betaHat", "beta.covariates", "X", "phi", 
		"psi", "varBetaHat", "Bounds")
	ret
}



######################################################################
####------------------- Start B simulations ----------------------####
######################################################################


for(b in 1:B){


		set.seed(seeds[b])


###---------------------X1 funct pred -----------------------###	
###----------------------- mean 0 ---------------------------###
###----------------- adapted from Goldsmith -----------------###

	funcs1.true <- matrix(0, nrow=n, ncol=length(s))
	for(i2 in 1:n){
 		for(j2 in 1:2){
			e=rnorm(2, 0, 1/j2^(2))
			funcs1.true[i2,]=funcs1.true[i2,]+e[1]*cos((2*pi/10)*s*j2) 
			funcs1.true[i2,]=funcs1.true[i2,]+e[2]*sin((2*pi/10)*s*j2) 
		}
	}	
	

	X1=funcs1.true
	

###---------------------X2 funct pred -----------------------###
###----------------------- mean 0 ---------------------------###
###-------- adapted from Brownian Bridge (BB)----------------###
###------- for BB, see Shorack (2000), p. 304 ---------------###




K<-4; M<-length(r); N<-n

Un<-matrix(0,nrow=K,ncol=N)
for(n in 1:N){Un[,n]<-rnorm(K,mean=0,sd=1)}

U<-matrix(0,nrow=K,ncol=N)
for(k in 1:K){U[k,]<-(Un[k,])/(k*pi)}

SIB<-matrix(0,nrow=K,ncol=M)
for(j in 1:K){ SIB[j,]<-2*(sqrt(2))*sin(1*pi*j*r/7) }		

BB<-t(SIB)%*%U		
	
funcs2.true=t(BB)

	
X2=funcs2.true


####################### generate outcomes #########################

L1 <- matrix(by, ncol=length(s), nrow=n)
L2 <- matrix(by, ncol=length(r), nrow=n)
LX1 <- L1*X1
LX2 <- L2*X2

eta <- LX1 %*% trueBeta1 + LX2 %*% trueBeta2
#################################################################################



errors=matrix(scale(rnorm(n*length(t),mean=0,sd=sqrt(VarE))), nr=n)

outcomes <- trueBeta0_mat+ eta + errors

SNR[b]<-sd(as.vector(trueBeta0_mat+eta))/sd(as.vector(errors))



################## X1 and X2 are observed with noise: V1 and V2 ################

V1=funcs1.true+matrix(rnorm(length(X1), mean=0, sd=sqrt(VarX1)),nrow=n,ncol=length(s))

V2=funcs2.true+matrix(rnorm(length(X2), mean=0, sd=sqrt(VarX2)),nrow=n,ncol=length(r))


#########  R code smoothing approach used in Goldsmith et al. (2011) 
#########  to reconstruct the trajectories for each functional predictor

########## SPARSE

NGRIDPTS=12							# number of sampled points for each function
N.fit=length(s) 					# number of points in the grid for the estimated functions


		## sparse data 

		N=NGRIDPTS		
		M=n
		J=1
		tstart=0.1
		tlength=5
		N.smooth.single=20
		N.smooth=10


		unique.average <- function(index, y) {
    
    		n <- sum(!is.na(y))
		    y.new <- rep(NA, n)
	    	index.new <- rep(NA, n)
    		count <- rep(0, n)
    
     
		    y.new[1] <- y[1]
    		index.new[1] <- index[1] 
		    count[1] <- 1
	    	current <- 1
    
		    if(n > 1) {
    		    for(i in 2:n) {
        		    if(index[i]==index[i-1]) {
            	    	count[i] <- count[i-1] + 1
                		y.new[current] <- ( y.new[current] * count[i-1] + y[i] )/count[i]
	        	        count[i-1] <- count[i]
                   
        	    }
    	        else {
		                current <- current + 1
    	            	y.new[current] <- y[i] 
        	    	    count[i] <- 1
            		    index.new[current] <- index[i] 
	    	        }
    		    }
	    	}
    	
	    	return(list(index=index.new, y=y.new))
		}


		sampled_s=t(sapply(1:n, function(u) sort(sample(1:length(s), size=NGRIDPTS  ))))
		tx=t(sapply(1:n, function(u) s[sampled_s[u,]]))
		y=t(sapply(1:n, function(u) V1[u,sampled_s[u,]]))

		



		t.knot <- seq(tstart, tstart + tlength-by, length=N.smooth.single)
    	t.fit<-s
    	
    	
		N.row <- dim(y)[1]
    
		y1 <- matrix(0, M, N)
		index <- matrix(0, M, N)
		t.fit.matrix <- matrix(rep(t.fit, each=N), nrow=N)


		for(m in 1:M ) {
			row.ind <- m
			index.temp <- apply( (t.fit.matrix - tx[row.ind, ])^2, 1, which.min )
			fit <- unique.average(index.temp, y[row.ind,] )
			y1[row.ind,] <- fit$y
			index[row.ind, ] <- fit$index
		}
    
		n.grid <- apply(!is.na(y1), 1, sum)


		library(SemiPar)
		mu.fit <- rep(0,  N.fit)
    
		mu.obs <- matrix(0, M, N)
    
		###  divide the interval into "N.fit" bins to reduce computational burden for smoothing 
		mu.sum <- rep(0, N.fit)
		mu.count <- rep(0, N.fit) 
		mu.mean <- rep(0, N.fit)
		for (m in 1:M) {
			row.ind <- m
			row.len <- n.grid[row.ind]
			mu.count[ index[row.ind, 1:row.len] ] <- mu.count[ index[row.ind, 1:row.len] ] + 1
			mu.sum[ index[row.ind, 1:row.len] ] <- mu.sum[ index[row.ind, 1:row.len] ] + y1[row.ind, 1:row.len ]
		} 
		mu.mean <- ifelse(mu.count==0, NA,  mu.sum/mu.count)    

   
    
		temp.data <- data.frame(y2 = mu.mean,  x2 = t.fit )
		attach(temp.data, pos=3)
		fit.mu <- spm( y2~f(x2,knots=t.knot) ,omit.missing=TRUE)
	
		newx <- data.frame(x2=t.fit)
		mu.fit <- predict(fit.mu,newx)
	
		 
		detach(pos=3)
		mu.obs <- matrix( mu.fit[ index ], ncol=N)

		resd <- matrix(0, nrow=M, ncol=N) 
		resd <- y - mu.obs

		###  divide the interval into "N.fit" bins to reduce computational burden for smoothing 
		cov.sum <- matrix(0, N.fit, N.fit)
		cov.count <- matrix(0, N.fit, N.fit) 
		cov.mean <-  matrix(0, N.fit, N.fit)
		for (m in 1:M) {
			row.ind <- m
			len.ind <- n.grid[row.ind]
			temp.ind <- index[row.ind, 1:len.ind]
			cov.count[ temp.ind, temp.ind ] <- cov.count[ temp.ind, temp.ind ] + 1
			cov.sum[ temp.ind, temp.ind ] <- cov.sum[ temp.ind, temp.ind ] + resd[row.ind, 1:len.ind ] %*% t( resd[row.ind, 1:len.ind])
		} 
		cov.mean <- ifelse(cov.count==0, NA,  cov.sum/cov.count)    

		cov.mean.nodiag <- cov.mean
		diag(cov.mean.nodiag) <- rep(NA, N.fit)
		resd.pair <- cbind( as.vector(cov.mean.nodiag), rep(t.fit, each=N.fit), rep(t.fit, N.fit) )
		resd.pair.diag <- cbind( diag(cov.mean), t.fit) 


		G <- matrix(0, N.fit, N.fit)

		### load the "SemiPar" package which contains functions for semiparametric smoothing 
		### using linear mixed models
		t.smooth <- seq(tstart, tstart + tlength-by ,length=N.smooth)
        
		myknots.t <- data.frame(x1t=rep(t.smooth, each=N.smooth), x2t=rep(t.smooth, N.smooth))         
		data.gt <- data.frame(gt=resd.pair[,1], x1t=resd.pair[,2], x2t=resd.pair[,3])
		attach(data.gt, pos=4)
		fit.t <- spm(gt ~ f(x1t, x2t, knots=myknots.t),omit.missing=TRUE)
		newdata.t <- data.frame(x1t=rep(t.fit, each=N.fit),x2t=rep(t.fit, N.fit))
		pred.t <- predict(fit.t, newdata.t)
		detach(pos=4)

		s.g <- matrix(pred.t, N.fit)
		G <- ( s.g + t(s.g) )/2

		eigenDecomp <- eigen(G)

		par(mfrow=c(2,3))
		for(i in 1:6) {
			plot(eigenDecomp$vectors[,i], type="l")
		}


		temp.data <- data.frame(var.y = as.vector(resd^2),  var.x = as.vector(tx) )
		attach(temp.data, pos=6)
		fit.var <- spm( var.y ~ f(var.x,knots=t.knot) ,omit.missing=TRUE)
		newx <- data.frame(var.x=t.fit)
		var.fit <- predict(fit.var,newx)          
		detach(pos=6)

		var.noise <- mean( var.fit - diag(matrix(pred.t,N.fit)) ) 

		var.noise
		
		
		###############################################
		
		if(var.noise<0){
				timesX1[b]<-var.noise
				var.noise=.002
		}
		
    	###############################################
	
    
	    fpca1.value <- eigenDecomp$values* tlength / N.fit

    	fpca1.value <- ifelse(fpca1.value>=0, fpca1.value, 0)

    
		K1=min(35, N)
		
		fpca1.vectors <- eigenDecomp$vectors[, 1:K1]*sqrt(N.fit/tlength)


		###     The eigenfunctions are unique only up to a change of signs.
		###     Select the signs of eigenfunctions so that the integration over the domain 
		###     is non-negative
    	for(i2 in 1:K1) {
	        v2 <- fpca1.vectors[,i2]
        	tempsign <- sum(v2)
    	    fpca1.vectors[,i2] <- ifelse(tempsign<0, -1,1) * v2
	    }


	
		lam1=fpca1.value[1:K1]
		phi1=fpca1.vectors



		###
		###     Estimatimate principal component scores and predict various functions
		###

		library(MASS) 

		index <- matrix(0, M, N)
		mu.hat <- matrix(0, M, N)    
		y.center <- matrix(0, M, N)

		phi1.hat <- array(0, dim=c( M, N, K1) )

		t.fit.matrix <- matrix(rep(t.fit, each=N), nrow=N)

		for(m in 1:M ) {
    	    row.ind <- m
        	index <- apply( (t.fit.matrix - tx[row.ind, ])^2, 1, which.min )
    	    y.center[row.ind, ] <- y[ row.ind ,] - mu.fit[ index ]
	        phi1.hat[row.ind, , ] <- phi1[index, ]
		}



		est.si1 <- matrix(0, nrow=M, ncol=K1)

		est.si1.sd <- matrix(0, nrow=M, ncol=K1)


		y.hat.center <- matrix(0, M, N.fit) 
		y.hat <- matrix(0, M, N.fit)
		y.hat.subject <- matrix(0, M, N.fit)

		y.hat.sd <- matrix(0, M, N.fit)
		y.hat.subject.sd <- matrix(0, M, N.fit)



		for(m in 1:M) {
			cov.y <- matrix(0, N, N)
			cov.xi.y <- matrix(0, K1,  N )

			ind1 <- m
			cov.y [ 1:N , 1:N ] <- phi1.hat[ind1,,] %*% diag( lam1) %*%  t( phi1.hat[ind1,,] )  + diag( rep(var.noise, N) )

			ind1 <- m
			cov.xi.y[1:K1, 1:N ] <- diag( lam1) %*% t( phi1.hat[ind1,,] )


			temp.mat <- cov.xi.y %*% ginv(cov.y)
			score <- temp.mat %*% as.vector( t(y.center[m,]) ) 
			var.score <- diag( c(lam1 ) ) -  temp.mat %*% t(cov.xi.y)

			est.si1[m ,] <- score[1:K1]

			y.hat.subject[m, ] <- mu.fit + phi1 %*% est.si1[m,] 
			y.hat.subject.sd[m,] <-  sqrt( diag( phi1 %*% var.score[1:K1, 1:K1] %*% t( phi1 ) ) )

			row.ind <- m
			y.hat.center[ row.ind ,] <-  phi1 %*% est.si1[m,]
			y.hat[row.ind,] <- y.hat.center[ row.ind ,]  +  mu.fit

			temp1 <- cbind( phi1)
			temp2 <- var.score[ c(1:K1) , c(1:K1) ]
			y.hat.sd[row.ind, ] <- sqrt( diag( temp1 %*% temp2 %*% t(temp1) ) )   
		}	
			
			SX1<-y.hat

		###############################################
		# end sparse data estimation procedure
		###############################################



################## sparse data plot


par(mfrow=c(1,2))
## X1
range_ylim =  c(-2.5,1.5)
power10 = floor(range_ylim[1])
axisx_range = seq(s[1]-0.1, s[length(s)], by = 1) 
axisx_label =axisx_range
axisy_range =  seq(from=range_ylim[1], to=range_ylim[2], length.out=4)
axisy_label = round(axisy_range,1)


for(j in 2:3){
	
plot(s,SX1[j,],type="l",col="red",ylim=c(-2.5,1.5),ylab="",main="",lwd=2,axes=F,cex.lab=1.5);
   axis(1,at=axisx_range ,col="grey35",labels=axisx_label,cex=1.2, font=2.5,cex.axis=1.3) 
   axis(2,at=axisy_range ,col="grey35",labels=axisy_label,cex=1.2, font=2.5,cex.axis=1.3)
lines(tx[j,],y[j,],type="p",col="red",pch=21, cex=1);
lines(s,X1[j,],type="l",col="black",lty=3,lwd=2)
}




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


NGRIDPTS=25							# number of sampled points for each function
N.fit=length(r) 					# number of points in the grid for the estimated functions




		## sparse data 

		N=NGRIDPTS		
		M=n
		J=1
		tstart=0.1
		tlength=7
		N.smooth.single=20
		N.smooth=20


		unique.average <- function(index, y) {
    
    		n <- sum(!is.na(y))
		    y.new <- rep(NA, n)
	    	index.new <- rep(NA, n)
    		count <- rep(0, n)
    
     
		    y.new[1] <- y[1]
    		index.new[1] <- index[1] 
		    count[1] <- 1
	    	current <- 1
    
		    if(n > 1) {
    		    for(i in 2:n) {
        		    if(index[i]==index[i-1]) {
            	    	count[i] <- count[i-1] + 1
                		y.new[current] <- ( y.new[current] * count[i-1] + y[i] )/count[i]
	        	        count[i-1] <- count[i]
                   
        	    }
    	        else {
		                current <- current + 1
    	            	y.new[current] <- y[i] 
        	    	    count[i] <- 1
            		    index.new[current] <- index[i] 
	    	        }
    		    }
	    	}
    	
	    	return(list(index=index.new, y=y.new))
		}


		sampled_r=t(sapply(1:n, function(u) sort(sample(1:length(r), size=NGRIDPTS  ))))
		tx=t(sapply(1:n, function(u) r[sampled_r[u,]]))
		y=t(sapply(1:n, function(u) V2[u,sampled_r[u,]]))

		

		t.knot <- seq(tstart, tstart + tlength-by, length=N.smooth.single)
		
    	t.fit<-r
    	
		N.row <- dim(y)[1]
    
		y1 <- matrix(0, M, N)
		index <- matrix(0, M, N)
		t.fit.matrix <- matrix(rep(t.fit, each=N), nrow=N)


		for(m in 1:M ) {
			row.ind <- m
			index.temp <- apply( (t.fit.matrix - tx[row.ind, ])^2, 1, which.min )
			fit <- unique.average(index.temp, y[row.ind,] )
			y1[row.ind,] <- fit$y
			index[row.ind, ] <- fit$index
		}
    
		n.grid <- apply(!is.na(y1), 1, sum)


		library(SemiPar)
		mu.fit <- rep(0,  N.fit)
    
		mu.obs <- matrix(0, M, N)
    
		###  divide the interval into "N.fit" bins to reduce computational burden for smoothing 
		mu.sum <- rep(0, N.fit)
		mu.count <- rep(0, N.fit) 
		mu.mean <- rep(0, N.fit)
		for (m in 1:M) {
			row.ind <- m
			row.len <- n.grid[row.ind]
			mu.count[ index[row.ind, 1:row.len] ] <- mu.count[ index[row.ind, 1:row.len] ] + 1
			mu.sum[ index[row.ind, 1:row.len] ] <- mu.sum[ index[row.ind, 1:row.len] ] + y1[row.ind, 1:row.len ]
		} 
		mu.mean <- ifelse(mu.count==0, NA,  mu.sum/mu.count)    

    
		temp.data <- data.frame(y2 = mu.mean,  x2 = t.fit )
		attach(temp.data, pos=3)
		fit.mu <- spm( y2~f(x2,knots=t.knot) ,omit.missing=TRUE)
		newx <- data.frame(x2=t.fit)
		mu.fit <- predict(fit.mu,newx)
		detach(pos=3)
		mu.obs <- matrix( mu.fit[ index ], ncol=N)

		resd <- matrix(0, nrow=M, ncol=N) 
		resd <- y - mu.obs

		###  divide the interval into "N.fit" bins to reduce computational burden for smoothing 
		cov.sum <- matrix(0, N.fit, N.fit)
		cov.count <- matrix(0, N.fit, N.fit) 
		cov.mean <-  matrix(0, N.fit, N.fit)
		for (m in 1:M) {
			row.ind <- m
			len.ind <- n.grid[row.ind]
			temp.ind <- index[row.ind, 1:len.ind]
			cov.count[ temp.ind, temp.ind ] <- cov.count[ temp.ind, temp.ind ] + 1
			cov.sum[ temp.ind, temp.ind ] <- cov.sum[ temp.ind, temp.ind ] + resd[row.ind, 1:len.ind ] %*% t( resd[row.ind, 1:len.ind])
		} 
		cov.mean <- ifelse(cov.count==0, NA,  cov.sum/cov.count)    

		cov.mean.nodiag <- cov.mean
		diag(cov.mean.nodiag) <- rep(NA, N.fit)
		resd.pair <- cbind( as.vector(cov.mean.nodiag), rep(t.fit, each=N.fit), rep(t.fit, N.fit) )
		resd.pair.diag <- cbind( diag(cov.mean), t.fit) 
		
		
		
		
		G <- matrix(0, N.fit, N.fit)

		### load the "SemiPar" package which contains functions for semiparametric smoothing 
		### using linear mixed models
		t.smooth <- seq(tstart, tstart + tlength-by ,length=N.smooth)
        
		myknots.t <- data.frame(x1t=rep(t.smooth, each=N.smooth), x2t=rep(t.smooth, N.smooth))         
		data.gt <- data.frame(gt=resd.pair[,1], x1t=resd.pair[,2], x2t=resd.pair[,3])
		attach(data.gt, pos=4)
		fit.t <- spm(gt ~ f(x1t, x2t, knots=myknots.t),omit.missing=TRUE)
		newdata.t <- data.frame(x1t=rep(t.fit, each=N.fit),x2t=rep(t.fit, N.fit))
		pred.t <- predict(fit.t, newdata.t)
		detach(pos=4)

		s.g <- matrix(pred.t, N.fit)
		G <- ( s.g + t(s.g) )/2

		eigenDecomp <- eigen(G)

		par(mfrow=c(2,3))
		for(i in 1:6) {
			plot(eigenDecomp$vectors[,i], type="l")
		}


		temp.data <- data.frame(var.y = as.vector(resd^2),  var.x = as.vector(tx) )
		attach(temp.data, pos=6)
		fit.var <- spm( var.y ~ f(var.x,knots=t.knot) ,omit.missing=TRUE)
		newx <- data.frame(var.x=t.fit)
		var.fit <- predict(fit.var,newx)          
		detach(pos=6)

		var.noise <- mean( var.fit - diag(matrix(pred.t,N.fit)) ) 


		
		###############################################
		if(var.noise<0){
		timesX2[b]<-var.noise
		var.noise=.002
		}
    	###############################################
		
    
    
	    fpca1.value <- eigenDecomp$values* tlength / N.fit

    	fpca1.value <- ifelse(fpca1.value>=0, fpca1.value, 0)

    
		K1=min(35, N)
		
		fpca1.vectors <- eigenDecomp$vectors[, 1:K1]*sqrt(N.fit/tlength)


		###     The eigenfunctions are unique only up to a change of signs.
		###     Select the signs of eigenfunctions so that the integration over the domain 
		###     is non-negative
    	for(i2 in 1:K1) {
	        v2 <- fpca1.vectors[,i2]
        	tempsign <- sum(v2)
    	    fpca1.vectors[,i2] <- ifelse(tempsign<0, -1,1) * v2
	    }


	
		lam1=fpca1.value[1:K1]
		phi1=fpca1.vectors



		###
		###     Estimatimate principal component scores and predict various functions
		###

		library(MASS) 

		index <- matrix(0, M, N)
		mu.hat <- matrix(0, M, N)    
		y.center <- matrix(0, M, N)

		phi1.hat <- array(0, dim=c( M, N, K1) )

		t.fit.matrix <- matrix(rep(t.fit, each=N), nrow=N)

		for(m in 1:M ) {
    	    row.ind <- m
        	index <- apply( (t.fit.matrix - tx[row.ind, ])^2, 1, which.min )
    	    y.center[row.ind, ] <- y[ row.ind ,] - mu.fit[ index ]
	        phi1.hat[row.ind, , ] <- phi1[index, ]
		}



		est.si1 <- matrix(0, nrow=M, ncol=K1)

		est.si1.sd <- matrix(0, nrow=M, ncol=K1)


		y.hat.center <- matrix(0, M, N.fit) 
		y.hat <- matrix(0, M, N.fit)
		y.hat.subject <- matrix(0, M, N.fit)

		y.hat.sd <- matrix(0, M, N.fit)
		y.hat.subject.sd <- matrix(0, M, N.fit)



		for(m in 1:M) {
			cov.y <- matrix(0, N, N)
			cov.xi.y <- matrix(0, K1,  N )

			ind1 <- m
			cov.y [ 1:N , 1:N ] <- phi1.hat[ind1,,] %*% diag( lam1) %*%  t( phi1.hat[ind1,,] )  + diag( rep(var.noise, N) )

			ind1 <- m
			cov.xi.y[1:K1, 1:N ] <- diag( lam1) %*% t( phi1.hat[ind1,,] )


			temp.mat <- cov.xi.y %*% ginv(cov.y)
			score <- temp.mat %*% as.vector( t(y.center[m,]) ) 
			var.score <- diag( c(lam1 ) ) -  temp.mat %*% t(cov.xi.y)

			est.si1[m ,] <- score[1:K1]

			y.hat.subject[m, ] <- mu.fit + phi1 %*% est.si1[m,] 
			y.hat.subject.sd[m,] <-  sqrt( diag( phi1 %*% var.score[1:K1, 1:K1] %*% t( phi1 ) ) )

			row.ind <- m
			y.hat.center[ row.ind ,] <-  phi1 %*% est.si1[m,]
			y.hat[row.ind,] <- y.hat.center[ row.ind ,]  +  mu.fit

			temp1 <- cbind( phi1)
			temp2 <- var.score[ c(1:K1) , c(1:K1) ]
			y.hat.sd[row.ind, ] <- sqrt( diag( temp1 %*% temp2 %*% t(temp1) ) )   
		}
			
			SX2<-y.hat

		###############################################
		# end sparse data estimation procedure
		###############################################




################## plot
par(mfrow=c(1,2))
## X2
range_ylim =  c(-2.5,1.5)
power10 = floor(range_ylim[1])
axisx_range = seq(r[1]-0.1, r[length(r)], by = 1) 
axisx_label =axisx_range
axisy_range =  seq(from=range_ylim[1], to=range_ylim[2], length.out=5)
axisy_label = round(axisy_range,1)


for(j in 2:3){
	
plot(r,SX2[j,],type="l",col="red",ylim=c(-2.5,1.5),ylab="",main="",lwd=2,axes=F,cex.lab=1.5);
   axis(1,at=axisx_range ,col="grey35",labels=axisx_label, cex=1.2, font=2.5,cex.axis=1.3) 
   axis(2,at=axisy_range ,col="grey35",labels=axisy_label,cex=1.2, font=2.5,cex.axis=1.3)
lines(tx[j,],y[j,],type="p",col="red",pch=21, cex=1);
lines(r,X2[j,],type="l",col="black",lty=3,lwd=2)
}

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




#############################################################
#########       PLOTS OF Y, X_i1, X_i2         ##############
#############################################################


clrs <- col2rgb("gray")
clrs <- rgb(t(clrs)/255, alpha=.25)


par(mfrow=c(1,3))
plot(t,outcomes[1,],type="l",ylim=range(outcomes),ylab="",main=expression(Y(t)),col=clrs)
for(i in 1:n){
	lines(t,outcomes[i,],col=clrs)
	}

for (i in 1:3){lines(t, outcomes[i,], col=i+1, lwd=1) }

	
plot(s,X1[1,],type="l",ylab="",ylim=range(X1),main=expression(X[i1](s)),col=clrs)
for(i in 1:n){
	lines(s,X1[i,],col=clrs)
	}			
		 
for (i in 1:3){lines(s, X1[i,], col=i+1, lwd=1);lines(s, V1[i,], col=i+1, lwd=1, lty=2)  }

plot(r,X2[1,],type="l",ylab="",ylim=range(X2),main=expression(X[i2](r)),col=clrs)
for(i in 1:n){
	lines(r,X2[i,],col=clrs)
	}	
for (i in 1:3){lines(r, X2[i,], col=i+1, lwd=1); lines(r, V2[i,], col=i+1, lwd=1, lty=2) }	


################################################################
##                         pffr                               ##
################################################################
	
	
yvec <- as.vector(t(outcomes))
 

tngrid=length(t)
sngrid=length(s)
rngrid=length(r)

### fit the model ###
tmat1 <- matrix(t, nrow=n*tngrid, nc=sngrid,byrow=FALSE)
tmat2 <- matrix(t, nrow=n*tngrid, nc=rngrid,byrow=FALSE)
smat <- matrix(s, nrow=n*tngrid, nc=sngrid, byrow=TRUE)
rmat <- matrix(r, nrow=n*tngrid, nc=rngrid, byrow=TRUE)

LSX1 <- L1*SX1
LSX2 <- L2*SX2

Lmat1 <- LSX1[rep(1:n,each=tngrid),]
Lmat2 <- LSX2[rep(1:n,each=tngrid),]

tvec <- matrix(t, nrow=n*tngrid, nc=1,byrow=FALSE)	


m100 <- bam(yvec ~s(tvec,bs="ps",k=12)+te(tmat1,smat,by=Lmat1, bs="ps")+te(tmat2,rmat,by=Lmat2, bs="ps",k=7),method="REML")



################################################################
##              for  METHOD 1    pffr                         ##
################################################################
	
###### obtain INTERCEPT 

term0 <- m100$smooth[[1]]

x=t

predDat <- data.frame(x)

names(predDat) <- c(term0$term)

PX0 <- PredictMat(term0, dat=predDat)


m100_hatBeta0_func<- matrix(PX0%*%m100$coefficients[term0$first.para:term0$last.para], ncol=length(t))

m100_hatBeta0[b,]<-m100_hatBeta0_func+m100$coefficients[1]



###### obtain BETA1

term1 <- m100$smooth[[2]]

x=rep(t, rep(length(s), length(t)))

predDat <- data.frame(x=rep(t, rep(length(s), length(t))), y=rep(s, length(t)), by=rep(1, length(t)*length(x)))

names(predDat) <- c(term1$term, term1$by)

PX <- PredictMat(term1, dat=predDat)


m100_hatBeta1_sparse[b,,] <- matrix(PX%*%m100$coefficients[term1$first.para:term1$last.para], nrow=length(s), ncol=length(t))

#######################################################################
######################### Pointwise CI for trueBeta0(t)  ##################
#######################################################################


Sigma_Beta0<-m100$Vp[term0$first.para:term0$last.para,term0$first.para:term0$last.para]
PX0<-as.matrix(PX0)

varBeta0<-matrix(0,nrow=length(t),ncol=1)
lBoundBeta0<-matrix(0,nrow=length(t),ncol=1)
uBoundBeta0<-matrix(0,nrow=length(t),ncol=1)

for(i in 1:(length(t))){
varBeta0[i]<-t(PX0[i,])%*%Sigma_Beta0%*%PX0[i,]
lBoundBeta0[i] <- as.vector(m100_hatBeta0[b,])[i] - 1.96*sqrt(varBeta0[i])
uBoundBeta0[i] <- as.vector(m100_hatBeta0[b,])[i] + 1.96*sqrt(varBeta0[i])
}

#######################################################################
######################### Pointwise CI for beta1(s,t)  ################
#######################################################################


Sigma_Beta1<-m100$Vp[term1$first.para:term1$last.para,term1$first.para:term1$last.para]
PX<-as.matrix(PX)

varBeta1<-matrix(0,nrow=length(s)*length(t),ncol=1)
lBoundBeta1<-matrix(0,nrow=length(s)*length(t),ncol=1)
uBoundBeta1<-matrix(0,nrow=length(s)*length(t),ncol=1)

for(i in 1:(length(s)*length(t))){
varBeta1[i]<-t(PX[i,])%*%Sigma_Beta1%*%PX[i,]
lBoundBeta1[i] <- as.vector(m100_hatBeta1_sparse[b,,])[i] - 1.96*sqrt(varBeta1[i])
uBoundBeta1[i] <- as.vector(m100_hatBeta1_sparse[b,,])[i] + 1.96*sqrt(varBeta1[i])
}
varBeta1<-matrix(varBeta1,nrow=length(s),ncol=length(t))
lBoundBeta1<-matrix(lBoundBeta1,nrow=length(s),ncol=length(t))
uBoundBeta1<-matrix(uBoundBeta1,nrow=length(s),ncol=length(t))



#########----------------    beta2(s,t) ---------------############



term2 <- m100$smooth[[3]]

x=rep(t, rep(length(r), length(t)))

predDat <- data.frame(x=rep(t, rep(length(r), length(t))), y=rep(r, length(t)), by=rep(1, length(t)*length(x)))

names(predDat) <- c(term2$term, term2$by)

PX <- PredictMat(term2, dat=predDat)


m100_hatBeta2_sparse[b,,] <- matrix(PX%*%m100$coefficients[term2$first.para:term2$last.para], nrow=length(r), ncol=length(t))


#######################################################################
######################### Pointwise CI for beta2(s,t)  ################
#######################################################################


Sigma_Beta2<-m100$Vp[term2$first.para:term2$last.para,term2$first.para:term2$last.para]
PX<-as.matrix(PX)

varBeta2<-matrix(0,nrow=length(r)*length(t),ncol=1)
lBoundBeta2<-matrix(0,nrow=length(r)*length(t),ncol=1)
uBoundBeta2<-matrix(0,nrow=length(r)*length(t),ncol=1)

for(i in 1:(length(r)*length(t))){
varBeta2[i]<-t(PX[i,])%*%Sigma_Beta2%*%PX[i,]
lBoundBeta2[i] <- as.vector(m100_hatBeta2_sparse[b,,])[i] - 1.96*sqrt(varBeta2[i])
uBoundBeta2[i] <- as.vector(m100_hatBeta2_sparse[b,,])[i] + 1.96*sqrt(varBeta2[i])
}
varBeta2<-matrix(varBeta2,nrow=length(r),ncol=length(t))
lBoundBeta2<-matrix(lBoundBeta2,nrow=length(r),ncol=length(t))
uBoundBeta2<-matrix(uBoundBeta2,nrow=length(r),ncol=length(t))


####################### Coverage and Width ########################

##------------------------- trueBeta0(t) -------------------------##


for(j in 1:(length(t))){
if(lBoundBeta0[j]<=trueBeta0[j]) if(uBoundBeta0[j]>=trueBeta0[j])
covBeta0[b,j]<-1
}

	
ICBeta0[b]<-mean(covBeta0[b,])	



WI0<-matrix(0,nrow=1,ncol=length(t))	

for(j in 1:(length(t))){
		
WI0[j]<-uBoundBeta0[j]-lBoundBeta0[j]
	
}
	
IWBeta0[b]<-mean(WI0)	#mean cov over all s and t



##------------------------- beta1(s,t)------------------------##


for(i in 1:(length(s))){
for(j in 1:(length(t))){
		
if(lBoundBeta1[i,j]<=trueBeta1[i,j]) if(uBoundBeta1[i,j]>=trueBeta1[i,j])
covBeta1[b,i,j]<-1
	
}
}	

ICBeta1[b]<-mean(covBeta1[b,,])	#mean cov over all s and t



WI1<-matrix(0,nrow=length(s),ncol=length(t))	#cov at each s and t

for(i in 1:(length(s))){
for(j in 1:(length(t))){
		
WI1[i,j]<-uBoundBeta1[i,j]-lBoundBeta1[i,j]
	
}
}
	
IWBeta1[b]<-mean(WI1)	#mean cov over all s and t


##------------------------- beta2(s,t)------------------------##


for(i in 1:(length(r))){
for(j in 1:(length(t))){

if(lBoundBeta2[i,j]<=trueBeta2[i,j]) if(uBoundBeta2[i,j]>=trueBeta2[i,j])
		
covBeta2[b,i,j]<-1
	
}
}	

ICBeta2[b]<-mean(covBeta2[b,,])	#mean cov over all s and t

WI2<-matrix(0,nrow=length(r),ncol=length(t))	#cov at each s and t

for(i in 1:(length(r))){
for(j in 1:(length(t))){
		
WI2[i,j]<-uBoundBeta2[i,j]-lBoundBeta2[i,j]
	
}
}
	
IWBeta2[b]<-mean(WI2)	#mean cov over all s and t


###############################################################################
###############################################################################
###############################################################################
################# pfr approach Goldsmith et al. (2011) ########################
###############################################################################
###############################################################################
###############################################################################


betaHat.MULT_Beta0<-matrix(nrow=1, ncol=length(t))
betaHat.MULT_Beta0_lb<-matrix(nrow=1, ncol=length(t))
betaHat.MULT_Beta0_ub<-matrix(nrow=1, ncol=length(t))
betaHat.MULT<-matrix(nrow=length(s) +length(r), ncol=length(t))
betaHat.MULT_lb<-matrix(nrow=length(s) +length(r), ncol=length(t))
betaHat.MULT_ub<-matrix(nrow=length(s) +length(r), ncol=length(t))
fitted_pfr<-matrix(nrow=n,ncol=length(t))

for(i in 1:(length(t))){
		fit_pfr<- pfr(Y=outcomes[,i], funcs=cbind(LSX1,LSX2), kz=30, kb=30)
		betaHat.MULT_Beta0[i]<-fit_pfr$beta.covariates
		betaHat.MULT[,i] <- fit_pfr$betaHat
		betaHat.MULT_ub[,i] <- fit_pfr$Bound[,1]
		betaHat.MULT_lb[,i] <- fit_pfr$Bound[,2]
		fitted_pfr[,i]<-fit_pfr$fitted.vals
		
		# save the gam fit in order to get access to Vp needed for CI for beta0
		gam_pfr<-fit_pfr$gam
		betaHat.MULT_Beta0_lb[i]<-fit_pfr$beta.covariates[1]-1.96*sqrt(gam_pfr$Vp[1,1])
		betaHat.MULT_Beta0_ub[i]<-fit_pfr$beta.covariates[1]+1.96*sqrt(gam_pfr$Vp[1,1])

} #end loop over t



### smooth Beta0



MULT_Beta0<-gam(as.vector(betaHat.MULT_Beta0)~s(t),method="REML")
term0 <- MULT_Beta0$smooth[[1]]
x=t
predDat <- data.frame(x)
names(predDat) <- c(term0$term)
PX0 <- PredictMat(term0, dat=predDat)
MULT_Beta0_hatBeta0<- matrix(PX0%*%MULT_Beta0$coefficients[term0$first.para:term0$last.para], ncol=length(t))
betaHat0.mult[b,]<-MULT_Beta0_hatBeta0+MULT_Beta0$coefficients[1]


### smooth lb for Beta0

MULT_Beta0<-gam(as.vector(betaHat.MULT_Beta0_lb)~s(t),method="REML")
term0 <- MULT_Beta0$smooth[[1]]
x=t
predDat <- data.frame(x)
names(predDat) <- c(term0$term)
PX0 <- PredictMat(term0, dat=predDat)
MULT_Beta0_hatBeta0_lb<- matrix(PX0%*%MULT_Beta0$coefficients[term0$first.para:term0$last.para], ncol=length(t))
betaHat0.mult_lb<-MULT_Beta0_hatBeta0_lb+MULT_Beta0$coefficients[1]

### smooth ub for Beta0

MULT_Beta0<-gam(as.vector(betaHat.MULT_Beta0_ub)~s(t),method="REML")
term0 <- MULT_Beta0$smooth[[1]]
x=t
predDat <- data.frame(x)
names(predDat) <- c(term0$term)
PX0 <- PredictMat(term0, dat=predDat)
MULT_Beta0_hatBeta0_ub<- matrix(PX0%*%MULT_Beta0$coefficients[term0$first.para:term0$last.para], ncol=length(t))
betaHat0.mult_ub<-MULT_Beta0_hatBeta0_ub+MULT_Beta0$coefficients[1]




#### Beta1
long_s<-rep(s,length(t))
t_mat<-matrix(0,nrow=length(t),ncol=length(s))
for(i in 1:(length(t))){
t_mat[i,]<-rep(t[i],length(s))		
}
long_t<-as.vector(t(t_mat))

MULT_te1<-gam(as.vector(betaHat.MULT[(1:length(s)),])~te(long_t,long_s),method="REML")

term_pfr_beta1 <- MULT_te1$smooth[[1]]

x=rep(t, rep(length(s), length(t)))

predDat <- data.frame(x=rep(t, rep(length(s), length(t))), y=rep(s, length(t)), by=rep(1, length(t)*length(x)))

names(predDat) <- c(term_pfr_beta1$term, term_pfr_beta1$by)

PX_pfr_beta1 <- PredictMat(term_pfr_beta1, dat=predDat)


betaHat1.mult[b,,] <- matrix(PX_pfr_beta1%*%MULT_te1$coefficients[term_pfr_beta1$first.para:term_pfr_beta1$last.para]+MULT_te1$coefficients[1], nrow=length(s), ncol=length(t))

### Beta1 lower bound

MULT_te1<-gam(as.vector(betaHat.MULT_lb[(1:length(s)),])~te(long_t,long_s),method="REML")

term_pfr_beta1 <- MULT_te1$smooth[[1]]

x=rep(t, rep(length(s), length(t)))

predDat <- data.frame(x=rep(t, rep(length(s), length(t))), y=rep(s, length(t)), by=rep(1, length(t)*length(x)))

names(predDat) <- c(term_pfr_beta1$term, term_pfr_beta1$by)

PX_pfr_beta1 <- PredictMat(term_pfr_beta1, dat=predDat)


betaHat1.mult_lb <- matrix(PX_pfr_beta1%*%MULT_te1$coefficients[term_pfr_beta1$first.para:term_pfr_beta1$last.para]+MULT_te1$coefficients[1], nrow=length(s), ncol=length(t))

### Beta1 upper bound

MULT_te1<-gam(as.vector(betaHat.MULT_ub[(1:length(s)),])~te(long_t,long_s),method="REML")

term_pfr_beta1 <- MULT_te1$smooth[[1]]

x=rep(t, rep(length(s), length(t)))

predDat <- data.frame(x=rep(t, rep(length(s), length(t))), y=rep(s, length(t)), by=rep(1, length(t)*length(x)))

names(predDat) <- c(term_pfr_beta1$term, term_pfr_beta1$by)

PX_pfr_beta1 <- PredictMat(term_pfr_beta1, dat=predDat)


betaHat1.mult_ub <- matrix(PX_pfr_beta1%*%MULT_te1$coefficients[term_pfr_beta1$first.para:term_pfr_beta1$last.para]+MULT_te1$coefficients[1], nrow=length(s), ncol=length(t))


#### Beta2


long_r<-rep(r,length(t))
t_mat_r<-matrix(0,nrow=length(t),ncol=length(r))
for(i in 1:(length(t))){
t_mat_r[i,]<-rep(t[i],length(r))		
}
long_t_r<-as.vector(t(t_mat_r))


MULT_te2<-gam(as.vector(betaHat.MULT[(length(s)+1):(length(s) +length(r)),])~te(long_t_r,long_r),method="REML")


term_pfr_beta2 <- MULT_te2$smooth[[1]]

x=rep(t, rep(length(r), length(t)))

predDat <- data.frame(x=rep(t, rep(length(r), length(t))), y=rep(r, length(t)), by=rep(1, length(t)*length(x)))

names(predDat) <- c(term_pfr_beta2$term, term_pfr_beta2$by)

PX_pfr_beta2 <- PredictMat(term_pfr_beta2, dat=predDat)


betaHat2.mult[b,,] <- matrix(PX_pfr_beta2%*%MULT_te2$coefficients[term_pfr_beta2$first.para:term_pfr_beta2$last.para]+MULT_te2$coefficients[1], nrow=length(r), ncol=length(t))


### Beta2 lower bound

MULT_te2<-gam(as.vector(betaHat.MULT_lb[(length(s)+1):(length(s) +length(r)),])~te(long_t_r,long_r),method="REML")


term_pfr_beta2 <- MULT_te2$smooth[[1]]

x=rep(t, rep(length(r), length(t)))

predDat <- data.frame(x=rep(t, rep(length(r), length(t))), y=rep(r, length(t)), by=rep(1, length(t)*length(x)))

names(predDat) <- c(term_pfr_beta2$term, term_pfr_beta2$by)

PX_pfr_beta2 <- PredictMat(term_pfr_beta2, dat=predDat)


betaHat2.mult_lb <- matrix(PX_pfr_beta2%*%MULT_te2$coefficients[term_pfr_beta2$first.para:term_pfr_beta2$last.para]+MULT_te2$coefficients[1], nrow=length(r), ncol=length(t))

### Beta 2 upper bound

MULT_te2<-gam(as.vector(betaHat.MULT_ub[(length(s)+1):(length(s) +length(r)),])~te(long_t_r,long_r),method="REML")


term_pfr_beta2 <- MULT_te2$smooth[[1]]

x=rep(t, rep(length(r), length(t)))

predDat <- data.frame(x=rep(t, rep(length(r), length(t))), y=rep(r, length(t)), by=rep(1, length(t)*length(x)))

names(predDat) <- c(term_pfr_beta2$term, term_pfr_beta2$by)

PX_pfr_beta2 <- PredictMat(term_pfr_beta2, dat=predDat)


betaHat2.mult_ub <- matrix(PX_pfr_beta2%*%MULT_te2$coefficients[term_pfr_beta2$first.para:term_pfr_beta2$last.para]+MULT_te2$coefficients[1], nrow=length(r), ncol=length(t))




##------------------------- Beta0(t) -------------------------##


covBeta0_pfr<-matrix(0,nrow=1,ncol=length(t))


for(j in 1:(length(t))){
if(betaHat0.mult_lb[j]<=trueBeta0[j]) if(betaHat0.mult_ub[j]>=trueBeta0[j])		
covBeta0_pfr[,j]<-1
}
	
ICBeta0_pfr[b]<-mean(covBeta0_pfr)	



WI0_pfr<-matrix(0,nrow=1,ncol=length(t))	

for(j in 1:(length(t))){
		
WI0_pfr[j]<-betaHat0.mult_ub[j]-betaHat0.mult_lb[j]
	
}
	
IWBeta0_pfr[b]<-mean(WI0_pfr)	#mean cov over all s and t



##------------------------- beta1(s,t)------------------------##

covBeta1_pfr<-matrix(0,nrow=length(s),ncol=length(t))
covBeta2_pfr<-matrix(0,nrow=length(r),ncol=length(t))

for(i in 1:(length(s))){
for(j in 1:(length(t))){
if(betaHat1.mult_lb[i,j]<=trueBeta1[i,j]) if(betaHat1.mult_ub[i,j]>=trueBeta1[i,j])	

covBeta1_pfr[i,j]<-1

}
}	

ICBeta1_pfr[b]<-mean(covBeta1_pfr)	#mean cov over all s and t

WI1_pfr<-matrix(0,nrow=length(s),ncol=length(t))	#cov at each s and t

for(i in 1:(length(s))){
for(j in 1:(length(t))){
		
WI1_pfr[i,j]<-betaHat1.mult_ub[i,j]-betaHat1.mult_lb[i,j]
	
}
}
	
IWBeta1_pfr[b]<-mean(WI1_pfr)	#mean cov over all s and t


##------------------------- beta2(s,t)------------------------##

for(i in 1:(length(r))){
for(j in 1:(length(t))){
if(betaHat2.mult_lb[i,j]<=trueBeta2[i,j]) if(betaHat2.mult_ub[i,j]>=trueBeta2[i,j])	
covBeta2_pfr[i,j]<-1	
}
}	

ICBeta2_pfr[b]<-mean(covBeta2_pfr)	#mean cov over all s and t

WI2_pfr<-matrix(0,nrow=length(r),ncol=length(t))	#cov at each s and t

for(i in 1:(length(r))){
for(j in 1:(length(t))){
		
WI2_pfr[i,j]<-betaHat2.mult_ub[i,j]-betaHat2.mult_lb[i,j]
	
}
}
	
	
IWBeta2_pfr[b]<-mean(WI2_pfr)	#mean cov over all s and t


######################### EMSE, BIAS2, VAR ####################################

ISE_Beta0[b]<-mean((m100_hatBeta0[b,]-trueBeta0)^2)
ISE_Beta1[b]<-mean((m100_hatBeta1_sparse[b,,]-trueBeta1)^2)
ISE_Beta2[b]<-mean((m100_hatBeta2_sparse[b,,]-trueBeta2)^2)


ISE_Beta0_pfr[b]<-mean((betaHat0.mult[b,]-trueBeta0)^2)
ISE_Beta1_pfr[b]<-mean((betaHat1.mult[b,,]-trueBeta1)^2)
ISE_Beta2_pfr[b]<-mean((betaHat2.mult[b,,]-trueBeta2)^2)



############################ R2 for pffr ######################################

meanoutcomes<-colSums(outcomes)/n

intvec<-rep(meanoutcomes,n)

yvec <- as.vector(t(outcomes))

IR2[b]<-1-crossprod(fitted(m100)-yvec)/crossprod(yvec-intvec)


############################ R2 for pfr ######################################


pfrBeta0_mat<-matrix(betaHat0.mult[b,],nrow=n,ncol=length(t),byrow=T)



all_pfr <- pfrBeta0_mat + LSX1 %*% betaHat1.mult[b,,] + LSX2 %*% betaHat2.mult[b,,]


pfr_fitted_outcomes <-  all_pfr 



yfitted_pfr <- as.vector(t(pfr_fitted_outcomes))


IR2_pfr[b]<-1-crossprod(yfitted_pfr-yvec)/crossprod(yvec-intvec)



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

}								#end B simulations#
#################################################################################
#################################################################################
#################################################################################
#################################################################################


MIR2<-mean(IR2)
MIR2_pfr<-mean(IR2_pfr)

########################## Beta 0 ################################
MISE_Beta0<-mean(ISE_Beta0)			# E over B simulations

##### ---------------------- VAR --------------------------#######

E_m100_hatBeta0<-matrix(0,nrow=1,ncol=length(t))	# E over B sims
for(j in 1:(length(t))){
E_m100_hatBeta0[1,j]<-mean(m100_hatBeta0[,j])
}


IVAR_Beta0<-matrix(nrow=B, ncol=length(t))

for(b in 1:B){
for(j in 1:(length(t))){
IVAR_Beta0[b,j]<-(m100_hatBeta0[b,j]-E_m100_hatBeta0[1,j])^2
}
}


MIVAR_Beta0<-mean(IVAR_Beta0)	#mean over t and B

##### ---------------------- BIAS -------------------------#######

MIBIAS_Beta0<-mean((E_m100_hatBeta0-trueBeta0)^2)



########################## Beta 1 ################################
MISE_Beta1<-mean(ISE_Beta1)

##### ---------------------- VAR --------------------------#######

E_m100_hatBeta1_sparse<-matrix(0,nrow=length(s),ncol=length(t))
for(i in 1:(length(s))){
for(j in 1:(length(t))){
E_m100_hatBeta1_sparse[i,j]<-mean(m100_hatBeta1_sparse[,i,j])
}
}


IVAR_Beta1<-array(0, dim=c(B,nrow=length(s), ncol=length(t)))

for(b in 1:B){
for(i in 1:(length(s))){
for(j in 1:(length(t))){
IVAR_Beta1[b,i,j]<-(m100_hatBeta1_sparse[b,i,j]-E_m100_hatBeta1_sparse[i,j])^2
}
}
}

MIVAR_Beta1<-mean(IVAR_Beta1)

##### ---------------------- BIAS -------------------------#######

MIBIAS_Beta1<-mean((E_m100_hatBeta1_sparse-trueBeta1)^2)

########################## Beta 2 ################################
MISE_Beta2<-mean(ISE_Beta2)

##### ---------------------- VAR --------------------------#######

E_m100_hatBeta2_sparse<-matrix(0,nrow=length(r),ncol=length(t))
for(i in 1:(length(r))){
for(j in 1:(length(t))){
E_m100_hatBeta2_sparse[i,j]<-mean(m100_hatBeta2_sparse[,i,j])
}
}


IVAR_Beta2<-array(0, dim=c(B,nrow=length(r), ncol=length(t)))

for(b in 1:B){
for(i in 1:(length(r))){
for(j in 1:(length(t))){
IVAR_Beta2[b,i,j]<-(m100_hatBeta2_sparse[b,i,j]-E_m100_hatBeta2_sparse[i,j])^2
}
}
}

MIVAR_Beta2<-mean(IVAR_Beta2)

##### ---------------------- BIAS -------------------------#######

MIBIAS_Beta2<-mean((E_m100_hatBeta2_sparse-trueBeta2)^2)






########################## Beta 0 pfr ################################
MISE_Beta0_pfr<-mean(ISE_Beta0_pfr)			# E over B simulations

##### ---------------------- VAR --------------------------#######

E_pfr_hatBeta0<-matrix(0,nrow=1,ncol=length(t))	# E over B sims
for(j in 1:(length(t))){
E_pfr_hatBeta0[1,j]<-mean(betaHat0.mult[,j])
}


IVAR_Beta0_pfr<-matrix(nrow=B, ncol=length(t))

for(b in 1:B){
for(j in 1:(length(t))){
IVAR_Beta0_pfr[b,j]<-(betaHat0.mult[b,j]-E_pfr_hatBeta0[1,j])^2
}
}


MIVAR_Beta0_pfr<-mean(IVAR_Beta0_pfr)	#mean over t and B

##### ---------------------- BIAS -------------------------#######

MIBIAS_Beta0_pfr<-mean((E_pfr_hatBeta0-trueBeta0)^2)





########################## Beta 1 pfr ################################
MISE_Beta1_pfr<-mean(ISE_Beta1_pfr)

##### ---------------------- VAR --------------------------#######

E_pfr_hatBeta1<-matrix(0,nrow=length(s),ncol=length(t))
for(i in 1:(length(s))){
for(j in 1:(length(t))){
E_pfr_hatBeta1[i,j]<-mean(betaHat1.mult[,i,j])
}
}


IVAR_Beta1_pfr<-array(0, dim=c(B,nrow=length(s), ncol=length(t)))

for(b in 1:B){
for(i in 1:(length(s))){
for(j in 1:(length(t))){
IVAR_Beta1_pfr[b,i,j]<-(betaHat1.mult[b,i,j]-E_pfr_hatBeta1[i,j])^2
}
}
}

MIVAR_Beta1_pfr<-mean(IVAR_Beta1_pfr)

##### ---------------------- BIAS -------------------------#######

MIBIAS_Beta1_pfr<-mean((E_pfr_hatBeta1-trueBeta1)^2)


########################## Beta 2 pfr ################################

MISE_Beta2_pfr<-mean(ISE_Beta2_pfr)

##### ---------------------- VAR --------------------------#######

E_pfr_hatBeta2<-matrix(0,nrow=length(r),ncol=length(t))
for(i in 1:(length(r))){
for(j in 1:(length(t))){
E_pfr_hatBeta2[i,j]<-mean(betaHat2.mult[,i,j])
}
}


IVAR_Beta2_pfr<-array(0, dim=c(B,nrow=length(r), ncol=length(t)))

for(b in 1:B){
for(i in 1:(length(r))){
for(j in 1:(length(t))){
IVAR_Beta2_pfr[b,i,j]<-(betaHat2.mult[b,i,j]-E_pfr_hatBeta2[i,j])^2
}
}
}

MIVAR_Beta2_pfr<-mean(IVAR_Beta2_pfr)

##### ---------------------- BIAS -------------------------#######

MIBIAS_Beta2_pfr<-mean((E_pfr_hatBeta2-trueBeta2)^2)

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


McovBeta0<-matrix(nrow=1,ncol=length(t))

for(j in 1:(length(t))){
McovBeta0[1,j]<-mean(covBeta0[,j])
}

McovBeta1<-matrix(nrow=length(s),ncol=length(t))

for(i in 1:(length(s))){
for(j in 1:(length(t))){
McovBeta1[i,j]<-mean(covBeta1[,i,j])
}
}

McovBeta2<-matrix(nrow=length(r),ncol=length(t))

for(i in 1:(length(r))){
for(j in 1:(length(t))){
McovBeta2[i,j]<-mean(covBeta2[,i,j])
}
}


#### pfr

McovBeta1_pfr<-matrix(nrow=length(s),ncol=length(t))

for(i in 1:(length(s))){
for(j in 1:(length(t))){
McovBeta1_pfr[i,j]<-mean(covBeta1_pfr[i,j])
}
}

McovBeta2_pfr<-matrix(nrow=length(r),ncol=length(t))

for(i in 1:(length(r))){
for(j in 1:(length(t))){
McovBeta2_pfr[i,j]<-mean(covBeta2_pfr[i,j])
}
}

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

IACBeta0<-mean(ICBeta0)
IACBeta1<-mean(ICBeta1)
IACBeta2<-mean(ICBeta2)


IACBeta0_pfr<-mean(ICBeta0_pfr)
IACBeta1_pfr<-mean(ICBeta1_pfr)
IACBeta2_pfr<-mean(ICBeta2_pfr)

IAWBeta0<-mean(IWBeta0)
IAWBeta1<-mean(IWBeta1)
IAWBeta2<-mean(IWBeta2)


IAWBeta0_pfr<-mean(IWBeta0_pfr)
IAWBeta1_pfr<-mean(IWBeta1_pfr)
IAWBeta2_pfr<-mean(IWBeta2_pfr)

##################################################################
Results<-matrix(c(MISE_Beta0,MIVAR_Beta0,MIBIAS_Beta0,
MISE_Beta1,MIVAR_Beta1,MIBIAS_Beta1,
MISE_Beta2,MIVAR_Beta2,MIBIAS_Beta2,
MISE_Beta0_pfr,MIVAR_Beta0_pfr,MIBIAS_Beta0_pfr,
MISE_Beta1_pfr,MIVAR_Beta1_pfr,MIBIAS_Beta1_pfr,
MISE_Beta2_pfr,MIVAR_Beta2_pfr,MIBIAS_Beta2_pfr),nrow=6,ncol=3,byrow=TRUE)

colnames(Results)<-c("MISE","MIVAR","MIBIAS")

rownames(Results)<-c("Beta0","Beta1","Beta2","Beta0_pfr","Beta1_pfr","Beta2_pfr")


CI_Results<-matrix(c(IACBeta0,IACBeta1,IACBeta2, IAWBeta0,IAWBeta1,IAWBeta2,IACBeta0_pfr, IACBeta1_pfr,IACBeta2_pfr,IAWBeta0_pfr,IAWBeta1_pfr,IAWBeta2_pfr),nrow=3,ncol=4,byrow=FALSE)

colnames(CI_Results)<-c("IAC","IAW","IAC_pfr","IAW_pfr")

rownames(CI_Results)<-c("Beta0","Beta1","Beta2")



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

round(sqrt(Results)*1000,2)
round(CI_Results,2)

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