R/testContourPt.r

In Weibull-R/contour: Likelihood Ratio Contours

testContourPt<-function(x,  dist="weibull", CL=0.9,DF=1,MLEfit=NULL, RadLimit=1e-5, theta)  {

## check basic parameters of x
	if(class(x)!="data.frame") {stop("mlefit takes a structured dataframe input, use mleframe")}
	if(ncol(x)!=3)  {stop("mlefit takes a structured dataframe input, use mleframe")}
	xnames<-names(x)
	if(xnames[1]!="left" || xnames[2]!="right"||xnames[3]!="qty")  {
		 stop("mlefit takes a structured dataframe input, use mleframe")  }
## test for any na's and stop, else testint below will be wrong

## It turns out that this code is general to all fitting methods:
	if(tolower(dist) %in% c("weibull","weibull2p","weibull3p")){
		fit_dist<-"weibull"
	}else{
		if(tolower(dist) %in% c("lnorm", "lognormal","lognormal2p", "lognormal3p")){
			fit_dist<-"lnorm"
		}else{
		## Note:  only lslr contains experimental support for "gumbel"
			stop(paste0("dist argument ", dist, "is not recognized for mle fitting"))
		}
	}

## initialize counts at zero, to be filled as found
	Nf=0
	Ns=0
	Nd=0
	Ni=0

## need this length information regardless of input object formation

	failNDX<-which(x$right==x$left)
	suspNDX<-which(x$right<0)
	Nf_rows<-length(failNDX)
	if(Nf_rows>0) {
		Nf<-sum(x[failNDX,3])
	}
	Ns_rows<-length(suspNDX)
	if(Ns_rows>0) {
		Ns<-sum(x[suspNDX,3])
	}
	discoveryNDX<-which(x$left==0)
	Nd_rows<-length(discoveryNDX)
	if(Nd_rows>0) {
		Nd<-sum(x[discoveryNDX,3])
	}
	testint<-x$right-x$left
	intervalNDX<-which(testint>0)
	interval<-x[intervalNDX,]
	intervalsNDX<-which(interval$left>0)
	Ni_rows<-length(intervalsNDX)
	if(Ni_rows>0) {
		Ni<-sum(x[intervalsNDX,3])
	}

## rebuild input vector from components, because this order is critical
	fsiq<-rbind(x[failNDX,], x[suspNDX,], x[discoveryNDX,], interval[intervalsNDX,])

## now form the arguments for C++ call
## fsdi is the time vector to pass into C++
	fsd<-NULL
	if((Nf+Ns)>0)  {
		fsd<-fsiq$left[1:(Nf_rows + Ns_rows)]
	}
	if(Nd>0) {
		fsd<-c(fsd,fsiq$right[(Nf_rows + Ns_rows + 1):(Nf_rows +  Ns_rows + Nd_rows)])
	}
	if(Ni>0)  {
		fsdi<-c(fsd, fsiq$left[(Nf_rows + Ns_rows + Nd_rows + 1):nrow(fsiq)],
		fsiq$right[(Nf_rows + Ns_rows + Nd_rows + 1):nrow(fsiq)])
	}else{
		fsdi<-fsd
	}

	q<-fsiq$qty
## third argument will be c(Nf,Ns,Nd,Ni)
	N<-c(Nf_rows,Ns_rows,Nd_rows,Ni_rows)

## establish distribution number
	if(fit_dist=="weibull"){
		dist_num=1
	}else{
		if(fit_dist=="lnorm"){
			dist_num=2
		}else{
			stop("distribution not resolved for mle fitting")
		}
	}

	MLEclassList<-list(fsdi=fsdi,q=q,N=N)

## start of main contour procedures
	if(is.null(MLEfit))  {
		require(WeibullR)
## in this case the x argument is already an mleframe
		MLEfit<-unname(mlefit(x))
	}else{
		unname(MLEfit)
	}
##		Beta_hat<-MLEfit[2]
##		Eta_hat<-MLEfit[1]
## par is provided as a vector c(shape, scale)
	par_hat <- c(MLEfit[2], MLEfit[1])
	MLLx<-MLEfit[3]

	FF<-1
	ratioLL  <-  MLLx- qchisq(CL,DF)/(2*FF)


## Test for successful identification of a contour point at angle theta
	testValue <- .Call("testContourPt", MLEclassList, par_hat, dist_num, MLLx, ratioLL, RadLimit, theta, package="contour")

	testValue
}