## modified by Marie-France Valois on August 7th, 2013
coxmfv.ns.ci.prg <- function(data, variable, fit, pos, df=NULL, knots=NULL, xref, x, typeout = "RR", level=0.95, na.action = T,link="log")
{
	#########################################################################################################
	#
	# The function cox.ns.ci.prg computes RR, logRR, or MPC and the 95% confidence limit
	# for a specified value x of the independent variable as compared to a specified reference value
	# xref .
	#
	# This function has 10 arguments which are described as follows:
	#
	# 	data:        Original S-Plus data set;
	# 	variable:    Name of the predictor variable;
	# 	fit:         The fitted object from the cox model;
	# 	pos: 	     the position where the first coefficient of the predictor variable begins in GLM summary;
	#  	df :         number of degrees of freedom (df) for the predictor variable. If no specified knots are
	#  	             given, the program will select the df + 1 knots (including boundary knots) at suitably
	#  	             chosen quantiles of the predictor variable;
	#  	knots:       knots that define the natural spline, including boundary points. When knots and
	#  	             df are both given, the number of knots is equal to df + 1;
	# 	xref :       the reference predictor value;
	#       x :          the value of the proposed variable to be compared to xref;
	#  	typeout:    You can choose "LOGRR", "RR", or "MPC" to output log Relative Risk, the
	#  	            Relative Risk, or the Mean Percent Change, respectively. The default is "RR"
	#       level:      The level for the confidence interval
        #  	na.action:  This parameter should be set to true (i.e., T) if the original regression model excluded 
	#                   observations that were flagged as missing by na.exclude. The default is T, which means to 
	#                   exclude missing values from the computations. If na.action=F, then na.fail will be applied,
        #                   which causes an error if missing values are present.
        #       link:       The link function used in the model
	############################################################################################################
	
        fitname <- substitute(fit)

        if(toupper(na.action) == F) {
		data = na.fail(data)
	}
        else { data = na.exclude(data)
        }
	# 	Get the variable column from data
	variable.vec <- data[, variable]
	# 	Get the range of variable
	variable.vec.range <- range(variable.vec, na.rm = TRUE)
	# 	if xref is out of range, a warning is given
	if((xref < variable.vec.range[1]) | (xref > variable.vec.range[2])) {
		warning(paste("the reference value xref is not in the range of data variabl"))
	}
	#  	if both df and knots are given, then the number of knots should be equal to df+1, otherwise
	# 	a warning will be given.
	##if(!missing(df) & !missing(knots)) {
	if(!is.null(df) & !is.null(knots)) {
		if((df + 1) != length(knots)) {
			warning(paste("df+1 is not equal to the number of knots"))
		}
	}
	# 	if knots is given instead of df, then df is equal to the number of knots - 1
	##if(missing(df) & (!missing(knots))) {
	if(is.null(df) & (!is.null(knots))) {
		df <- length(knots) - 1
	}
	# 	if df is given instead of knots, then the program will select the df + 1 knots (including boundary knots) at suitably
	#  	chosen quantiles of the predictor variable;
	##if(!missing(df) & (missing(knots))) {
	if(!is.null(df) & (is.null(knots))) {
		knots <- quantile(variable.vec, seq(0, 1, length = (df + 1)), na.rm = TRUE)
	}
	# 	if either df or knots is given, a warning is given.
	##if(missing(df) & missing(knots)) {
	if(is.null(df) & is.null(knots)) {
		warning(paste("please specify df or knots"))
	}
	# 	evaluate the basis functions at the reference point xref
	##basis0 <- ns(xref, knots = knots, Boundary.knots = TRUE)
	basis0 <- ns(xref, knots = knots[-c(1,df+1)], Boundary.knots = c(variable.vec.range[1],variable.vec.range[2]))
        # 	get the begining position of basis functions
	begin.pos <- pos
	# 	get the end position of basis functions
	end.pos <- pos + df - 1
	# 	get the estimated coefficients
	coeff <- fit$coefficients[begin.pos:end.pos]
	# 	get the covariance matrix of the estimated coefficients
	var <- fit$var[begin.pos:end.pos,begin.pos:end.pos]
	length.x <- length(x)
        for (i in 1:length.x)
        {
        # 	if x[i] is out of range, a warning is given
	if((x[i] < variable.vec.range[1]) | (x[i] > variable.vec.range[2])) {
		warning(paste("the specified value x is not in the range of data variabl"))
	}
	#	 evaluate the basis functions at the proposed point x[i]
	##basis1 <- ns(x[i], knots = knots, Boundary.knots = TRUE)
	basis1 <- ns(x[i], knots = knots[-c(1,df+1)], Boundary.knots = c(variable.vec.range[1],variable.vec.range[2]))
	# 	evaluate the difference of the bases functions at the 
	# 	proposed point x[i] from the reference point xref 
	diff <- basis1 - basis0
		# evaluate the variance of the log relative risk
	##varlogrr <- t(diff) %*% var %*% diff
	varlogrr <- t(as.vector(diff)) %*% var %*% as.vector(diff)
	# evaluate the log relative risk
	##logrr <- t(diff) %*% coeff
	logrr <- t(as.vector(diff)) %*% coeff
	# get the confidence interval for the log relative risk
        normstat <- qnorm(1-(1-level)/2)
	lower.logrr <- logrr - normstat * sqrt(varlogrr)
	upper.logrr <- logrr + normstat * sqrt(varlogrr)
        if (toupper(link)=="LOG"){
	# output mean percent change and the associated confidence interval if typeout is not "RR" or "LOGRR"
	# output log relative risk and the associated confidence interval if typeout = "LOGRR"
	# output relative risk and the associated confidence interval if typeout = "RR"
	if(toupper(typeout) == "RR") {
		rr <- exp(logrr)
		lower.rr <- exp(lower.logrr)
		upper.rr <- exp(upper.logrr)
		cat("BOR", "\n")
		cat("Fitted Model is : ", as.name(fitname), "\n")
		cat("Effect is : ", variable, "\n")
                cat("When the value is  ", x[i], "\n")
		cat("RR  =  ", rr, "\n")
		cat("Lower CI of RR = ", lower.rr, "\n")
		cat("Upper CI of RR = ", upper.rr, "\n")
		cat("EOR", "\n\n\n")
	}
	else if(toupper(typeout) == "LOGRR") {
		cat("BOR", "\n")
		cat("Fitted Model is : ", as.name(fitname), "\n")
		cat("Effect is : ", variable, "\n")
                cat("When the value is  ", x[i], "\n")
		cat("LOG RR  =  ", logrr, "\n")
		cat("Lower CI of LOG RR = ", lower.logrr, "\n")
		cat("Upper CI of LOG RR = ", upper.logrr, "\n")
		cat("EOR", "\n\n\n")
	}
	else {
		mpc <- (exp(logrr) - 1) * 100
		lower.mpc <- (exp(lower.logrr) - 1) * 100
		upper.mpc <- (exp(upper.logrr) - 1) * 100
		cat("BOR", "\n")
		cat("Fitted Model is : ", as.name(fitname), "\n")
		cat("Effect is : ", variable, "\n")
                cat("When the value is  ", x[i], "\n")
		cat("MPC  =  ", mpc, "%\n")
		cat("Lower CI = ", lower.mpc, "\n")
		cat("Upper CI = ", upper.mpc, "\n")
		cat("EOR", "\n\n\n")
	}
        }
        else {
                cat("BOR", "\n")
		cat("Fitted Model is : ", as.name(fitname), "\n")
		cat("Effect is : ", variable, "\n")
                cat("When the value is  ", x[i], "\n")
		cat("RR  =  ", logrr, "\n")
		cat("Lower CI of RR = ", lower.logrr, "\n")
		cat("Upper CI of RR = ", upper.logrr, "\n")
		cat("EOR", "\n\n\n")
       }
       }
}


