R/plot.thregI.R

In thregI: Threshold Regression for Interval-Censored Data with a Cure Rate Option

### The following is the definition of the plot.thregI function
"plot.thregI" <-
function (x,var,scenario,graph,nolegend=0,nocolor=0,...)
{
	para <- match.call(expand.dots = FALSE)
  indx <- match(c("object", "var", "timevalue", "scenario"), names(para), nomatch=0)

	if (missing(x)) stop("An x argument is required")
	if (missing(var)) stop("A var argument is required")
	if (missing(graph)) stop("A graph argument is required")

  timevalue<- sort(model.extract(x$mf, "response")[,1])
	names(timevalue)<-NULL
	m_scenario<-match(c("scenario"), names(para), 0)
	m_hr<-match(c("var"), names(para), 0)
	m_graph<-match(c("graph"), names(para), 0)
  scenario_value<-NULL
  scenario_covariate<-NULL
  if(m_scenario!=0) {
	cur_scenario_string<-para[[m_scenario]]
	while(length(cur_scenario_string)>=2) {
		if(length(cur_scenario_string)==3) {
		 	current_scenario_covariate_string<-cur_scenario_string[[3]]
				if(length(current_scenario_covariate_string)==2) {
				  current_scenario_covariate<-current_scenario_covariate_string[[1]]
				  current_covariate_value<-current_scenario_covariate_string[[2]]
				  scenario_covariate<-c(current_scenario_covariate,scenario_covariate)
				  scenario_value<-c(current_covariate_value,scenario_value)
				}
		else {
				  stop("wrong scenario specification")
		  	}
				cur_scenario_string<-cur_scenario_string[[2]]
			}
			else if(length(cur_scenario_string)==2) {
				current_scenario_covariate_string<-cur_scenario_string
				current_scenario_covariate<-current_scenario_covariate_string[[1]]
				current_covariate_value<-current_scenario_covariate_string[[2]]

				scenario_covariate<-c(current_scenario_covariate,scenario_covariate)
				scenario_value<-c(current_covariate_value,scenario_value)

				current_scenario_covariate_string<-NULL
				cur_scenario_string<-NULL
				current_scenario_covariate<-NULL
				current_covariate_value<-NULL
			}
		}

        }
  #add intercept term
	scenario_covariate<-c("(Intercept)",scenario_covariate)
	scenario_value<-c(1,scenario_value)

	names(scenario_value)<-scenario_covariate
  m_hr<-match(c("var"), names(para), 0)
	hr_var<-as.character(para[[m_hr]])
  if (!is.factor(x$mf[[hr_var]])) stop("The variable for the hazard ratio calculation should be a factor variable")
  if (length(levels(x$mf[[hr_var]]))<2) stop("The variable for the hazard ratio calculation should have more than one level")
  if (!any(names(x$mf)==hr_var)) stop("The variable ",as.character(para[[m_hr]]) , " for the hazard ratio calculation should be included in the model")

  scenario_value_hr<-rbind(rep(0,length(levels(x$mf[[hr_var]]))-1),diag(rep(1,length(levels(x$mf[[hr_var]]))-1)))
	#names(scenario_value_hr)<-paste(hr_var,levels(x$mf[[hr_var]])[-1],sep="")
  colnames(scenario_value_hr) <- paste(hr_var,levels(x$mf[[hr_var]])[-1],sep="")
  if (any(matrix(rep(names(scenario_value),times=length(dimnames(scenario_value_hr)[[2]])),nrow=length(dimnames(scenario_value_hr)[[2]]),byrow=TRUE)==dimnames(scenario_value_hr)[[2]])) stop("Don't specify the scenario value of the dummy variable for the hazard ratio calculation!")
  scenario_value<-matrix(rep(scenario_value,times=length(levels(x$mf[[hr_var]]))),nrow=length(levels(x$mf[[hr_var]])),byrow=TRUE)
  dimnames(scenario_value)[[2]] <- as.list(scenario_covariate)
	scenario_value<-cbind(scenario_value,scenario_value_hr)

  #judge if any covariate for lny0 lacks of scenario value
	judgematrix_lny0<-matrix(rep(x$lny0,times=dim(scenario_value)[2]),ncol=dim(scenario_value)[2])==matrix(rep(dimnames(scenario_value)[[2]],times=length(x$lny0)),nrow=length(x$lny0),byrow=TRUE)
	if (any(!apply(judgematrix_lny0,1,any)))
  stop(paste("scenario value for",  x$lny0[!apply(judgematrix_lny0,1,any)][1], "is required!"))

	#choose the right position of scenario values to combine
	positionpick_lny0=judgematrix_lny0%*%c(1:dim(scenario_value)[2])

  lny0<-scenario_value[,positionpick_lny0]%*%x$coef[1:length(x$lny0)]
	y0<-exp(lny0)


  #judge if any covariate for mu lacks of scenario value
	judgematrix_mu<-matrix(rep(x$mu,times=dim(scenario_value)[2]),ncol=dim(scenario_value)[2])==matrix(rep(dimnames(scenario_value)[[2]],times=length(x$mu)),nrow=length(x$mu),byrow=TRUE)

  if (any(!apply(judgematrix_mu,1,any)))
  stop(paste("scenario value for",  x$mu[!apply(judgematrix_mu,1,any)][1], "is required!"))

	#choose the right position of scenario values to combine
	positionpick_mu=judgematrix_mu%*%c(1:dim(scenario_value)[2])

  mu<-scenario_value[,positionpick_mu]%*%x$coef[(length(x$lny0)+1):(length(x$lny0)+length(x$mu))]
	dim(y0)<-NULL
	dim(lny0)<-NULL
	dim(mu)<-NULL
  y0<-c(rep(y0,length(timevalue)))
	lny0<-c(rep(lny0,length(timevalue)))
	mu<-c(rep(mu,length(timevalue)))
	lenth_timevalue<-length(timevalue)

	#replicate values for matrix calculation
	timevalue_rep<-rep(timevalue,each=length(dimnames(scenario_value_hr)[[2]])+1)

	f<-exp((lny0-.5*(log(2*pi*(timevalue_rep^3))+(y0+mu*timevalue_rep)^2/timevalue_rep)))
  S<-exp(log(pnorm((mu*timevalue_rep+y0)/sqrt(timevalue_rep))-exp(-2*y0*mu)*pnorm((mu*timevalue_rep-y0)/sqrt(timevalue_rep))))
	h<-f/S

	if(m_graph!=0) 	graph_type<-para[[m_graph]]
	if(graph_type=="ds"){
  dim(f)<-c(length(dimnames(scenario_value_hr)[[2]])+1,lenth_timevalue)
	dt_graph<-cbind(timevalue,t(f))
	y_label="Estimated f(t)"
	y_lim=c(0,max(f, na.rm = T))
	legend_pos="topright"
	}
	else if(graph_type=="sv"){
	        dim(S)<-c(length(dimnames(scenario_value_hr)[[2]])+1,lenth_timevalue)
        	dt_graph<-cbind(timevalue,t(S))
		y_label="Estimated S(t)"
		y_lim=c(0,max(S, na.rm = T)+0.05)
		legend_pos="topright"
		#legend_pos="topleft"
	}
	else if(graph_type=="hz"){
	        dim(h)<-c(length(dimnames(scenario_value_hr)[[2]])+1,lenth_timevalue)
        	dt_graph<-cbind(timevalue,t(h))
		y_label="Estimated h(t)"
		y_lim=c(0,max(h, na.rm = T))
		legend_pos="topright"
		#legend_pos="topleft"
	}
	if(nocolor==0) {
		matplot(dt_graph[order(timevalue),1],dt_graph[order(timevalue),2:(length(dimnames(scenario_value_hr)[[2]])+2)],type = "l",lty=1:(length(dimnames(scenario_value_hr)[[2]])+1),xlab="time",ylab=y_label,ylim=y_lim)
	}
	else {
		matplot(dt_graph[order(timevalue),1],dt_graph[order(timevalue),2:(length(dimnames(scenario_value_hr)[[2]])+2)],type = "l",col=1,lty=1:(length(dimnames(scenario_value_hr)[[2]])+1),xlab="time",ylab=y_label,ylim=y_lim)
	}

	if(nolegend==0) {
	if(nocolor==0) {
	legend(legend_pos,paste(hr_var,"=",as.character(levels(x$mf[[hr_var]]))),col=1:(length(levels(x$mf[[hr_var]]))),lty=1:(length(levels(x$mf[[hr_var]]))),bty="n")
		}
		else {
		legend(legend_pos,paste(hr_var,"=",as.character(levels(x$mf[[hr_var]]))),col=1,lty=1:(length(levels(x$mf[[hr_var]]))),bty="n")
		}
	}
}