% Example for iterGFLM.m

%%%%%%%%%%%%%%%%% Generate data %%%%%%%%%%%%%%%%%%%%%%%%
% Predictors are balanced data simulated from 
% X_i(t)=\sum_{k=1}^4 r_{ik} phi_k(t) 
% where i=1,...,100, each predictor are measure at 20 equally-spaced time
% points in [0,1], r_{ik} are from N(0,lamba_k) with variance lambda_k
% equal to 1, 0.5, 0.25, 0.125 and phi_k, k=1,...,4 are chosen from fourier
% basis. 
% Response are generated by Y_i=binomial(1,pi) where 
% logit(pi)=1 + \int beta*X_i 
% and beta is a linear combination of phi_k.
% iterGFLM can be used to estimate beta and pi.
N=100;
M=20;
phi={};
for i=1:2
    phi{2*i-1}=@(t)sqrt(2)*sin(2*pi*i*t);
    phi{2*i}=@(t)sqrt(2)*cos(2*pi*i*t);
end
T=0:1/(M-1):1;
phibasis=[phi{1}(T);phi{2}(T);phi{3}(T);phi{4}(T)];
r=randn(4,N);
r=diag(sqrt([1,0.5,0.25,0.125]))*r;
X=r'*phibasis;
eta0=1+sqrt([1/3,1/3,1/6,1/6])*r;
beta0=sqrt([1/3,1/3,1/6,1/6])*[phi{1}(T);phi{2}(T);phi{3}(T);phi{4}(T)];
pi=1./(1+exp(-eta0));
Y=random('Binomial',ones(1,N),pi);

% fit the parametric model using logit link function and binomial variance function
linktyp=4;
vartyp=2;
% the link function and variance will be plotted by default, use fig=0 to
% cancel that and we will plot them together later 
fig=0;
[eta,mu,der,sigma2,beta,alpha,tgrid,vbeta,dis,optbw,ggrid,gfct,vgrid,vfct] = iterGFLM(Y,X,T,[],[],2,linktyp,vartyp,ones(1,N),fig,[]);

% use predGFLM to predic Y for new X
newX=phibasis'*diag(sqrt([1,0.5,0.25,0.125]))*rand(4,10);
newY=predGFLM(newX',T,beta,tgrid,alpha,ggrid,gfct);

% use devGFLM to calculate the Pearson Chi-square and deviance
[D1,D2]=devGFLM(Y,mu,sigma2,vartyp);

% fit the semiparametric model using unknown link function and unknown variance function
linktyp=0;
vartyp=0;
fig=0;
[eta2,mu2,der2,sigma22,beta2,alpha2,vbeta2,pval2,dis2,optbw2,ggrid2,gfct2,vgrid2,vfct2] = iterGFLM(Y,X,T,[],[],2,linktyp,vartyp,ones(1,N),fig,[]);

% compare the true model, the fitted parametric model and the fitted nonparametric model
% plot the link function from the three models
figure
subplot(3,1,1)
plot(eta0,Y,'.r')
hold on
plot(sort(eta0),sort(pi),'-b')
hold off
xlabel('eta')
ylabel('observed Y')
xlim([-2,2])
ylim([0,1])
title('logit link function')
subplot(3,1,2)
plot(eta,Y,'.r')
hold on
plot(ggrid,gfct,'-b')
hold off
xlabel('estimated eta')
ylabel('observed Y')
xlim([-2,2])
ylim([0,1])
title('logit link function')
subplot(3,1,3)
plot(eta2+1,Y,'.r')
hold on
plot(ggrid2+1,gfct2,'-b')
hold off
xlabel('estimated eta (recentered)')
ylabel('observed Y')
xlim([-2,2])
ylim([0,1])
title('estimated link function')
% plot the variance function from the three models
figure
subplot(3,1,1)
plot(pi,(Y-pi).^2,'.r')
hold on
tmpgrid=0:0.1:1';
plot(tmpgrid,tmpgrid.*(1-tmpgrid),'-b')
hold off
xlabel('mu')
ylabel('squared residual')
xlim([0,1])
ylim([0,1])
title('binomial variance function')
subplot(3,1,2)
plot(mu,(Y-mu).^2,'.r')
hold on
plot(vgrid,vfct,'-b')
hold off
xlabel('estimated mu')
ylabel('estimated squared residual')
xlim([0,1])
ylim([0,1])
title('binomial variance function')
subplot(3,1,3)
plot(mu2,(Y-mu2).^2,'.r')
hold on
plot(vgrid2,vfct2,'-b')
hold off
xlabel('estimated mu')
ylabel('estimated squared residual')
xlim([0,1])
ylim([0,1])
title('estimated variance function')