## modified by Kashima 31 Jan. 2011.
ns.ci.prg <- function(data_input, variable, fit, pos, df=NULL, knots=NULL, xref, x, typeout = "RR", level=0.95, na.action = na.exclude,linkf="log")
{
	#########################################################################################################
	# modified by Saori.Kashima on 31 Jan 2011
	#
	# The function ns.ci.prg computes RR, logRR, or MPC and the 95% confidence limit
	# for a specified value x of the independent variable as compared to a specified reference value
	# xref .
	#
	# This function has 10 arguments which are described as follows:
	#
	# 	data:        Original S-Plus data set;
	# 	variable:    Name of the predictor variable;
	# 	fit:         The fitted object from the generalized linear model;
	# 	pos: 	     the position where the first coefficient of the predictor variable begins in GLM summary;
	#  	df :         number of degrees of freedom (df) for the predictor variable. If no specified knots are
	#  	             given, the program will select the df + 1 knots (including boundary knots) at suitably
	#  	             chosen quantiles of the predictor variable;
	#  	knots:       knots that define the natural spline, including boundary points. When knots and
	#  	             df are both given, the number of knots is equal to df + 1;
	# 	xref :       the reference predictor value;
	#       x :          the value of the proposed variable to be compared to xref;
	#  	typeout:    You can choose "LOGRR", "RR", or "MPC" to output log Relative Risk, the
	#  	            Relative Risk, or the Mean Percent Change, respectively. The default is "RR"
	#       level:      The level for the confidence interval
	#  	na.action:  This parameter should be set to na.exclude if the original regression model excluded 
	#                   observations that were flagged as missing by na.exclude. The default is na.exclude, which means to 
	#                   exclude missing values from the computations. If na.action=na.fail, then na.fail will be applied,
	#                   which causes an error if missing values are present.
	#       linkf:       The link function used in the model
	############################################################################################################
	## modified by Kashima 31 Jan. 2011.
	########################################################################################################

	## modified by Kashima
	##data = na.action(data)
	data = data_input
	
	# 	Get the variable column from data
	variable.vec <- data[, variable]
	# 	Get the range of variable
	variable.vec.range <- range(variable.vec, na.rm = TRUE)
	# 	if xref is out of range, a warning is given
	if((xref < variable.vec.range[1]) | (xref > variable.vec.range[2])) {
		warning(paste("the reference value xref is not in the range of data variabl"))
	}
	#  	if both df and knots are given, then the number of knots should be equal to df+1, otherwise
	# 	a warning will be given.
	##if(!missing(df) & !missing(knots)) {
	if(!is.null(df) & !is.null(knots)) {
		if((df + 1) != length(knots)) {
			warning(paste("df+1 is not equal to the number of knots"))
		}
	}
	# 	if knots is given instead of df, then df is equal to the number of knots - 1
	##if(missing(df) & (!missing(knots))) {
	if(is.null(df) & (!is.null(knots))) {
		df <- length(knots) - 1
	}
	# 	if df is given instead of knots, then the program will select the df + 1 knots (including boundary knots) at suitably
	#  	chosen quantiles of the predictor variable;
	## modified by Kashima
	##if(!missing(df) & (missing(knots))) {
	if(!is.null(df) & (is.null(knots))) {
		knots <- quantile(variable.vec, seq(0, 1, length = (df + 1)), na.rm = TRUE)
	}
	# 	if either df or knots is given, a warning is given.
	##if(missing(df) & missing(knots)) {
	if(is.null(df) & is.null(knots)) {
		warning(paste("please specify df or knots"))
	}
	# 	evaluate the basis functions at the reference point xref
	## modified by Kashima 31 Jan. 2011.
	##basis0 <- ns(xref, knots = knots, Boundary.knots = TRUE)
	basis0 <- ns(xref, knots = knots[-c(1,df+1)], Boundary.knots = c(variable.vec.range[1],variable.vec.range[2]))
	# 	get the begining position of basis functions
	begin.pos <- pos
	# 	get the end position of basis functions
	end.pos <- pos + df - 1
	# 	get the estimated coefficients
	coeff <- summary(fit)$coefficients[begin.pos:end.pos, 1]
	#	 get the standard deviations of the estimated coefficients
	std <- summary(fit)$coefficients[begin.pos:end.pos, 2]
	# 	get the correlation matrix of the estimated coefficients      
	##  modified by Kashima.
	##corr <- summary(fit)$correlation[begin.pos:end.pos, begin.pos:end.pos]
	corr <- summary(fit, correlation=TRUE)$correlation[begin.pos:end.pos, begin.pos:end.pos]
	# 	compute the covariance matrix of the estimated coefficients
	var <- matrix(0, df, df)
	for(i in 1:df) {
		for(j in i:df) {
			var[i, j] <- std[i] * std[j] * corr[i, j]
		}
	}
	for(i in 1:df) {
		for(j in 1:(i - 1)) {
			var[i, j] <- var[j, i]
		}
	}
        length.x <- length(x)
        for (i in 1:length.x)
	{
	# 	if x[i] is out of range, a warning is given
	if((x[i] < variable.vec.range[1]) | (x[i] > variable.vec.range[2])) {
		warning(paste("the specified value x is not in the range of data variabl"))
	}
	#	 evaluate the basis functions at the proposed point x[i]
	## modified by Kashima 31 Jan. 2011.
	##basis1 <- ns(x[i], knots = knots, Boundary.knots = TRUE)
	basis1 <- ns(x[i], knots = knots[-c(1,df+1)], Boundary.knots = c(variable.vec.range[1],variable.vec.range[2]))
	# 	evaluate the difference of the bases functions at the 
	# 	proposed point x[i] from the reference point xref 
	diff <- basis1 - basis0
	# evaluate the variance of the log relative risk
	## modified by Kashima.
	## varlogrr <- t(diff) %*% var %*% diff
	varlogrr <- t(as.vector(diff)) %*% var %*% as.vector(diff)
	# evaluate the log relative risk
	## modified by Kashima.
	##logrr <- t(diff) %*% coeff
	logrr <- t(as.vector(diff)) %*% coeff
	# get the confidence interval for the log relative risk
	normstat <- qnorm(1-(1-level)/2)
	lower.logrr <- logrr - normstat * sqrt(varlogrr)
	upper.logrr <- logrr + normstat * sqrt(varlogrr)
	if (toupper(linkf)=="LOG"){
	# output mean percent change and the associated confidence interval if typeout is not "RR" or "LOGRR"
	# output log relative risk and the associated confidence interval if typeout = "LOGRR"
	# output relative risk and the associated confidence interval if typeout = "RR"
	if(toupper(typeout) == "RR") {
		rr <- exp(logrr)
		lower.rr <- exp(lower.logrr)
		upper.rr <- exp(upper.logrr)
		cat("BOR", "\n")
		cat("Effect is : ", variable, "\n")
                cat("When the value is  ", x[i], "\n")
		cat("RR  =  ", rr, "\n")
		cat("Lower CI of RR = ", lower.rr, "\n")
		cat("Upper CI of RR = ", upper.rr, "\n")
		cat("EOR", "\n\n\n")
	}
	else if(toupper(typeout) == "LOGRR") {
		cat("BOR", "\n")
		cat("Effect is : ", variable, "\n")
		cat("When the value is  ", x[i], "\n")
		cat("LOG RR  =  ", logrr, "\n")
		cat("Lower CI of LOG RR = ", lower.logrr, "\n")
		cat("Upper CI of LOG RR = ", upper.logrr, "\n")
		cat("EOR", "\n\n\n")
	}
	else {
		mpc <- (exp(logrr) - 1) * 100
		lower.mpc <- (exp(lower.logrr) - 1) * 100
		upper.mpc <- (exp(upper.logrr) - 1) * 100
		cat("BOR", "\n")
		cat("Effect is : ", variable, "\n")
		cat("When the value is  ", x[i], "\n")
		cat("MPC  =  ", mpc, "%\n")
		cat("Lower CI = ", lower.mpc, "%\n")
		cat("Upper CI = ", upper.mpc, "%\n")
		cat("EOR", "\n\n\n")
	}
	}
	else {
                cat("BOR", "\n")
		cat("Effect is : ", variable, "\n")
		cat("When the value is  ", x[i], "\n")
		cat("MD  =  ", logrr, "\n")
		cat("Lower CI of MD = ", lower.logrr, "\n")
		cat("Upper CI of MD = ", upper.logrr, "\n")
		cat("EOR", "\n\n\n")
	}
	}
}


