R/code.r

Defines functions cable.change.conf find.fisher grad alph.vect alpha stick.ar.0 cable.ar.p.resid cable.Sig.ar.p cable.ar.p.iter pick.peak cable.fit.known.change cable.ar.0.fit bentcable.ar cable.diagnose.resid bentcable.dev.plot cable.lines cable.ar.p.diag cable.ar.p.plot cable.dev loglik.surf.ar.p loglik.given.g.ar.p cable.loop.ar.p sse.ar.p.given.g sse.ar.p.given.t.g sse.ar.p.given.beta.phi sse.ar.p sse.ar.p.core ar.p.resid.core is.stationary loglik.surf loglik.given.g cable.loop cable.design.mat fullcable.t basic.cable

Documented in alpha alph.vect ar.p.resid.core basic.cable bentcable.ar bentcable.dev.plot cable.ar.0.fit cable.ar.p.diag cable.ar.p.iter cable.ar.p.plot cable.ar.p.resid cable.change.conf cable.design.mat cable.dev cable.diagnose.resid cable.fit.known.change cable.lines cable.loop cable.loop.ar.p cable.Sig.ar.p find.fisher fullcable.t grad is.stationary loglik.given.g loglik.given.g.ar.p loglik.surf loglik.surf.ar.p pick.peak sse.ar.p sse.ar.p.core sse.ar.p.given.beta.phi sse.ar.p.given.g sse.ar.p.given.t.g stick.ar.0

#updated apr8, 2015

#-----------------------------------------
# regular cable fit deviance surface code
#-----------------------------------------

basic.cable<-function(z, g){

	if(g > 0)
		((z + g)^2/(4 * g)) * (z >=  - g & z <= g) + z * (z > g)
	else z * (z > 0)
}

fullcable.t<-function(t, b0, b1, b2, tau, gamma)
{
	b0 + b1 * t + b2 * basic.cable(t - tau, gamma)
}

cable.design.mat<-function(n,t.vect,tau,gamm)
{
	cbind(int=rep(1,n),t=t.vect,cable.t=basic.cable(t.vect-tau,gamm))
}


cable.loop<-function(i, j, tau.vect, gamm.vect, y.vect, design, max.flag = TRUE)
{
        t.hat <- tau.vect[i]
        g.hat <- gamm.vect[j]
        design[,3] <- basic.cable(design[,2] - t.hat, g.hat)
        fit <- try(lm(formula = y.vect ~ t + cable.t, data = design))

	#if(length(fit) == 12)  
   if(!inherits(fit,"try-error"))
       return( - dim(design)[1] * log(summary(fit)$sigma^2))
	else
		return(NA)
}

loglik.given.g<-function(i.vect, j, tau.vect, gamm.vect, y.vect,
	design, max.flag = TRUE){

		sapply(i.vect, cable.loop, j = j, tau.vect=tau.vect, gamm.vect=gamm.vect,
			y.vect=y.vect, design = design, max.flag = max.flag)
}

loglik.surf<-function(i.vect, j.vect, tau.vect, gamm.vect, y.vect,
	design, max.flag = TRUE){

	sapply(j.vect, loglik.given.g, i.vect = i.vect, tau.vect=tau.vect,
		gamm.vect=gamm.vect, y.vect=y.vect, design = design, max.flag = max.flag)
}

#---------------------------------------
# AR(p) cable fit deviance surface code
#---------------------------------------

is.stationary<-function(phi.vect){

	###############
	# stand-alone #
	########################################################
	# This function checks if the provided AR coefficients #
	# correspond to a stationary AR time series.           #
	########################################################

	#len<-length(try(arima.sim(n=2,list(ar=phi.vect)),silent=TRUE))

	#if(len==2)
   #   return(TRUE)
	#else
	#	return(FALSE)

    trial<-try(arima.sim(n=2,list(ar=phi.vect)),silent=TRUE)

    return(!inherits(trial,"try-error"))

}

ar.p.resid.core<-function(b0,b1,b2,tau,gamm,phi.vect,y.vect,t.vect,n,pp){

	if(pp==1)
		phi.vect<-as.matrix(phi.vect)

	err.vect<-y.vect-fullcable.t(t.vect,b0,b1,b2,tau,gamm)

	v.mat<-matrix(nrow=(n-pp),ncol=(pp+1))

	for(j in pp:0)
		v.mat[,(pp-j+1)]<-err.vect[(j+1):(n-pp+j)]
		
	innov<-as.vector(v.mat[,1]-v.mat[,-1]%*%phi.vect)

	return(list(resid=err.vect,innov=innov))

}

sse.ar.p.core<-function(b0,b1,b2,tau,gamm,phi.vect,y.vect,t.vect,n,pp){

	return(sum(
		ar.p.resid.core(b0,b1,b2,tau,gamm,phi.vect,y.vect,t.vect,n,pp)$innov^2
	))

}

sse.ar.p<-function(theta,y.vect,t.vect=NA,n=NA,pp=NA,stick=FALSE){

	# `optim()' doesn't like the variable `p'!

	b0<-theta[1]
	b1<-theta[2]
	b2<-theta[3]
	tau<-theta[4]

	if(!stick){

		gamm<-theta[5]
		phi.vect<-theta[-(1:5)]
	}

	else{

		gamm<-0
		phi.vect<-theta[-c(1:4)]
	}	

	if(is.na(n)){

		n<-length(y.vect)
		t.vect<-c(0:(n-1))

	}

	if(is.na(pp))
		pp<-length(phi.vect)

	return(sse.ar.p.core(b0,b1,b2,tau,gamm,phi.vect,y.vect,t.vect,n,pp))
}

sse.ar.p.given.beta.phi<-function(theta,beta.vect,phi.vect,y.vect,t.vect,n,pp){

	tau<-theta[1]
	gamm<-theta[2]
	b0<-beta.vect[1]
	b1<-beta.vect[2]
	b2<-beta.vect[3]

	return(sse.ar.p.core(b0,b1,b2,tau,gamm,phi.vect,y.vect,t.vect,n,pp))
}

sse.ar.p.given.t.g<-function(theta,tau,gamm,y.vect,t.vect,n,pp){

	b0<-theta[1]
	b1<-theta[2]
	b2<-theta[3]
	phi.vect<-theta[-(1:3)]

	return(sse.ar.p.core(b0,b1,b2,tau,gamm,phi.vect,y.vect,t.vect,n,pp))
}

