### An example
### Author: Huixia Judy Wang
###            Department of Statistics, North Carolina State University
### Last Modified On: Thursday, August 12, 2010 at 16:59
### Created Date:  July 25, 2010


    ########################################
    ######## Generate a data set     #######
    ########################################
library(mvtnorm)
set.seed(12368908)
tau=0.5
nt = 13  # number of repeated times
n = 30
N = n*nt
beta=c(0,1)  # effects of z1 and z2
eta=1 # degree of variation of alpha3

# intra-subject correlation matrix
sigma = matrix(0.8, nrow=nt, ncol=nt)
diag(sigma)=1
# generate data for all subjects assuming no missing visits
No = rep(0:(nt-1),n)
time = No + runif(N, -0.5, 0.5)
# covariates that have time-varying effects
x1 = rexp(N, 1)
x2 = rnorm(N, 0, 1)
x3 = runif(N, time/10, 2+time/10)
# covariates that has constant effects
z1 = rep(rbinom(n, 1, 0.5), each=nt)
z2 = rep(rbinom(n, 1, 0.5), each=nt)
eij = rmvnorm(n, mean=rep(0,nt), sigma=sigma)
eij = as.vector(t(eij))
eij = eij - qnorm(tau)
ID = rep(1:n, each=nt)
# 20% of cases missed the scheduled visits ########################
ru=NULL
for(i in 1:n) ru=c(ru, rbinom(nt, size=1, prob=0.8))
ru = matrix(ru, byrow=T, nrow=n)
ru[,1]=1
ru = as.vector(ru)
ind = which(ru==1)
time = time[ind]
x1 = x1[ind]
x2 = x2[ind]
x3 = x3[ind]
z1 = z1[ind]
z2 = z2[ind]
eij = eij[ind]
ID = ID[ind]  # subject id
No = No[ind]  # the index of time at which the measurement was observed, varying from 0 to nt
alpha0 = 15 + 20*sin(time*pi/20)
alpha1 = -4 + eta*(20-3*time)^3/1000
alpha2 = 2 -3 * cos((3*time-25)*pi/15)
alpha3 = 6 - 0.6 *time
y = alpha0 + alpha1*x1 + alpha2*x2 + alpha3*x3 + beta[1]*z1 + beta[2]*z2 + eij
design.VC = cbind(1, x1, x2, x3) # design matrix of covariates with time-varying effects
design.C = cbind(z1, z2)   # design matrix of covariates with constant effects
dat = data.frame(y=y, time=time, design.VC=design.VC, design.C = design.C, ID=ID,No=No)


        ########################################
        #    Analyze the simulated data set    #
        ########################################
        
# load the developed functions
source("functions-PLVC.R")

# test the constancy of the varying effects
# intercept functional: alpha1(t)
tau=0.5; test.VC="VC"; test.idx=1; degree=3; se="iid"; restricted="yes"
tt=QRS.PLVC(y, time, design.VC, design.C, ID, tau=0.5, test.VC, test.idx, degree=3, se, restrcited="yes")

# estimated coef curves
par(mfrow=c(2,2))
plot(tt$coef.VC[,1]~time); points(alpha0~time, col=2, cex=0.5)
plot(tt$coef.VC[,2]~time); points(alpha1~time, col=2, cex=0.5)
plot(tt$coef.VC[,3]~time); points(alpha2~time, col=2, cex=0.5)
plot(tt$coef.VC[,4]~time); points(alpha3~time, col=2, cex=0.5)

# alpha2(t)
tau=0.5; test.VC="VC"; test.idx=2; degree=3; se="iid"; restricted="yes"
QRS.PLVC(y, time, design.VC, design.C, ID, tau=0.5, test.VC, test.idx, degree=3, se, restrcited="yes")


# test the constant effects
tau=0.5; test.VC="C"; test.idx=1; degree=3; se="iid"; restricted="yes"
QRS.PLVC(y, time, design.VC, design.C, ID, tau=0.5, test.VC, test.idx, degree=3, se, restrcited="yes")


tau=0.5; test.VC="C"; test.idx=2; degree=3; se="iid"; restricted="yes"
QRS.PLVC(y, time, design.VC, design.C, ID, tau=0.5, test.VC, test.idx, degree=3, se, restrcited="yes")