NSCILME <- function(data_input, variable, fit, pos, df=NULL, knots=NULL, xref, x, typeout = "MC", level=0.95, na.action = na.exclude)
{
	#########################################################################################################
	#
	# The function ns.ci.lme.prg computes the MC (Mean Change) and the 95% confidence limit
	# for a specified value x of the independent variable as compared to a specified reference value
	# xref .
	#
	# This function has 10 arguments which are described as follows:
	#
	# 	data:        Original S-Plus data set;
	# 	variable:    Name of the predictor variable;
	# 	fit:         The fitted object from the mixed model;
	# 	pos: 	     the position where the first coefficient of the predictor variable begins in Mixed model summary;
	#  	df :         number of degrees of freedom (df) for the predictor variable. If no specified knots are
	#  	             given, the program will select the df + 1 knots (including boundary knots) at suitably
	#  	             chosen quantiles of the predictor variable;
	#  	knots:       knots that define the natural spline, including boundary points. When knots and
	#  	             df are both given, the number of knots is equal to df + 1;
	# 	xref :       the reference predictor value;
	#       x :          the value of the proposed variable to be compared to xref;
	#  	typeout:    You can only choose "MC" to output the Mean Change. The default is "MC"
	#       level:      The level for the confidence interval
        #  	na.action:  This parameter should be set to na.exclude if the original regression model excluded 
	#                   observations that were flagged as missing by na.exclude. The default is na.exclude, which means to 
	#                   exclude missing values from the computations. If na.action=na.fail, then na.fail will be applied,
        #                   which causes an error if missing values are present.
	############################################################################################################
	
        ##data = na.action(data)
	data = data_input
	# 	Get the variable column from data
	variable.vec <- data[, variable]
	# 	Get the range of variable
	variable.vec.range <- range(variable.vec, na.rm = TRUE)
	# 	if xref is out of range, a warning is given
	if((xref < variable.vec.range[1]) | (xref > variable.vec.range[2])) {
		warning(paste("the reference value xref is not in the range of data variable"))
	}
	#  	if both df and knots are given, then the number of knots should be equal to df+1, otherwise
	# 	a warning will be given.
	##if(!missing(df) & !missing(knots)) {
	if(!is.null(df) & !is.null(knots)) {
		if((df + 1) != length(knots)) {
			warning(paste("df+1 is not equal to the number of knots"))
		}
	}
	# 	if knots is given instead of df, then df is equal to the number of knots - 1
	##if(missing(df) & (!missing(knots))) {
	if(is.null(df) & (!is.null(knots))) {
		df <- length(knots) - 1
	}
	# 	if df is given instead of knots, then the program will select the df + 1 knots (including boundary knots) at suitably
	#  	chosen quantiles of the predictor variable;
	##if(!missing(df) & (missing(knots))) {
	if(!is.null(df) & (is.null(knots))) {
		knots <- quantile(variable.vec, seq(0, 1, length = (df + 1)), na.rm = TRUE)
	}
	# 	if either df or knots is given, a warning is given.
	##if(missing(df) & missing(knots)) {
	if(is.null(df) & is.null(knots)) {
		warning(paste("please specify df or knots"))
	}
	# 	evaluate the basis functions at the reference point xref
	##basis0 <- ns(xref, knots = knots, Boundary.knots = TRUE)
	basis0 <- ns(xref, knots = knots[-c(1,df+1)], Boundary.knots = c(variable.vec.range[1],variable.vec.range[2]))
        # 	get the begining position of basis functions
	begin.pos <- pos
	# 	get the end position of basis functions
	end.pos <- pos + df - 1
	# 	get the estimated coefficients
	coeff <- summary(fit)$tTable[begin.pos:end.pos, 1]
	#	 get the standard deviations of the estimated coefficients
	std <- summary(fit)$tTable[begin.pos:end.pos, 2]
	# 	get the correlation matrix of the estimated coefficients      
	##corr <- summary(fit)$corFixed[begin.pos:end.pos, begin.pos:end.pos]
	corr <- summary(fit)$corFixed[begin.pos:end.pos, begin.pos:end.pos]
	# 	compute the covariance matrix of the estimated coefficients
	var <- matrix(0, df, df)
	for(i in 1:df) {
		for(j in i:df) {
			var[i, j] <- std[i] * std[j] * corr[i, j]
		}
	}
	for(i in 1:df) {
		for(j in 1:(i - 1)) {
			var[i, j] <- var[j, i]
		}
	}

        length.x <- length(x)
        for (i in 1:length.x)
        {
        # 	if x[i] is out of range, a warning is given
	if((x[i] < variable.vec.range[1]) | (x[i] > variable.vec.range[2])) {
		warning(paste("the specified value x is not in the range of data variable"))
	}
	#	 evaluate the basis functions at the proposed point x[i]
	##basis1 <- ns(x[i], knots = knots, Boundary.knots = TRUE)
	basis1 <- ns(x[i], knots = knots[-c(1,df+1)], Boundary.knots = c(variable.vec.range[1],variable.vec.range[2]))
	# 	evaluate the difference of the bases functions at the 
	# 	proposed point x[i] from the reference point xref 
	diff <- basis1 - basis0
		# evaluate the variance of the log relative risk
	##varmc <- t(diff) %*% var %*% diff
	varmc <- t(as.vector(diff)) %*% var %*% as.vector(diff)

	# evaluate the log relative risk
	##mc <- t(diff) %*% coeff
	mc <- t(as.vector(diff)) %*% coeff
	# get the confidence interval for the log relative risk
        normstat <- qnorm(1-(1-level)/2)
	lower.mc <- mc - normstat * sqrt(varmc)
	upper.mc <- mc + normstat * sqrt(varmc)

	# output mean change and the associated confidence interval
	cat("BOR", "\n")
	cat("Effect is : ", variable, "\n")
        cat("When the value is  ", x[i], "\n")
	cat("MC  =  ", mc, "\n")
	cat("Lower CI = ", lower.mc, "\n")
	cat("Upper CI = ", upper.mc, "\n")
	cat("EOR", "\n\n\n")
}
}