sse.ar.p.given.g<-function(theta,gamm,y.vect,t.vect,n,pp){

	b0<-theta[1]
	b1<-theta[2]
	b2<-theta[3]
	tau<-theta[4]
	phi.vect<-theta[-(1:4)]

	return(sse.ar.p.core(b0,b1,b2,tau,gamm,phi.vect,y.vect,t.vect,n,pp))
}

cable.loop.ar.p<-function(tau,gamm,y.vect,t.vect,n,pp,fit=FALSE,verbose=FALSE)
{
	xreg<-cable.design.mat(n,t.vect,tau,gamm)[,-1]

	if(verbose)
		message(paste("tau=",tau,", gamma=",gamm))

	fitted<-arima(y.vect,order=c(pp,0,0),xreg=xreg,
		method="CSS",optim.control=list(maxit=5000))

	phi.hat<-fitted$coef[1:pp]
	beta.hat<-fitted$coef[-c(1:pp)]

	if(!fit){

		sse<-sse.ar.p.given.t.g(c(beta.hat,phi.hat),tau,gamm,y.vect,t.vect,n,pp)

		return(-(n-pp)*log(sse/(n-pp)))
	}

	else{

		return(fitted)

	}
}

loglik.given.g.ar.p<-function(tau.vect,gamm,y.vect,t.vect,n,pp){

	sapply(tau.vect,cable.loop.ar.p,gamm=gamm,
		y.vect=y.vect,t.vect=t.vect,n=n,pp=pp)
}

loglik.surf.ar.p<-function(
	tau.vect,gamma.vect,y.vect,t.vect=NA,n=NA,pp=NA){

	if(is.na(n)){

		n<-length(y.vect)
		t.vect<-c(0:(n-1))

	}

	sapply(gamma.vect,loglik.given.g.ar.p,
		tau.vect=tau.vect,y.vect=y.vect,t.vect=t.vect,n=n,pp=pp)

}

cable.dev<-function(tau.vect,gamma.vect,y.vect,t.vect=NULL,p=0){

	###############
	# stand-alone #
	####################################################################
	# 'tau.vect' and 'gamma.vect' specify a tau-gamma grid;            #
	# 'p' is the AR order for bent-cable residuals.                    #
	# This function computes the bent-cable profile deviance surface   #
	# matrix over the specified tau-gamma grid. It is intended to be   #
	# used for plotting a contour or perspective plot of the deviance  #
	# surface, in which case 'tau.vect' and 'gamma.vect' should have   #
	# length >=2 so that returned object is a matrix. If               #
	# such a plot is not required, then 'tau.vect' and 'gamma.vect'    #
	# can be scalar.                                                   #
	#                                                                  #
	# WARNING:                                                         #
	# When p>0, time series data are assumed, and 't.vect' should be   #
	# equidistant with unit increments.                                #
	####################################################################

	n<-length(y.vect)

	if(is.null(t.vect))
		t.vect<-c(0:(n-1))

	if(p==0){

		design<-data.frame(cable.design.mat(n,t.vect,0,0))

		loglik.mat<-loglik.surf(c(1:length(tau.vect)),c(1:length(gamma.vect)),
			tau.vect, gamma.vect, y.vect, design,TRUE)
	}

	else{

		message("Please be patient...")

		loglik.mat<-loglik.surf.ar.p(tau.vect,gamma.vect,y.vect,t.vect,n,p)
	}

	return(loglik.mat-max(na.omit(loglik.mat)))

}

#---------------------------
# plots code
#---------------------------

cable.ar.p.plot<-function(ar.p.fit,xlab="time",ylab="",
	main=NULL,ctp.ci=NULL){

	###############
	# stand-alone #
	#####################################################
	# This function plots the time series data for a    #
	# 'cable.ar.p.iter()' object provided by the user   #
	# as 'ar.p.fit'. This object must be obtained for   #
	# AR(p) data with p>0. Superimposed on the data is  #
	# the bent-cable fit from the object. The transition#
	# tau +/- gamma and tau are marked in red.          #
	# An optional 'cable.change.conf()' object as       #
	# 'ctp.ci' is also superimposed in blue.            #
	#####################################################


	y.vect<-ar.p.fit$y
	p<-ar.p.fit$p
	t.vect<-ar.p.fit$t

	main<-paste(main,": AR",p,"fit")

	plot(t.vect,y.vect,type="l",xlab=xlab,ylab=ylab,main=main)

	if(!ar.p.fit$stick)
		cable.lines(t.vect,ar.p.fit$est[1:5],col="red",lty=1)
	else
		cable.lines(t.vect,c(ar.p.fit$est[1:4],0),col="red",lty=1)

	if(!is.null(ctp.ci))
		abline(v=c(ctp.ci$change.hat,ctp.ci$interval),lty=2,col="blue")

}

