% ============
% Description:
% ============
%
%             This is the main program for Functional Quasi-Likelihood
%             Regression. The response is functional (densely observed) and 
%             the predictor is a vector. PACE is first used to  obtain
%             principal component scores and then SPQR (semiparametric  
%             quasi-likelihood regression) is run to fit their dependence 
%             on the predictor. The details are specified in Chiou and M\"uller (2003). 
%             
%             
%             
%
% ======
% Usage:
% ======
%
% function [ X, W ] = FQR( y, t, z, kernel,bw,linktyp,vartyp,theta,fig, p)
% 
% Input:
%
%      y:          1*n cell array, y{i} is the vector of measurements for the ith subject,
%                  i=1,...,n.
%      t:          1*n cell array, t{i} is the vector of time points for the ith subject on which
%                  corresponding measurements y{i} are taken, i=1,...,n.
%      z:          n*p matrix, z(:,j) is the jth predictor for the ith
%                  subject
%      p:          a struct obtained from setOptions.m sets the rest of arguments for PCA.m
%                  ex:
%                  >> p = setOptions();
%
%      kernel:     kernel function to be used for smoothing link and
%                  variance functions
%                  'epan' - Epanechnikov: 0.75*(1-x^2) on [-1,1] 
%                  'rect' - rectangle:    0.5 on [-1,1]
%                  'gauss'- gaussian:     exp(-x^2/2)/sqrt(2*pi) on [-4,4]
%                  Default is 'epan'.
%      bw:         a vector of length 3 if specified, correspondingly
%                  giving the bandwidth for smoothing the link function,
%                  the derivative of link function and the variance
%                  function 
%                  GCV will be used if not specified
%      linktyp:    an integer indicating the type of link function
%                  0 - unknown; 1 - identity; 2 - power 2; 3 - log; 4 -
%                  logit; 5 - cloglog; 6 - inverse; 7 - sqrt
%                  Default is 0.
%      vartyp:     an integer indicating the variance function
%                  0 - unknown; 1 - constant; 2 - binomial; 3 - poisson; 4
%                  - gamma
%                  Default is 0.
%      theta:      a vector of length n needed only if vartyp is 2 or 4;
%                  theta(i) givs the sample size ni for the ith subject
%                  from binomial B(ni,pi) or gamma G(ni,lambdai) distribution.
%      

% Output:
%      X:         a cell array that contains all returned values from PCA.m
%                 where the last element contains the names of X

%      W:         a K*1 cell array that contains all returned values
%                 from iterGLM.m for K principal scores
%    Output in W:
%
%      eta:       1*n vector where eta(i) is the estimated f(E(Y|mu)) for
%                 the ith subject (f is the link function) 
%      mu:        1*n vector where mu(i) is the estimated g(eta) (g=f^{-1})
%                 for the ith subject
%      der:       1*n vector where der(i) is the estimated g'(eta) for the
%                 ith subject
%      sigma2:    1*n vector where sigma2(i) is the estimated var(Y|mu) for
%                 the ith subject
%      beta:      1*m vector of the regression coefficients; if the link
%                 function is unknown, beta is normalized to have norm 1
%                 with first component nonnegative (without intercept)
%      alpha:     estimated intercept (0 if link function is unknown) 
%      vbeta:     covariance matrix of beta (overdispersion is included,
%                 for non-overdispersed case, divide by output dis 
%                 below); if link function is known, this is the
%                 covariance matrix of [alpha,beta] 
%      pval:      1*m vector of p-values for testing beta=0; if link
%                 function is known, this is 1*(m+1) vector of p-values for
%                 testing [alpha,beta]=0
%      dis:       estimated Pearson overdispersion parameter (1 if link
%                 function is unknown)
%      optbw:     1*3 vector giving the final bandwidths for smoothing the
%                 link function, the derivative of link function and the
%                 variance function from GCV; NaN if the corresponding
%                 function is known. 
%      ggrid:     1*100 vector giving a dense time grid for gfct 
%      gfct:      1*100 vector giving the evaluated link function on time
%                 grid ggrid 
%      vgrid:     1*100 vector giving a dense time grid for vfct 
%      vfct:      1*100 vector giving the evaluated variance function on time
%                 grid vgrid