## Modified by Marie-France Valois on August 7th, 2013
plot.cox.ns.ci.prg <- function(data, variable, fit, pos, df=NULL, knots=NULL, xref, typeout = "RR", step = 0.1, xtitle = toupper(variable),ytitle = toupper(typeout), level=0.95, na.action = T, col=1,lty="l",lwd=1,link="log")
{
	#####################################################################################################
	#	
	# The function plot.cox.ns.ci.prg produces a publication-quality graph of RR, logRR, or
	# MPC and the associated 95% confidence limits for an independent variable with respect to a
	# specified reference value xref . 
	#
	# This function has 12 arguments which are described as follows:
	#
	# 	data:       Original S-Plus data set;
	# 	variable:   Name of the predictor variable;
	# 	fit:        The fitted object from the generalized linear model;
	# 	pos:        the position where the first coefficient of the predictor variable begins in GLM summary;
	#  	df :        number of degrees of freedom (df) for the predictor variable. If no specified knots are
	#  	            given, the program will select the df + 1 knots (including boundary knots) at suitably
	#  	            chosen quantiles of the predictor variable;
	#  	knots:      knots that define the natural spline, including boundary points. When knots and
	#  	            df are both given, the number of knots is equal to df + 1;
	# 	xref :      the reference predictor value;
        #     level:      The level for the confidence interval
	#  	typeout:    You can choose "LOGRR", "RR", or "MPC" to output log Relative Risk, the
	#  	            Relative Risk, or the Mean Percent Change, respectively. The default is "RR"
	#  	step:       the step length to generate a vector of predictor values to get their RR, or logRR, 
	#                   or MPC as compared to xref . The smaller the step, the more continuous the graph
	#  	            looks, but the longer it will take to run the program. The default is 0.1;
	#  	xtitle:     the title of the x-axis. The default is the predictor variable (in uppercase);
	# 	ytitle:     the title of the y-axis. The default is the argument typeout (in uppercase);
	#  	na.action:  This parameter should be set to true (i.e., T) if the original regression model excluded 
	#                   observations that were flagged as missing by na.exclude. The default is T, which means to 
	#                   exclude missing values from the computations. If na.action=F, then na.fail will be applied,
        #                   which causes an error if missing values are present.
        #       col:        The color of lines
        #       lty:        The types of lines
        #       lwd:        The width of lines
	#
	########################################################################################################
     
	if(toupper(na.action) == F) {
		data = na.fail(data)
	}
        else { data = na.exclude(data)
        }
        
	# 	Get the variable column from data
	variable.vec <- data[, variable]
	# 	Get the range of variable
	variable.vec.range <- range(variable.vec, na.rm = TRUE)
	# 	if xref is out of range, a warning is given	
	if((xref < variable.vec.range[1]) | (xref > variable.vec.range[2])) {
		warning(paste("the reference value xref is not in the range of data variable"))
	}
	# Generate a sequence of points for which the Log RR as compared to the reference point xref will be calculated.
	xvec = seq(variable.vec.range[1], variable.vec.range[2], by = step)
	#  	if both df and knots are given, then the number of knots should be equal to df+1, otherwise
	# 	a warning will be given.
	##if(!missing(df) & !missing(knots)) {
	if(!is.null(df) & !is.null(knots)) {
		if((df + 1) != length(knots)) {
			warning(paste("df+1 is not equal to the number of knots"))
		}
	}
	# 	if knots is given instead of df, then df is equal to the number of knots - 1
	##if(missing(df) & (!missing(knots))) {
	if(is.null(df) & (!is.null(knots))) {
		df <- length(knots) - 1
	}
	# 	if df is given instead of knots, then the program will select the df + 1 knots (including boundary knots) at suitably
	#  	chosen quantiles of the predictor variable;
	##if(!missing(df) & (missing(knots))) {
	if(!is.null(df) & (is.null(knots))) {
		knots <- quantile(variable.vec, seq(0, 1, length = (df + 1)), na.rm = TRUE)
	}
	# 	if either df or knots is given, a warning is given.
	##if(missing(df) & missing(knots)) {
	if(is.null(df) & is.null(knots)) {
		warning(paste("please specify df or knots"))
	}
	# 	get the begining position of basis functions
	begin.pos <- pos
	# 	get the end position of basis functions
	end.pos <- pos + df - 1
	# 	get the estimated coefficients
        coeff <- fit$coefficients[begin.pos:end.pos]
	# 	get the covariance matrix of the estimated coefficients
	var <- fit$var[begin.pos:end.pos,begin.pos:end.pos]
	# evaluate the basis functions at the reference point xref
	##basis0 <- ns(xref, knots = knots, Boundary.knots = TRUE)
	basis0 <- ns(xref, knots = knots[-c(1,df+1)], Boundary.knots = c(variable.vec.range[1],variable.vec.range[2]))
	# extract the length of the vector needed to save logRR.
	nx <- length(xvec)
	# Initialize the vectors of the log RR and lower and upper bound of confidence interval of log RR
	logrr <- rep(0, nx)
	lower.logrr <- rep(0, nx)
	upper.logrr <- rep(0, nx)
	# calculate the log RR and confidence interval for each point
	for(i in 1:nx) {
		# evaluate the basis functions at the proposed point x	
		##basis1 <- ns(xvec[i], knots = knots, Boundary.knots = TRUE)
		basis1 <- ns(xvec[i], knots = knots[-c(1,df+1)], Boundary.knots = c(variable.vec.range[1],variable.vec.range[2]))
		# evaluate the difference of the bases functions at the 
		# proposed point x from the reference point xref 
		diff <- basis1 - basis0
		# evaluate the variance of log rate ratio
		##varlogrr <- t(diff) %*% var %*% diff
		varlogrr <- t(as.vector(diff)) %*% var %*% as.vector(diff)
		# evaluate the log rate ratio
		##logrr[i] <- t(diff) %*% coeff
		logrr[i] <- t(as.vector(diff)) %*% coeff
		# calculate the lower and upper bound of confidence interval of log RR
                normstat <- qnorm(1-(1-level)/2)
		lower.logrr[i] <- logrr[i] - normstat* sqrt(varlogrr)
		upper.logrr[i] <- logrr[i] + normstat* sqrt(varlogrr)
	}
        if (toupper(link)=="LOG"){
	# output the graph of mean percent change and the associate confidence interval if typeout is not "RR" or "LOGRR"
	# output the graph of log relative risk and the associate confidence interval if typeout = "LOGRR"
	# output the graph of relative risk and the associate confidence interval if typeout = "RR"
	if(toupper(typeout) == "RR") {
		plot(xvec, exp(logrr), ylim = c(exp(min(lower.logrr)), exp(max(upper.logrr))), type = lty, xlab = xtitle, ylab = ytitle, col=col,lwd=lwd)
		lines(xvec, exp(lower.logrr), type = lty,col=col,lwd=lwd)
		lines(xvec, exp(upper.logrr), type = lty,col=col,lwd=lwd)
		abline(h = 1)
		abline(v = xref)
	}
	else if(toupper(typeout) == "LOGRR") {
		plot(xvec, logrr, ylim = c(min(lower.logrr), max(upper.logrr)), type = lty, xlab = xtitle, ylab = ytitle, col=col,lwd=lwd)
		lines(xvec, lower.logrr, type = lty,col=col,lwd=lwd)
		lines(xvec, upper.logrr, type = lty,col=col,lwd=lwd)
		abline(h = 0)
		abline(v = xref)
	}
	else {
		plot(xvec, (exp(logrr) - 1) * 100, ylim = c((exp(min(lower.logrr)) - 1) * 100, (exp(max(upper.logrr)) - 1) * 100), type = lty, xlab = xtitle, 
                    ylab = ytitle, , col=col,lwd=lwd)
		lines(xvec, (exp(lower.logrr) - 1) * 100, type = lty,col=col,lwd=lwd)
		lines(xvec, (exp(upper.logrr) - 1) * 100, type = lty,col=col,lwd=lwd)
		abline(h = 0)
		abline(v = xref)
	}
        }
        else {
        plot(xvec, logrr, ylim = c(min(lower.logrr), max(upper.logrr)), type = lty, xlab = xtitle, ylab = ytitle, col=col,lwd=lwd)
		lines(xvec, lower.logrr, type = lty,col=col,lwd=lwd)
		lines(xvec, upper.logrr, type = lty,col=col,lwd=lwd)
		abline(h = 0)
		abline(v = xref)
        }
}