cable.ar.p.diag<-function(ar.p.fit,resid.type="p",xlab="time",ylab="",
	main=NULL,main.all=NULL,ctp.ci=NULL){

	###############
	# stand-alone #
	#####################################################
	# This function produces (P)ACF diagnostic plots    #
	# for a 'cable.ar.p.iter()' object provided by the  #
	# user as 'ar.p.fit'. It plots the time series data #
	# from the object, and superimposes on it the bent- #
	# cable fit from the object. The transition tau +/- #
	# gamma and tau are marked in red.                  #
	# An optional 'cable.change.conf()' object as       #
	# 'ctp.ci' is also superimposed in blue. Plot title #
	# 'main' is for the time series plot, and           #
	# 'main.all' is for the entire set of plots.        #
	#####################################################

	par(mfcol=c(4,2))

	y.vect<-ar.p.fit$y
	t.vect<-ar.p.fit$t
	stick<-ar.p.fit$stick
	p<-ar.p.fit$p

	if(p==0){
		stop("This function cannot handle p=0.")
		return()
	}

	if(!stick)
		gamm<-ar.p.fit$est[5]
	else
		gamm<-0

	bend<-ar.p.fit$est[4]+c(-1,0,1)*gamm

	plot(1:2,type="n",bty="n",xaxt="n",yaxt="n",xlab="",ylab="")
	text(1.5,1.5,paste(main.all,"\n",main),cex=3)

	cable.ar.p.plot(ar.p.fit,xlab=xlab,ylab=ylab,main=main,ctp.ci=ctp.ci)

	fit.resid<-cable.ar.p.resid(ar.p.fit)

	plot(t.vect,fit.resid$resid,xlab=xlab,ylab="",
		main="fitted residuals",type=resid.type)
	abline(h=0,col="red")
	abline(v=bend,col="red")
	if(!is.null(ctp.ci))
		abline(v=c(ctp.ci$change.hat,ctp.ci$interval),lty=2,col="blue")

	plot(t.vect,c(rep(NA,p),fit.resid$innov),xlab=xlab,ylab="",
		main="fitted innovations",type=resid.type)
	abline(h=0,col="red")
	abline(v=bend,col="red")
	if(!is.null(ctp.ci))
		abline(v=c(ctp.ci$change.hat,ctp.ci$interval),lty=2,col="blue")

	acf(fit.resid$resid,main="ACF for Fitted Residuals")
	pacf(fit.resid$resid,main="PACF for Fitted Residuals")

	acf(fit.resid$innov,main="ACF for Fitted Innovations")
	pacf(fit.resid$innov,main="PACF for Fitted Innovations")

	par(mfrow=c(1,1))

}

cable.lines<-function(x,theta,col="black",lwd=1,lty=2,fit.lty=1){

	###############
	# stand-alone #
	####################################################
	# This function is meant to superimpose a fitted   #
	# bent cable over the plot of response-covariate   #
	# data the exhibit a bent cable structure.         #
	#                                                  #
	# 'x' contains the design points for the existing  #
	# plot, or the range of the design points.         #
	# This function then adds to it a bent cable with  #
	# coefficients 'theta'.                            #
	# Superimposed on the cable is the transition      #
	# marked by tau and tau +/- gamma. 'lwd' is for    #
	# the entire plot, 'lty' is for the transition,    #
	# and 'fit.lty' is for the cable.                  #
	####################################################

	changes<-theta[4]+c(-1,0,1)*theta[5]

	segments(min(x)-1,fullcable.t(min(x)-1,theta[1],theta[2],theta[3],
		theta[4],theta[5]),
		changes[1],fullcable.t(changes[1],theta[1],theta[2],theta[3],
		theta[4],theta[5]),col=col,lwd=lwd,lty=fit.lty)

	segments(changes[3],fullcable.t(changes[3],theta[1],theta[2],theta[3],
		theta[4],theta[5]),
		max(x)+1,fullcable.t(max(x)+1,theta[1],theta[2],theta[3],
		theta[4],theta[5]),col=col,lwd=lwd,lty=fit.lty)

	middle<-seq(changes[1],changes[3],range(changes)%*%c(-1,1)*.01)

	lines(middle,fullcable.t(middle,theta[1],theta[2],theta[3],
		theta[4],theta[5]),col=col,lwd=lwd,lty=fit.lty)

	abline(v=changes,col=col,lwd=lwd,lty=lty)

}


bentcable.dev.plot<-function(tau.vect,gamma.vect=NULL,y.vect,t.vect=NULL,stick=FALSE,p=0){

	##################
	# main interface #
	####################################################################
	# 'tau.vect' and 'gamma.vect' specify a tau-gamma grid.            #
	# 'p' is the AR order for bent-cable residuals.                    #
	# for 'stick=F': this function plots the contours of and returns   #
	#	the bent-cable profile deviance surface matrix over the         #
	#	specified tau-gamma grid - 'tau.vect' and 'gamma.vect'          #
	#	must have the same length                                       #
	# for 'stick=T': it plots and returns the broken-stick profile     #
	#	deviance function over the specified tau-domain while           #
	#	taking gamma=0 and thus ignores 'gamma.vect'                    #
	####################################################################

	n<-length(y.vect)

	if(is.null(t.vect))
		t.vect<-c(0:(n-1))

	if(!stick){

		dev<-cable.dev(tau.vect,gamma.vect,y.vect,t.vect,p)

		contour(tau.vect,gamma.vect,dev)
	}

	else{

		if(!is.null(gamma.vect))
			message("Ignoring 'gamma.vect'...")

		gamma.vect<-0

		dev<-cable.dev(tau.vect,gamma.vect,y.vect,t.vect,p)

		plot(tau.vect,dev,type="l",main="Broken-Stick Deviance",xlab="tau")

	}


	return(list(dev=dev,tau=tau.vect,gamma=gamma.vect))

}


cable.diagnose.resid<-function(fit){

	#plots PACF for residuals from "fit"
	#returns AR order based on AIC

	fit.resid<-resid(fit)

	pacf(fit.resid)

	return(c(p.aic=length(ar(fit.resid)$ar)))

}



#---------------------------
# AR(p) iterative fit code
#---------------------------