########################
## modified by Kashima 31 Jan. 2011.
plot.ns.ci.prg <- function(data_input, variable, fit, pos=NULL, df, knots=NULL, xref, typeout = "RR", step = 0.1, xtitle = toupper(variable),ytitle = toupper(typeout), level=0.95, na.action = T, col=1,lty="l",lwd=1,linkf="log", ylim=NULL, xlim=NULL)
{
	#####################################################################################################
	#
	# The function plot.ns.ci.prg produces a publication-quality graph of RR, logRR, or
	# MPC and the associated 95% confidence limits for an independent variable with respect to a
	# specified reference value xref . 
	#
	# This function has 12 arguments which are described as follows:
	#
	# 	data:       Original S-Plus data set;
	# 	variable:   Name of the predictor variable;
	# 	fit:        The fitted object from the generalized linear model;
	# 	pos:        the position where the first coefficient of the predictor variable begins in GLM summary;
	#  	df :        number of degrees of freedom (df) for the predictor variable. If no specified knots are
	#  	            given, the program will select the df + 1 knots (including boundary knots) at suitably
	#  	            chosen quantiles of the predictor variable;
	#  	knots:      knots that define the natural spline, including boundary points. When knots and
	#  	            df are both given, the number of knots is equal to df + 1;
	# 	xref :      the reference predictor value;
	#  	typeout:    You can choose "LOGRR", "RR", or "MPC" to output log Relative Risk, the
	#  	            Relative Risk, or the Mean Percent Change, respectively. The default is "RR"
	#  	step:       the step length to generate a vector of predictor values to get their RR, or logRR, 
	#                   or MPC as compared to xref . The smaller the step, the more continuous the graph
	#  	            looks, but the longer it will take to run the program. The default is 0.1;
	#  	xtitle:     the title of the x-axis. The default is the predictor variable (in uppercase);
	# 	ytitle:     the title of the y-axis. The default is the argument typeout (in uppercase);
	#       level:      The level for the confidence interval
	#  	na.action:  This parameter should be set to na.exclude if the original regression model excluded 
	#                   observations that were flagged as missing by na.exclude. The default is na.exclude, which means to 
	#                   exclude missing values from the computations. If na.action=na.fail, then na.fail will be applied,
	#                   which causes an error if missing values are present.
	#       col:        The color of lines
	#       lty:        The types of lines
	#       lwd:        The width of lines
	#       linkf:       The link function used in the model
	########################################################################################################
	## modified by Kashima 31 Jan. 2011.
	########################################################################################################
	
	## modified by Kashima.
	##data = na.action(data)
	data = data_input
	# 	Get the variable column from data
	variable.vec <- data[, variable]
	# 	Get the range of variable
	variable.vec.range <- range(variable.vec, na.rm = TRUE)
	# 	if xref is out of range, a warning is given	
	if((xref < variable.vec.range[1]) | (xref > variable.vec.range[2])) {
		warning(paste("the reference value xref is not in the range of data variable"))
	}
	# Generate a sequence of points for which the Log RR as compared to the reference point xref will be calculated.
	xvec = seq(variable.vec.range[1], variable.vec.range[2], by = step)
	#  	if both df and knots are given, then the number of knots should be equal to df+1, otherwise
	# 	a warning will be given.
	##if(!missing(df) & !missing(knots)) {
	if(!is.null(df) & !is.null(knots)) {
		if((df + 1) != length(knots)) {
			warning(paste("df+1 is not equal to the number of knots"))
		}
	}
	# 	if knots is given instead of df, then df is equal to the number of knots - 1
	##if(missing(df) & (!missing(knots))) {
	if(is.null(df) & (!is.null(knots))) {
		df <- length(knots) - 1
	}
	# 	if df is given instead of knots, then the program will select the df + 1 knots (including boundary knots) at suitably
	#  	chosen quantiles of the predictor variable;
	##if(!missing(df) & (missing(knots))) {
	if(!is.null(df) & (is.null(knots))) {
		knots <- quantile(variable.vec, seq(0, 1, length = (df + 1)), na.rm = TRUE)
	}
	# 	if either df or knots is given, a warning is given.
	##if(missing(df) & missing(knots)) {
	if(is.null(df) & is.null(knots)) {
		warning(paste("please specify df or knots"))
	}
	# 	get the begining position of basis functions
	begin.pos <- pos
	# 	get the end position of basis functions
	end.pos <- pos + df - 1
	# 	get the estimated coefficients
	coeff <- summary(fit)$coefficients[begin.pos:end.pos, 1]
	# 	get the correlation matrix of the estimated coefficients      
	## modified by Kashima.
	##corr <- summary(fit)$correlation[begin.pos:end.pos, begin.pos:end.pos]
	corr <- summary(fit, correlation=TRUE)$correlation[begin.pos:end.pos, begin.pos:end.pos]
	#	 get the standard deviations of the estimated coefficients
	std <- summary(fit)$coefficients[begin.pos:end.pos, 2]
	# 	 compute the  covariance matrix of estimated coefficients of basis 
	#        natural splines.
	var <- matrix(0, df, df)
	for(i in 1:df) {
		for(j in i:df) {
			var[i, j] <- std[i] * std[j] * corr[i, j]
		}
	}
	for(i in 1:df) {
		for(j in 1:(i - 1)) {
			var[i, j] <- var[j, i]
		}
	}
	# evaluate the basis functions at the reference point xref
	## modified by Kashima.
	##basis0 <- ns(xref, knots = knots, Boundary.knots = TRUE)
	basis0 <- ns(xref, knots = knots[-c(1,df+1)], Boundary.knots = c(variable.vec.range[1],variable.vec.range[2]))
	# extract the length of the vector needed to save logRR.
	nx <- length(xvec)
	# Initialize the vectors of the log RR and lower and upper bound of confidence interval of log RR
	logrr <- rep(0, nx)
	lower.logrr <- rep(0, nx)
	upper.logrr <- rep(0, nx)
	# calculate the log RR and confidence interval for each point
	for(i in 1:nx) {
		# evaluate the basis functions at the proposed point x	
		## modified by Kashima 31 Jan. 2011.
		##basis1 <- ns(xvec[i], knots = knots, Boundary.knots = TRUE)
		basis1 <- ns(xvec[i], knots = knots[-c(1,df+1)], Boundary.knots = c(variable.vec.range[1],variable.vec.range[2]))
		# evaluate the difference of the bases functions at the 
		# proposed point x from the reference point xref 
		diff <- basis1 - basis0
		# evaluate the variance of log rate ratio
		## modified by Kashima.
		##varlogrr <- t(diff) %*% var %*% diff
		varlogrr <- t(as.vector(diff)) %*% var %*% as.vector(diff)
		# evaluate the log rate ratio
		## modified by Kashima.
		##logrr[i] <- t(diff) %*% coeff
		logrr[i] <- t(as.vector(diff)) %*% coeff
		# calculate the lower and upper bound of confidence interval of log RR
 		normstat <- qnorm(1-(1-level)/2)
		lower.logrr[i] <- logrr[i] - normstat * sqrt(varlogrr)
		upper.logrr[i] <- logrr[i] + normstat * sqrt(varlogrr)
	}
	if (toupper(linkf)=="LOG"){
	# output the graph of mean percent change and the associate confidence interval if typeout is not "RR" or "LOGRR"
	# output the graph of log relative risk and the associate confidence interval if typeout = "LOGRR"
	# output the graph of relative risk and the associate confidence interval if typeout = "RR"
	if(toupper(typeout) == "RR") {
		## modified by Kashima.
		if (is.null(ylim)){
			ylim=c(exp(min(lower.logrr)), exp(max(upper.logrr)))
		}
                ##plot(xvec, exp(logrr), ylim = c(exp(min(lower.logrr)), exp(max(upper.logrr))), type = lty, xlab = xtitle, ylab = ytitle, col=col,lwd=lwd)
                plot(xvec, exp(logrr), ylim = ylim, xlim = xlim, type = lty, xlab = xtitle, ylab = ytitle, col=col,lwd=lwd)
		##
		abline(h = 1, col="gray")
		abline(v = xref, col="gray")
		par(new=T)
		## modified by Kashima.
		##plot(xvec, exp(logrr), ylim = c(exp(min(lower.logrr)), exp(max(upper.logrr))), type = lty, xlab = xtitle, ylab = ytitle, col=col,lwd=lwd)
		plot(xvec, exp(logrr), ylim = ylim, xlim = xlim, type = lty, xlab = xtitle, ylab = ytitle, col=col,lwd=lwd)
		lines(xvec, exp(lower.logrr), type = lty,lwd=lwd, lty=2)
		lines(xvec, exp(upper.logrr), type = lty,lwd=lwd, lty=2)

	}
	else if(toupper(typeout) == "LOGRR") {
		## modified by Kashima.
		if (is.null(ylim)){
			ylim=c(min(lower.logrr), max(upper.logrr))
		}
		## modified by Kashima.
		##plot(xvec, logrr, ylim = c(min(lower.logrr), max(upper.logrr)), type = lty, xlab = xtitle, ylab = ytitle, col=col,lwd=lwd)
                plot(xvec, logrr, ylim = ylim, xlim = xlim, type = lty, xlab = xtitle, ylab = ytitle, col=col,lwd=lwd)
		abline(h = 0, col="gray")
		abline(v = xref, col="gray")
		par(new=T)
		## modified by Kashima.
		##plot(xvec, logrr, ylim = c(min(lower.logrr), max(upper.logrr)), type = lty, xlab = xtitle, ylab = ytitle, col=col,lwd=lwd)
		plot(xvec, logrr, ylim = ylim, xlim = xlim, type = lty, xlab = xtitle, ylab = ytitle, col=col,lwd=lwd)
		lines(xvec, lower.logrr, type = lty,lwd=lwd, lty=2)
		lines(xvec, upper.logrr, type = lty,lwd=lwd, lty=2)
	}
	else {
		## modified by Kashima.
		if (is.null(ylim)){
			ylim=c((exp(min(lower.logrr)) - 1) * 100, (exp(max(upper.logrr)) - 1) * 100)
		}
		## modified by Kashima.
		##plot(xvec, (exp(logrr) - 1) * 100, ylim = c((exp(min(lower.logrr)) - 1) * 100, (exp(max(upper.logrr)) - 1) * 100), type = lty, xlab = xtitle, 
                ##    ylab = ytitle, , col=col,lwd=lwd)
		plot(xvec, (exp(logrr) - 1) * 100, ylim = ylim, xlim = xlim, type = lty, xlab = xtitle, 
                    ylab = ytitle, , col=col,lwd=lwd)
		abline(h = 0, col="gray")
		abline(v = xref, col="gray")
		par(new=T)
		## modified by Kashima.
		##plot(xvec, (exp(logrr) - 1) * 100, ylim = c((exp(min(lower.logrr)) - 1) * 100, (exp(max(upper.logrr)) - 1) * 100), type = lty, xlab = xtitle, 
		##    ylab = ytitle, , col=col,lwd=lwd)
		plot(xvec, (exp(logrr) - 1) * 100, ylim = ylim, xlim = xlim, type = lty, xlab = xtitle, 
                    ylab = ytitle, , col=col,lwd=lwd)
		lines(xvec, (exp(lower.logrr) - 1) * 100, type = lty,lwd=lwd, lty=2)
		lines(xvec, (exp(upper.logrr) - 1) * 100, type = lty,lwd=lwd, lty=2)
	}
        }
        else {
		## modified by Kashima.
		if (is.null(ylim)){
			ylim=c(min(lower.logrr), max(upper.logrr))
		}
		## modified by Kashima.
		##plot(xvec, logrr, ylim = c(min(lower.logrr), max(upper.logrr)), type = lty, xlab = xtitle, ylab = ytitle, col=col,lwd=lwd)
		plot(xvec, logrr, ylim = ylim, xlim=xlim, type = lty, xlab = xtitle, ylab = ytitle, col=col,lwd=lwd)
		abline(h = 0, col="gray")
		abline(v = xref, col="gray")
		par(new=T)
		## modified by Kashima.
		##plot(xvec, logrr, ylim = c(min(lower.logrr), max(upper.logrr)), type = lty, xlab = xtitle, ylab = ytitle, col=col,lwd=lwd)
        	plot(xvec, logrr, ylim = ylim, xlim=xlim, type = lty, xlab = xtitle, ylab = ytitle, col=col,lwd=lwd)
		lines(xvec, lower.logrr, type = lty,lwd=lwd, lty=2)
		lines(xvec, upper.logrr, type = lty,lwd=lwd, lty=2)
        }
}