bentcable.ar<-function(y.vect,tgdev=NULL,p=0,stick=FALSE,t.vect=NULL,
	init.cable=NULL,init.phi=NULL,tol=1e-4,
	method0="css",method1="yw",ci.level=.95,main=NULL){

	##################
	# main interface #
	##################	
	
	# If 'init.cable' is missing then 'dev.mat'	
	# must be provided so as to produce an initial AR(p) fit.

	# 'method0' and 'method1' can be "css" (conditional sum-of-squares),
	# "mle", or "yw" (Yule-Walker). 'method0' is the default fitting method;
	# if 'method0' fails, then the fit is repeated using 'method1'.

	#####################################################################
	# This function fits a bent cable to response data 'y.vect' with    #
	# covariate values in 't.vect' (optional). For time-series data, a  #
	# stationary AR(p) autocorrelation structure is assumed. By         #
	# default, time points start at t=0, with increments of 1. Although #
	# the user can supply an alternative 't.vect', caution should be    #
	# taken in this case (see below).                                   #
	#                                                                   #
	# The fitted bent cable and estimated transition (tau and tau +/-   #
	# gamma) are plotted in red over the original response data. (P)ACF #
	# diagnostics are also plotted for fits with p>0. A confidence      #
	# interval for the critical time point (CTP) is also computed -     #
	# assuming the default time scale - and plotted in blue. By         #
	# definition, the CTP exists only when a change of sign in the      #
	# cable's slope takes place during transition. If the CTP (almost)  #
	# does not exist, then the CTP computation will fail. Failure can   #
	# also occur for a user-supplied time vector - even if successfully #
	# computed, the resulting CTP and confidence interval would be      #
	# meaningless as the computation assumes the default time scale.    #
	# The user is advised to rely on the default, and, if necessary,    #
	# transform the results to the preferred time scale after the model #
	# has been fitted.                                                  #
	#                                                                   #
	# Exception to this nuisance is in a non-time-series context where  #
	# the response data are independent, and the covariate is not time  #
	# and/or the design points are not equidistant. Then the "initial   #
	# AR(0) cable" is actually the max. likelihood fit for the original #
	# iid bent cable model (Chiu et al. 2006). No confidence intervals  #
	# are computed in this case.                                        #
	#                                                                   #
	# WARNING:                                                          #
	# The bent-cable likelihood surface often has multiple peaks. The   #
	# user is advised to examine the likelihood value to decide if the  #
	# returned fit corresponds to a local maximum. For a conditional    #
	# max. likelihood fit, the conditional sum-of-squares error (SSE)   #
	# can replace the likelihood for this purpose. If a 'css' fit is    #
	# successful, the conditional SSE is stored in                      #
	# '$cable$ar.p.fit$value' from the returned object. If 'css' fails  #
	# and the function refits using 'yw' or 'mle', then the conditional #
	# SSE is not directly retrievable. For model comparisons purposes,  #
	# the user can use the 'yw' or 'mle' estimate as the starting value #
	# for a subsequent 'css' fit, and the conditional SSE will appear   #
	# in the screen output as "initial value" while the 'css' algorithm #
	# iterates.                                                         #
	#####################################################################

	if(!is.null(tgdev) & (!is.null(init.cable) | !is.null(init.phi))){
		stop("'tgdev' cannot be used together with 'init.cable' or
			'init.phi'.")
		return()
	}

	if(p>0){

		if(is.null(init.phi))
			init.phi<-rep(c(.5,-.5),p)[1:p]

		else
			if(length(init.phi)!=p){

				message("****************************************")
				message("Your 'init.phi' and 'p' don't match:")
				message("ignoring 'init.phi', using given 'p'...")
				message("****************************************")

				init.phi<-rep(c(.5,-.5),p)[1:p]

			}

	}

	else

		if(length(init.phi)>0){

			message("****************************************")
			message("Your 'init.phi' and 'p' don't match:")
			message("ignoring 'p', using given 'init.phi'...")
			message("****************************************")

			p<-length(init.phi)

		}

	n<-length(y.vect)

	if(is.null(t.vect))
		t.vect<-c(0:(n-1))

	wrong.t<-FALSE

	if(sum(t.vect!=c(0:(n-1)))>0){

		message("*******************************")
		message("WARNING:")
		message("Cannot estimate CTP correctly")
		message("unless 't.vect' is c(0,1,2,...)")
		message("*******************************")

		wrong.t<-TRUE

	}

	if(!stick){ # fit bent cable:

		if(is.null(init.cable)){ # just getting initial values

			if(is.null(tgdev)){
				stop("Must supply 'init.cable' or 'tgdev'!")
				return()
			}

			message("*****************************************")
			message(paste("Finding initial values for AR(",p,") cable...",sep=""))
			message("*****************************************")

			if(p==0){

				message("*******************************")
				message("Finding best AR(0) cable fit...")
				message("*******************************")

				initfit<-try(cable.ar.0.fit(y.vect,t.vect,tgdev$tau,tgdev$gamma,tgdev$dev))

            #if(length(initfit$fit)!=6){
            if(inherits(initfit$fit,"try-error")){

					message("*****************************************************")
					message("AR(0) fit failed (deviance surface too irregular).")
					message("You may still try to use the printed values above")
					message("(if any) as initial values for subsequent AR(p) fits.")
					message("*****************************************************")
					return()

				}

				message("******************************************")
				message("If you have time-series data, then")
				message("choose your p according to AIC and PACF...")
				message("******************************************")

				print(cable.diagnose.resid(initfit$fit))

				message("************************************************")
				message("Please apply 'coef()[1:5]' to the '$fit' element")
				message("of the returned AR(0) fit as initial values for")
				message("subsequent AR(p) fits.")
				message("************************************************")

				return(initfit)

			}

			else{

				initfit<-try(cable.fit.known.change(y.vect,t.vect,n,tgdev$tau,tgdev$gamma,tgdev$dev,p))

				#if(length(initfit$fit)!=13){
            if(inherits(initfit$fit,"try-error")){

					message("*************************")
					message(paste("AR(",p,") fit failed.",sep=""))
					message("*************************")
					return()
				}

				message("**************************************")
				message("Please use the '$init' element of the")
				message("returned fit as initial values for")
				message(paste("subsequent AR(",p,") fits.",sep=""))
				message("**************************************")

				return(initfit)
				
			}
		}

		else{ # ready for real fit

			if(length(init.cable)>5){

				message("Using only 'init.cable[1:5]'...")
				init.cable<-init.cable[1:5]
			}

			init.vect<-c(init.cable,init.phi)

			fit<-cable.ar.p.iter(init.vect,y.vect,t.vect,n,tol,method0,method1)

			if(p==0){

				plot(t.vect,y.vect,xlab="",ylab="")
				title(paste("AR(0) bent cable:\n",main))
				cable.lines(t.vect,coef(fit$fit)[1:5],col="red")

			}

			if(!wrong.t)
				ctp.ci<-try(cable.change.conf(fit,ci.level))
			else
				ctp.ci<-NULL

			if(length(ctp.ci)!=3){

				message("******************************************")
				message("Failed to compute CTP confidence interval!")
				message("******************************************")

				if(p>0)
					cable.ar.p.diag(fit,main=main)

				return(list(cable=fit))
			}	

			else{

				if(p>0)
					cable.ar.p.diag(fit,main=main,ctp.ci=ctp.ci)
				else
					abline(v=c(ctp.ci$change,ctp.ci$int),lty=2,col="blue")

				return(list(cable=fit,ctp=ctp.ci))
			}
		}
	}

	# fit broken stick:


	if(is.null(init.cable)){ # just getting intial values

		if(is.null(tgdev)){
			stop("Must supply 'init.cable' or 'tgdev'!")
			return()
		}

		message("*****************************************")
		message(paste("Finding initial values for AR(",p,") stick...",sep=""))
		message("*****************************************")

		if(length(tgdev$gamma)>1 | tgdev$gamma[1]!=0){
			stop("Your 'tgdev' did not come from a stick fit.")
			return()
		}

		if(p==0){

			message("*******************************")
			message("Finding best AR(0) stick fit...")
			message("*******************************")

			initfit<-try(cable.ar.0.fit(y.vect,t.vect,tgdev$tau,tgdev$gamma,tgdev$dev,TRUE))
         
         #if(length(initfit$fit)!=6){
			if(inherits(initfit$fit,"try-error")){

				message("*****************************************************")
				message("AR(0) fit failed (deviance too irregular).")
				message("You may still try to use the printed values above")
				message("(if any) as initial values for subsequent AR(p) fits.")
				message("*****************************************************")
				return()

			}

			message("******************************************")
			message("If you have time-series data, then")
			message("choose your p according to AIC and PACF...")
			message("******************************************")

			message(cable.diagnose.resid(initfit$fit))

			message("************************************************")
			message("Please apply 'coef()[1:4]' to the '$fit' element")
			message("of the returned AR(0) fit as initial values for")
			message("subsequent AR(p) fits.")
			message("************************************************")

			return(initfit)

		}

		initfit<-try(cable.fit.known.change(y.vect,t.vect,n,tgdev$tau,tgdev$gamma,tgdev$dev,p))
      
      #if(length(initfit$fit)!=13){
		if(inherits(initfit$fit,"try-error")){

			message("*************************")
			message(paste("AR(",p,") fit failed.",sep=""))
			message("*************************")
			return()
		}

		message("**************************************")
		message("Please use the '$init' element of the")
		message("returned fit as initial values for")
		message(paste("subsequent AR(",p,") fits.",sep=""))
		message("**************************************")

		return(initfit)

	}

	# ready for real fit

	if(length(init.cable)>4){

		message("Using only 'init.cable[1:4]'...")
		init.cable<-init.cable[1:4]
	}

	if(p==0){

		fit<-stick.ar.0(init.cable,y.vect,t.vect,n)

		plot(t.vect,y.vect,xlab="",ylab="")
		cable.lines(t.vect,c(coef(fit$fit)[1:4],0),col="red")
		title(paste("AR(0) broken stick:\n",main))

	}

	else{	# fit AR broken stick:

		init.vect<-c(init.cable,init.phi)

		fit<-cable.ar.p.iter(init.vect,y.vect,t.vect,n,tol,method0,method1,stick=TRUE)
	}

	if(!wrong.t)
		ctp.ci<-try(cable.change.conf(fit,ci.level))
	else
		ctp.ci<-NULL

	if(length(ctp.ci)!=3){

		message("******************************************")
		message("Failed to compute CTP confidence interval!")
		message("******************************************")

		if(p>0)
			cable.ar.p.diag(fit,main=main)

		return(list(cable=fit))
	}	

	else{

		if(p>0)
			cable.ar.p.diag(fit,main=main,ctp.ci=ctp.ci)
		else
			abline(v=c(ctp.ci$change,ctp.ci$int),lty=2,col="blue")

		return(list(cable=fit,ctp=ctp.ci))
	}

}


cable.ar.0.fit<-function(y.vect,t.vect=NULL,tau.vect,gamma.vect,dev.mat,stick=FALSE){

	###############
	# stand-alone #
	#####################################################
	# This function fits a bent cable to response data  #
	# 'y.vect' over optional covariate value 't.vect',  #
	# assuming iid errors. Unlike 'bentcable.ar()',     #
	# this function does not require explicit initial   #
	# parameter values, but determines them via         #
	# 'tau.vect', 'gamma.vect', and 'dev.mat', which    #
	# should be obtained from 'bentcable.dev.plot()'.   #
   #####################################################


	n<-length(y.vect)

	if(is.null(t.vect))
		t.vect<-c(0:(n-1))

	if(!is.matrix(dev.mat))
		dev.mat<-matrix(dev.mat)

	fit0<-cable.fit.known.change(y.vect,t.vect,n,tau.vect,gamma.vect,dev.mat)

	if(!stick)
		return(cable.ar.p.iter(fit0$init,y.vect,t.vect,n,tol=NA))

	else
		return(stick.ar.0(fit0$init[1:4],y.vect,t.vect,n))

}


cable.fit.known.change<-function(y.vect,t.vect=NULL,n=NA,tau.vect,gamma.vect,dev.mat,p=0){

	###############
	# stand-alone #
	######################################################
	# This function fits a bent cable to response data   #
	# 'y.vect' over optional covariate values 't.vect',  #
	# conditioned on a known transition (i.e. known tau  #
	# and gamma). The user supplies 'tau.vect',          #
	# 'gamma.vect', and 'dev.mat' as a result of         #
	# 'bentcable.dev.plot()'. Using these, this function #
	# locates the peak of the deviance surface over the  #
	# specified tau-gamma grid, and computes the bent-   #
	# cable fit at this grid point. If multiple peaks    #
	# exist (such as along a ridge), then only that at   #
	# the smallest tau and smallest gamma is used.       #
	#                                                    #
	# For 'p=0', iid errors are assumed, and 't.vect'    #
	# can be non-equidistant. For p>0, AR(p) errors are  #
	# assumed, and 't.vect' must be equidistant with     #
	# unit increments, or the returned values will be    #
	# meaningless.                                       #
	######################################################


	peak<-pick.peak(dev.mat==0,nrow(dev.mat),ncol(dev.mat))

	t0<-tau.vect[peak$row]
	g0<-gamma.vect[peak$col]

	if(is.na(n))
		n<-length(y.vect)

	if(is.null(t.vect))
		t.vect<-c(0:(n-1))

	if(p==0){

		design<-data.frame(cable.design.mat(n,t.vect,t0,g0))

		fit<-lm(formula=y.vect~t+cable.t,data=design)

		return(list(fit=fit,init=c(coef(fit),t0,g0)))
	}

	else{

		fit<-cable.loop.ar.p(t0,g0,y.vect,t.vect,n,p,TRUE)

		init<-coef(fit)
		init<-c(init[-c(1:p)],t0,g0,init[1:p])

		names(init)[1:4]<-c("b0","b1","b2","tau")

		if(g0>0)
			names(init)[5]<-"gamm"
		else
			init<-init[-5]

		return(list(fit=fit,init=init))

	}

}

pick.peak<-function(dev.mat,nrow,ncol){

	pick.row<-dev.mat%*%rep(1,ncol)
	pick.row<-c(1:nrow)[pick.row>0][1]	# take smallest value in case of ridge

	pick.col<-rep(1,nrow)%*%dev.mat
	pick.col<-c(1:ncol)[pick.col>0][1]	# take smallest value in case of ridge

	return(list(row=pick.row,col=pick.col))
}


cable.ar.p.iter<-function(init,y.vect,t.vect=NULL,n=NA,
	tol,method0="css",method1="yw",stick=FALSE){

	###############
	# stand-alone #
	#########################################################################
	# for 'stick=F':                                                        #
	#	'init' must be in the order: b0, b1, b2, tau, gamma, phi.vect;       #
	#	if phi is present (i.e. AR(p) with p>0), then 't.vect' must be       #
	#	equidistant with unit increments                                     #
	#                                                                       #
	# for 'stick=T':                                                        #
	#	'init' must be in the order: b0, b1, b2, tau, phi.vect; phi          #
	#	must be present and 't.vect' must be equidistant with unit           #
	#       increments                                                      #
	#                                                                       #
	# 'method0' and 'method1' can be "css" (conditional sum-of-squares),    #
	#   "mle" or "yw" (Yule-Walker)                                         #
	#                                                                       #
	# 'tol' is the tolerance for determining convergence                    #
	#                                                                       #
	# This function fits a bent cable to the response data 'y.vect' with    #
	# optional covariate values in 't.vect', assuming AR errors. The AR     #
	# order is determined from the dimension of 'init'.                     #
	#########################################################################

	if(is.na(n))
		n<-length(y.vect)

	if(is.null(t.vect))
		t.vect<-c(0:(n-1))

	if(!stick)
		pp<-length(init[-c(1:5)])
	else
		pp<-length(init[-c(1:4)])

	if(pp==0){

		if(stick){
			stop("This function cannot handle AR(0) broken sticks.")
			return()
		}

		message("Trying 'nls()' Gauss-Newton algorithm...")

		fit0<-try(nls(formula=y.vect~fullcable.t(t.vect,b0,b1,b2,tau,gamm),trace=TRUE,
			start=list(b0=init[1],b1=init[2],b2=init[3],tau=init[4],gamm=init[5])))

      #if(length(fit0)!=6){
		if(inherits(fit0,"try-error")){

			message("Trying 'nls()' Golub-Pereyra algorithm...")

			fit0<-try(nls(formula=y.vect~fullcable.t(t.vect,b0,b1,b2,tau,gamm),trace=TRUE,
				algorithm="plinear",
				start=list(b0=init[1],b1=init[2],b2=init[3],tau=init[4],
				gamm=init[5])))

         #if(length(fit0)!=6){
			if(inherits(fit0,"try-error")){

				message("Trying 'nls()' Port algorithm...")

				fit0<-try(nls(formula=y.vect~fullcable.t(t.vect,b0,b1,b2,tau,gamm),trace=TRUE,
					algorithm="port",
					start=list(b0=init[1],b1=init[2],b2=init[3],tau=init[4],
					gamm=init[5])))

            #if(length(fit0)!=6){
            if(inherits(fit0,"try-error")){
					stop("Initial fit unsuccessful. Please try alternative initial values.")
					return()
				}
			}		
		}

		message("Converged!")

		return(list(fit=fit0,y=y.vect,t=t.vect,n=n,p=pp,stick=stick))

	}

	fit<-optim(init,sse.ar.p,y.vect=y.vect,t.vect=t.vect,n=n,pp=pp,stick=stick,
		control=list(trace=TRUE,abstol=tol,maxit=5000),method="BFGS")

	if(method0=="css"){

		if(!stick)
			phi.hat<-fit$par[-c(1:5)]
		else
			phi.hat<-fit$par[-c(1:4)]

		if(is.stationary(phi.hat)){

			message("CSS successful")
			return(list(estimate=fit$par,ar.p.fit=fit,y=y.vect,t=t.vect,n=n,p=pp,
				stick=stick,method=method0))
		}

		else{

			message("Non-stationary CSS fit")
			method0<-method1
		}

	}

	if(method0!="css"){

		message(paste("Applying method",method0,"..."))

		if(!stick)
			start<-fit$par
		else
			start<-c(fit$par[1:4],0,fit$par[-c(1:4)])

		b0<-start[1]
		b1<-start[2]
		b2<-start[3]
		tau<-start[4]
		gamm<-start[5]

		error<-999
		ct<-1

		while(error>tol){

			fit0.resid<-y.vect-fullcable.t(t.vect,b0,b1,b2,tau,gamm)

			Sig.ar.p.out<-cable.Sig.ar.p(fit0.resid,pp,n=n,method=method0)
			Sig<-Sig.ar.p.out$Sig
			phi.vect<-Sig.ar.p.out$phi

			xreg<-cable.design.mat(n,t.vect,tau,gamm)
			X.tran.Sig.inv<-t(xreg)%*%solve(Sig)

			beta.given.change<-solve(X.tran.Sig.inv%*%xreg)%*%
				X.tran.Sig.inv%*%y.vect

			b0<-beta.given.change[1]
			b1<-beta.given.change[2]
			b2<-beta.given.change[3]

			change.given.beta<-optim(c(tau,gamm),sse.ar.p.given.beta.phi,
				beta.vect=beta.given.change,phi.vect=phi.vect,
				y.vect=y.vect,t.vect=t.vect,n=n,pp=pp,
				control=list(trace=FALSE,abstol=tol,maxit=5000),method="BFGS")$par
				# conditional least squares

			tau<-change.given.beta[1]

			if(!stick)
				gamm<-change.given.beta[2]
			else
				gamm<-0

			if(gamm<0){

				message(paste("gamma hat =",gamm,"< 0: forcing to 2*tol"))
				gamm<-2*tol

			}


			est<-c(b0,b1,b2,tau,gamm,phi.vect)

			error<-sum(abs(start-est))
			message(paste(ct,": vector difference ",round(error,digits=ceiling(abs(log10(tol)))+2)))

			start<-est

			ct<-ct+1

		}

		names(est)[1:5]<-c("b0","b1","b2","tau","gamm")

		if(stick)
			est<-est[-5]

		return(list(estimate=est,ar.p.fit=Sig.ar.p.out$fit,y=y.vect,t=t.vect,n=n,p=pp,
			stick=stick,method=method0))
	}
}

cable.Sig.ar.p<-function(w.vect,p,n=NA,method="yw"){

	if(is.na(n))
		n<-length(w.vect)

	ar.p.fit<-ar(w.vect,aic=FALSE,order.max=p,demean=TRUE,method=method)
		# demean in case w.vect is from a fit that doesn't go through data

	#--------------------------------------------------------------
	# WARNING: ar() with method="mle" can give convergence problems
	# w/in iteration(s), although overall iterations converge
	#--------------------------------------------------------------

	phi.vect<-ar.p.fit$ar

	switch(method,

		yw={

			autocor.vect<-acf(w.vect,plot=FALSE,lag.max=p,demean=TRUE)$acf
		},

		mle={

			if(p==1)
				autocor.vect<-c(1,phi.vect[1])

			else{

				autocor.vect<-rep(0,p-1) # initialize from lags 1 to p-1
				Phi.mat<-matrix(nrow=(p-1),ncol=(p-1))
					# Phi.mat%*%autocor.vect = -phi.vect

				Phi.mat[1,1]<-phi.vect[2]-1
				if(p>2)
					Phi.mat[1,-1]<-phi.vect[-c(1:2)]

				phi.vect<-c(phi.vect,rep(0,p)) # 0's added for use in loop only

				for(h in 2:p-1){

					Phi.mat[h,h]<-phi.vect[2*h]-1
					Phi.mat[h,h-1]<-phi.vect[1]+phi.vect[2*h-1]

					if(h>2)
						Phi.mat[h,1:(h-2)]<-phi.vect[(h-1):2]+phi.vect[(h+1):(2*h-2)]

					if(2*h<p)
						Phi.mat[h,(h+1):(p-h)]<-phi.vect[(2*h+1):p]
				
				}

				Phi.mat<-replace(Phi.mat,is.na(Phi.mat),0)

				phi.vect<-phi.vect[1:p] # revert to original phi.vect

				autocor.vect<-solve(Phi.mat,-phi.vect[-p]) 
					# Phi.mat%*%autocor.vect = -phi.vect

				autocor.vect<-c(1,autocor.vect) # add back lag 0

				autocor.vect<-c(autocor.vect,phi.vect%*%autocor.vect[p:1]) 
					# add back lag p
			}

		}
	)

	for(i in (p+2):n)
		autocor.vect<-c(autocor.vect,phi.vect%*%autocor.vect[(i-1):(i-p)])

	Sig<-diag(n)

	for(i in 1:(n-1))
		Sig[i,((i+1):n)]<-autocor.vect[2:(n-i+1)]

	return(list(phi=phi.vect,Sig=(Sig+t(Sig)-diag(n)),fit=ar.p.fit))

}

cable.ar.p.resid<-function(ar.p.fit){

	###############
	# stand-alone #
	##################################################
	# This function computes fitted residuals and    #
	# innovations for a 'cable.ar.p.iter()' object   #
	# or the '$cable' element of a 'bentcable.ar()'  #
	# object, where p>0.                             #
	##################################################

	b0<-ar.p.fit$est[1]
	b1<-ar.p.fit$est[2]
	b2<-ar.p.fit$est[3]
	tau<-ar.p.fit$est[4]

	if(ar.p.fit$stick){

		phi.vect<-ar.p.fit$est[-c(1:4)]
		gamm<-0
	}

	else {

		gamm<-ar.p.fit$est[5]
		phi.vect<-ar.p.fit$est[-c(1:5)]
	}

	y.vect<-ar.p.fit$y
	t.vect<-ar.p.fit$t
	n<-ar.p.fit$n
	pp<-ar.p.fit$p

	return(ar.p.resid.core(b0,b1,b2,tau,gamm,phi.vect,y.vect,t.vect,n,pp))
}

stick.ar.0<-function(init.vect,y.vect,t.vect=NULL,n=NA){

	if(is.na(n))
		n<-length(y.vect)

	if(is.null(t.vect))
		t.vect<-c(0:(n-1))

	message("Trying 'nls()' Gauss-Newton algorithm...")

	fit<-try(nls(formula=y.vect~fullcable.t(t.vect,b0,b1,b2,tau,0),
		start=list(b0=init.vect[1],b1=init.vect[2],b2=init.vect[3],
		tau=init.vect[4]),trace=TRUE))
	
   #if(length(fit)!=6){
   if(inherits(fit,"try-error")){

		message("Trying 'nls()' Golub-Pereyra algorithm...")

		fit<-try(nls(formula=y.vect~fullcable.t(t.vect,b0,b1,b2,tau,0),
			algorithm="plinear",
			start=list(b0=init.vect[1],b1=init.vect[2],b2=init.vect[3],
			tau=init.vect[4]),trace=TRUE))

      #if(length(fit)!=6){
		if(inherits(fit,"try-error")){

			message("Trying 'nls()' Port algorithm...")

			fit<-try(nls(formula=y.vect~fullcable.t(t.vect,b0,b1,b2,tau,0),
				algorithm="port",
				start=list(b0=init.vect[1],b1=init.vect[2],b2=init.vect[3],
				tau=init.vect[4]),trace=TRUE))

         #if(length(fit)!=6){
			if(inherits(fit,"try-error")){
				stop("Initial fit unsuccessful. Please try alternative initial values.")
				return()
			}

		}
	}

	message("Converged!")

	return(list(fit=fit,y=y.vect,t=t.vect,n=n,p=0,stick=TRUE))
}


#----------------------------------------------------------------
# find confidence interval for point at which slope changes sign
#
# =================================================
# CANNOT HANDLE TIME VECTOR THAT'S NOT c(0,1,2,...)
# =================================================
#----------------------------------------------------------------

alpha<-function(theta.vect,summand,what){

	# computes alpha value in gradient formulas

	# `summand' is time point t
	# `what' is index of alpha

	what<-as.character(what)

	tau<-theta.vect[4]

	if(length(theta.vect)>4)
		gamm<-theta.vect[5]
	else
		gamm<-0

	return(
		switch(what,
			"1"=(summand>tau+gamm),
			"2"={
				if(gamm==0)
					0
				else
					(abs(summand-tau)<=gamm)*(summand-tau+gamm)/(2*gamm)
			},
			"3"={
				if(gamm==0)
					0
				else
					(abs(summand-tau)<=gamm)*(1-((summand-tau)/gamm)^2)/4
			},
			"4"=summand-tau
		)
	)
}

alph.vect<-function(theta.vect,p,summand,position){

	# produces alpha vector in gradient formulas

	# `summand' is time point t
	# `position' is variable wrt which derivative is taken
	# `position' can only be 3,4,5

	position<-as.character(position)

	b2<-theta.vect[3]

	if(length(theta.vect)>4)
		gamm<-theta.vect[5]
	else
		gamm<-0

	a.vect<-NULL

	switch(position,
		"3"={
			for(i in 0:p){
				a.vect<-c(a.vect,
					alpha(theta.vect,summand-i,1)*
					alpha(theta.vect,summand-i,4)+
					gamm*alpha(theta.vect,summand-i,2)^2)
			}	
		},

		"4"={
			for(i in 0:p){
				a.vect<-c(a.vect,
					alpha(theta.vect,summand-i,1)+
					alpha(theta.vect,summand-i,2))
			}	
		},

		"5"={
			for(i in 0:p){
				a.vect<-c(a.vect,
					alpha(theta.vect,summand-i,3))
			}
		}
	)

	return(a.vect)
}

grad<-function(theta.vect,phi.vect,summand,position){

	# takes gradient of score wrt theta

	# `summand' is time point t
	# `position' is variable wrt which derivative is taken

	position<-as.character(position)

	p<-length(phi.vect)

	phi.vect<-c(-1,phi.vect)
	b2<-theta.vect[3]

	return(
		switch(position,
			"1"=-sum(phi.vect),
			"2"=-sum(phi.vect*(summand-c(0:p))),
			"3"=-sum(phi.vect*alph.vect(theta.vect,p,summand,3)),
			"4"=b2*sum(phi.vect*alph.vect(theta.vect,p,summand,4)),
			"5"=-b2*sum(phi.vect*alph.vect(theta.vect,p,summand,5))
		)
	)

}

find.fisher<-function(n,theta.vect,phi.vect,sigma.sq){

	fisher.dim<-length(theta.vect)

	fisher<-matrix(0,ncol=fisher.dim,nrow=fisher.dim) # initialize
	working.fisher<-fisher

	for(summand in 0:(n-1)){
		for(j in 1:fisher.dim){
			for(k in j:fisher.dim){
				working.fisher[j,k]<-grad(theta.vect,phi.vect,summand,j)*
					grad(theta.vect,phi.vect,summand,k)
			}
		}
		fisher<-fisher+working.fisher
	}

	for(j in 2:fisher.dim)
		for(k in 1:(j-1))
			fisher[j,k]<-fisher[k,j]

	return(fisher/sigma.sq)
}



cable.change.conf<-function(ar.p.fit,level){

	###############
	# stand-alone #
	#####################################################
	# This function computes an approximate confidence  #
	# interval for the CTP of an AR(p) bent cable.      #
	#                                                   #
	# 'ar.p.fit' is a 'cable.ar.p.iter()' object, or    #
	# the '$cable' element of a 'bentcable.ar()'        #
	# object.                                           #
	#                                                   #
	# `level' is between 0 and 1.                       #
	#                                                   #
	# WARNING:                                          #
	# This function can only handle cases where time    #
	# steps are 0, 1, 2, ...                            #
	#####################################################

	stick<-ar.p.fit$stick

	if(!stick)
		k<-5
	else{
		k<-4
		gamm<-0
	}

	n<-ar.p.fit$n
	p<-ar.p.fit$p

	if(p>0){
	
		theta.vect<-ar.p.fit$est[1:k]
		phi.vect<-ar.p.fit$est[-c(1:k)]

		if(k==5)
			gamm<-theta.vect[5]

		if(ar.p.fit$method=="css")
			sigma.sq<-ar.p.fit$ar.p.fit$val/(n-p)
		else{
#			sigma.sq<-ar.p.fit$ar.p.fit$var
			sigma.sq<-sse.ar.p.core(theta.vect[1],theta.vect[2],theta.vect[3],
				theta.vect[4],gamm,phi.vect,ar.p.fit$y,ar.p.fit$t,n,p)/(n-p)
		}
	}

	else{

		theta.vect<-coef(ar.p.fit$fit)[1:k]
		phi.vect<-NULL

		fit.summ<-summary(ar.p.fit$fit)
		sigma.sq<-fit.summ$sig^2*fit.summ$df[2]/n

	}

	fisher.mat<-find.fisher(n,theta.vect,phi.vect,sigma.sq)

	tau<-theta.vect[4]
	
	if(!stick){

		b1<-theta.vect[2]
		b2<-theta.vect[3]
		gamm<-theta.vect[5]

		change.hat<-tau-gamm-2*b1*gamm/b2

		alpha.vect<-c(0,-2*gamm/b2,2*b1*gamm/b2^2,1,-(2*b1+b2)/b2)

		v<-alpha.vect%*%solve(fisher.mat)%*%alpha.vect
	}

	else{

		v<-solve(fisher.mat)[4,4]

		if(is.na(v)|v==0)
			return()

		change.hat=tau

	}	

	z<-qnorm(level+(1-level)/2)

	return(list(change.hat=as.numeric(change.hat),var=v,
		interval=change.hat+c(-1,1)*z*sqrt(v)))

}

Try the bentcableAR package in your browser

Any scripts or data that you put into this service are public.

bentcableAR documentation built on May 30, 2017, 2:11 a.m.