PLOTNSCILME <- function(data_input, variable, fit, pos, df=NULL, knots=NULL, xref, typeout = "MC", step = 0.1, xtitle = upperCase(variable),ytitle = upperCase(typeout), level=0.95, na.action = na.exclude, col=1,lty="l",lwd=1)
{
	#####################################################################################################
	#	
	# The function plot.ns.ci.lme.prg produces a publication-quality graph of MC (Mean Change) and 
        # the associated 95% confidence limits for an independent variable with respect to a
	# specified reference value xref . 
	#
	# This function has 12 arguments which are described as follows:
	#
	# 	data:       Original S-Plus data set;
	# 	variable:   Name of the predictor variable;
	# 	fit:        The fitted object from the generalized linear model;
	# 	pos:        the position where the first coefficient of the predictor variable begins in GLM summary;
	#  	df :        number of degrees of freedom (df) for the predictor variable. If no specified knots are
	#  	            given, the program will select the df + 1 knots (including boundary knots) at suitably
	#  	            chosen quantiles of the predictor variable;
	#  	knots:      knots that define the natural spline, including boundary points. When knots and
	#  	            df are both given, the number of knots is equal to df + 1;
	# 	xref :      the reference predictor value;
        #  	typeout:    You can choose "MC" to output the Mean Change. The default is "MC"
	#  	step:       the step length to generate a vector of predictor values to get their MC 
	#                   as compared to xref . The smaller the step, the more continuous the graph
	#  	            looks, but the longer it will take to run the program. The default is 0.1;
	#  	xtitle:     the title of the x-axis. The default is the predictor variable (in uppercase);
	# 	ytitle:     the title of the y-axis. The default is the argument typeout (in uppercase);
	#       level:      The level for the confidence interval
        #  	na.action:  This parameter should be set to na.exclude if the original regression model excluded 
	#                   observations that were flagged as missing by na.exclude. The default is na.exclude, which means to 
	#                   exclude missing values from the computations. If na.action=na.fail, then na.fail will be applied,
        #                   which causes an error if missing values are present.
        #       col:        The color of lines
        #       lty:        The types of lines
        #       lwd:        The width of lines
	########################################################################################################
     
	##data = na.action(data)
        data = data_input

	# 	Get the variable column from data
	variable.vec <- data[, variable]
	# 	Get the range of variable
	variable.vec.range <- range(variable.vec, na.rm = TRUE)
	# 	if xref is out of range, a warning is given	
	if((xref < variable.vec.range[1]) | (xref > variable.vec.range[2])) {
		warning(paste("the reference value xref is not in the range of data variable"))
	}
	# Generate a sequence of points for which the Log RR as compared to the reference point xref will be calculated.
	xvec = seq(variable.vec.range[1], variable.vec.range[2], by = step)
	#  	if both df and knots are given, then the number of knots should be equal to df+1, otherwise
	# 	a warning will be given.
	##if(!missing(df) & !missing(knots)) {
	if(!is.null(df) & !is.null(knots)) {
		if((df + 1) != length(knots)) {
			warning(paste("df+1 is not equal to the number of knots"))
		}
	}
	# 	if knots is given instead of df, then df is equal to the number of knots - 1
	##if(missing(df) & (!missing(knots))) {
	if(is.null(df) & (!is.null(knots))) {
		df <- length(knots) - 1
	}
	# 	if df is given instead of knots, then the program will select the df + 1 knots (including boundary knots) at suitably
	#  	chosen quantiles of the predictor variable;
	##if(!missing(df) & (missing(knots))) {
	if(!is.null(df) & (is.null(knots))) {
		knots <- quantile(variable.vec, seq(0, 1, length = (df + 1)), na.rm = TRUE)
	}
	# 	if either df or knots is given, a warning is given.
	##if(missing(df) & missing(knots)) {
	if(is.null(df) & is.null(knots)) {
		warning(paste("please specify df or knots"))
	}
	# 	get the begining position of basis functions
	begin.pos <- pos
	# 	get the end position of basis functions
	end.pos <- pos + df - 1
	# 	get the estimated coefficients
	coeff <- summary(fit)$tTable[begin.pos:end.pos, 1]
	# 	get the correlation matrix of the estimated coefficients      
	corr <- summary(fit)$corFixed[begin.pos:end.pos, begin.pos:end.pos]
	#	 get the standard deviations of the estimated coefficients
	std <- summary(fit)$tTable[begin.pos:end.pos, 2]
	# 	 compute the  covariance matrix of estimated coefficients of basis 
	#        natural splines.
	var <- matrix(0, df, df)
	for(i in 1:df) {
		for(j in i:df) {
			var[i, j] <- std[i] * std[j] * corr[i, j]
		}
	}
	for(i in 1:df) {
		for(j in 1:(i - 1)) {
			var[i, j] <- var[j, i]
		}
	}
	# evaluate the basis functions at the reference point xref
	##basis0 <- ns(xref, knots = knots, Boundary.knots = TRUE)
	basis0 <- ns(xref, knots = knots[-c(1,df+1)], Boundary.knots = c(variable.vec.range[1],variable.vec.range[2]))
	# extract the length of the vector needed to save MC.
	nx <- length(xvec)
	# Initialize the vectors of the MC and lower and upper bound of confidence interval of MC
	mc <- rep(0, nx)
	lower.mc <- rep(0, nx)
	upper.mc <- rep(0, nx)
	# calculate the MC and confidence interval for each point
	for(i in 1:nx) {
		# evaluate the basis functions at the proposed point x	
		##basis1 <- ns(xvec[i], knots = knots, Boundary.knots = TRUE)
		basis1 <- ns(xvec[i], knots = knots[-c(1,df+1)], Boundary.knots = c(variable.vec.range[1],variable.vec.range[2]))
		# evaluate the difference of the bases functions at the 
		# proposed point x from the reference point xref 
		diff <- basis1 - basis0
		# evaluate the variance of mean change
		##varmc <- t(diff) %*% var %*% diff
		varmc <- t(as.vector(diff)) %*% var %*% as.vector(diff)
		# evaluate the mean change
		##mc[i] <- t(diff) %*% coeff
		mc[i] <- t(as.vector(diff)) %*% coeff
		# calculate the lower and upper bound of confidence interval of MC
                normstat <- qnorm(1-(1-level)/2)
		lower.mc[i] <- mc[i] - normstat * sqrt(varmc)
		upper.mc[i] <- mc[i] + normstat * sqrt(varmc)
	}

	# output the graph of mean change and the associate confidence interval if typeout is "MC"
        plot(xvec, mc, ylim = c(min(lower.mc), max(upper.mc)), type = lty, xlab = xtitle, ylab = ytitle, col=col,lwd=lwd)
	lines(xvec, lower.mc, type = lty,lwd=lwd)
	lines(xvec, upper.mc, type = lty,lwd=lwd)
	abline(h = 0)
	abline(v = xref)
}