version 1.0

DESCRIPTION

@@ -0,0 +1,13 @@
+Package: Kpart
+Type: Package
+Title: Spline Fitting
+Version: 1.0
+Date: 2012-08-02
+Author: Eric Golinko
+Depends: leaps
+Maintainer: <egolinko@gmail.com>
+Description: Kpart spline fitting
+License: GPL (>= 2)
+Packaged: 2012-08-06 12:20:32 UTC; tf2
+Repository: CRAN
+Date/Publication: 2012-08-06 18:42:00

MD5

@@ -0,0 +1,12 @@
+ae77ff834325c20eea4c679fd78a9b20 *DESCRIPTION
+682aa6143db5599c48a55c2c627a6f6d *NAMESPACE
+77f662926661ebeff90a3a1b413c7051 *R/Kpart.example.R
+54417a55a4b1fcf5648257776d504863 *R/Kpart.knots.R
+c11c45b66b5d93baaf3ba360327caba3 *R/Kpart.plot.R
+24632d59a6ce136f0cb6360d6e25dd29 *R/lm.Kpart.R
+9489e67ee4ef536779777b2ab4f2aecd *inst/extdata/example_data.txt
+7964959f7a09cc72b6bd015cba55342c *man/Kpart-package.Rd
+95f21b8d4f078bab3a39d96a1cc4cecf *man/Kpart.example.Rd
+2ddd746c19ad688a680126ca1c47bf8a *man/Kpart.knots.Rd
+2fd46e3f24d575946304fa11aa0de26a *man/Kpart.plot.Rd
+7b95d60a5b2862bd5d04d562695f8227 *man/lm.Kpart.Rd

NAMESPACE

@@ -0,0 +1,7 @@
+export(
+Kpart.knots,
+lm.Kpart,
+Kpart.plot,
+Kpart.example
+)
+exportPattern("^[[:alpha:]]+")

R/Kpart.example.R

@@ -0,0 +1,103 @@
+Kpart.example <-
+function(data=NULL,K=4)
+{
+
+	if(is.null(data)) {
+		p = paste(system.file("extdata", package="Kpart"), "/example_data.txt", sep="")
+		data = read.table(p, header=T)
+	}
+
+	#variable assignment
+
+	x=data[,2]
+	x.2=x^2
+	x.3=x^3
+	y=data[,1]
+	L=floor(length(y)/K)
+
+	#Min/Max algorithm to find absolute maximum deviate in the kth partition
+
+	m=matrix(nrow=K,ncol=1)
+	for (k in 1:K) m[k]=abs(max(y[(L*(k-1)+1):(k*L)])-mean(y[(L*(k-1)+1):(k*L)]))
+
+	W=matrix(nrow=L, ncol=K)
+	for (l in 1:L) for (k in 1:K) W[l,k]=ifelse(m[k] != abs(y[(L*(k-1)+l)]-mean(y[(L*(k-1)+1):(k*L)])), 0, x[(L*(k-1)+l)])
+	pops=colSums(W)
+
+	#matrix of potential spline knots
+
+	TT=matrix(nrow=length(x), ncol=K)
+	for (i in 1:length(x)) for (k in 1:K) TT[i,k]=(ifelse(0 > (x[i]-pops[k])^3 , 0 ,(x[i]-pops[k])^3))
+	bigT=data.frame(TT)
+	for(i in 1:length(pops))
+	colnames(bigT)[i]=paste("X",as.character(pops[i]),sep="")
+
+	#model selection process
+
+	X=cbind(x,x.2,x.3,bigT)
+	r=regsubsets(X,y,,method="exhaustive",nbest=1,nmax=(K+4))
+	a=(summary(r)[c("bic")])
+	AA=data.frame(a)
+	V=cbind(1:length(t(AA)),AA)
+	aa=matrix(nrow=length(t(AA)),ncol=1)
+	for (i in 1:length(t(AA))) aa[i,1]=ifelse(min(V[,2])==V[i,2],V[i,1],0)
+	aaa=data.frame(aa)
+	aaaa=colSums(aaa)
+	c=coef(r,id=aaaa)
+	cc=data.frame(c)
+	C=matrix(nrow=2,ncol=length(c))
+	C[2,]=cc[,1]
+	C[1,]=names(c)
+
+	#matrix of spline knots
+
+	t=matrix(nrow=length(C[1,]),ncol=length(names(bigT)))
+	for (i in 1:length(C[1,])) 
+	for (j in 1:length(names(bigT))) 
+	t[i,j]=(ifelse(C[1,i]==(names(bigT)[j]),(names(bigT)[j]),NA))
+	tt=na.omit(c(t))
+
+	P=matrix(nrow=K, ncol=1)
+	for (i in 1:length(tt)) P[i]=ifelse(tt[i]=="NA", NA, paste("bigT",tt[i], sep="$"))
+	PP=na.omit(as.character(P))
+
+	#matrix of x,x^2, and x^3 terms 
+
+	qq=matrix(nrow=3,ncol=3)
+	for (i in 2:4) qq[(i-1),1]=ifelse(C[1,i]=="x","x",NA)
+	for (i in 2:4) qq[(i-1),2]=ifelse(C[1,i]=="x.2","x.2",NA)
+	for (i in 2:4) qq[(i-1),3]=ifelse(C[1,i]=="x.3","x.3",NA)
+	QQ=na.omit(as.character(qq))
+	Q=c(QQ)
+
+	EE=c(Q,PP)
+	E=paste(EE,collapse="+")
+
+	#paste of terms to be estimated by lm()
+
+	FF=as.formula(paste("y~ ", E))
+	f=lm(FF)
+
+	#plot of x and y axis and fitted values with vertical lines at knots 
+
+	plot(x,y,xlab="your_x",ylab="your_y")
+	lines(x,f$fitted.values, col="red")
+
+	lineyears=matrix(nrow=length(PP),ncol=length(colSums(W)))
+	for(i in 1:length(PP))
+	for(j in 1:length(colSums(W)))
+	lineyears[i,j]=ifelse(data.frame(strsplit(PP[i],"X"))[2,1]==colSums(W)[j],colSums(W)[j],0)
+
+	spline.knots=colSums(lineyears)[which(colSums(lineyears)!=0)]
+
+	par(new=TRUE)
+	for (i in 1:length(spline.knots)) 
+	abline(v=spline.knots[i], lty=2)
+	#prints spline knots used in model
+
+	print(spline.knots)
+	#returns linear model which summary(lm(.)) can be used
+
+	return(f)
+
+}

R/Kpart.knots.R

@@ -0,0 +1,81 @@
+Kpart.knots <-
+function(data,K)
+{
+
+	#library(leaps)
+
+	#variable assignment
+	x=data[,2]
+	x.2=x^2
+	x.3=x^3
+	y=data[,1]
+	L=floor(length(y)/K)
+
+	#Min/Max algorithm to find absolute maximum deviate in the kth partition
+	m=matrix(nrow=K,ncol=1)
+	for (k in 1:K) m[k]=abs(max(y[(L*(k-1)+1):(k*L)])-mean(y[(L*(k-1)+1):(k*L)]))
+
+	W=matrix(nrow=L, ncol=K)
+	for (l in 1:L) for (k in 1:K) W[l,k]=ifelse(m[k] != abs(y[(L*(k-1)+l)]-mean(y[(L*(k-1)+1):(k*L)])), 0, x[(L*(k-1)+l)])
+	pops=colSums(W)
+
+	#matrix of potential spline knots
+	TT=matrix(nrow=length(x), ncol=K)
+	for (i in 1:length(x)) for (k in 1:K) TT[i,k]=(ifelse(0 > (x[i]-pops[k])^3 , 0 ,(x[i]-pops[k])^3))
+	bigT=data.frame(TT)
+	for(i in 1:length(pops))
+	colnames(bigT)[i]=paste("X",as.character(pops[i]),sep="")
+
+	#model selection process
+
+	X=cbind(x,x.2,x.3,bigT)
+	r=regsubsets(X,y,,method="exhaustive",nbest=1,nmax=(K+4))
+	a=(summary(r)[c("bic")])
+	AA=data.frame(a)
+	V=cbind(1:length(t(AA)),AA)
+	aa=matrix(nrow=length(t(AA)),ncol=1)
+	for (i in 1:length(t(AA))) aa[i,1]=ifelse(min(V[,2])==V[i,2],V[i,1],0)
+	aaa=data.frame(aa)
+	aaaa=colSums(aaa)
+	c=coef(r,id=aaaa)
+	cc=data.frame(c)
+	C=matrix(nrow=2,ncol=length(c))
+	C[2,]=cc[,1]
+	C[1,]=names(c)
+
+	#matrix of spline knots
+
+	t=matrix(nrow=length(C[1,]),ncol=length(names(bigT)))
+	for (i in 1:length(C[1,])) 
+	for (j in 1:length(names(bigT))) 
+	t[i,j]=(ifelse(C[1,i]==(names(bigT)[j]),(names(bigT)[j]),NA))
+	tt=na.omit(c(t))
+
+	P=matrix(nrow=K, ncol=1)
+	for (i in 1:length(tt)) P[i]=ifelse(tt[i]=="NA", NA, paste("bigT",tt[i], sep="$"))
+	PP=na.omit(as.character(P))
+
+	#matrix of x,x^2, and x^3 terms 
+
+	qq=matrix(nrow=3,ncol=3)
+	for (i in 2:4) qq[(i-1),1]=ifelse(C[1,i]=="x","x",NA)
+	for (i in 2:4) qq[(i-1),2]=ifelse(C[1,i]=="x.2","x.2",NA)
+	for (i in 2:4) qq[(i-1),3]=ifelse(C[1,i]=="x.3","x.3",NA)
+	QQ=na.omit(as.character(qq))
+	Q=c(QQ)
+
+	EE=c(Q,PP)
+	E=paste(EE,collapse="+")
+
+	#paste of terms to be estimated by lm()
+
+	FF=as.formula(paste("y~ ", E))
+	f=lm(FF)
+	#returns spline knots used in model
+	lineyears=matrix(nrow=length(PP),ncol=length(colSums(W)))
+	for(i in 1:length(PP))
+	for(j in 1:length(colSums(W)))
+	lineyears[i,j]=ifelse(data.frame(strsplit(PP[i],"X"))[2,1]==colSums(W)[j],colSums(W)[j],0)
+	spline.knots=colSums(lineyears)[which(colSums(lineyears)!=0)]
+	return(spline.knots)
+}

R/Kpart.plot.R

@@ -0,0 +1,93 @@
+Kpart.plot <-
+function(data,K)
+{
+
+	#library(leaps)
+
+	#variable assignment
+
+	x=data[,2]
+	x.2=x^2
+	x.3=x^3
+	y=data[,1]
+	L=floor(length(y)/K)
+
+	#Min/Max algorithm to find absolute maximum deviate in the kth partition
+
+	m=matrix(nrow=K,ncol=1)
+	for (k in 1:K) m[k]=abs(max(y[(L*(k-1)+1):(k*L)])-mean(y[(L*(k-1)+1):(k*L)]))
+
+	W=matrix(nrow=L, ncol=K)
+	for (l in 1:L) for (k in 1:K) W[l,k]=ifelse(m[k] != abs(y[(L*(k-1)+l)]-mean(y[(L*(k-1)+1):(k*L)])), 0, x[(L*(k-1)+l)])
+	pops=colSums(W)
+
+	#matrix of potential spline knots
+
+	TT=matrix(nrow=length(x), ncol=K)
+	for (i in 1:length(x)) for (k in 1:K) TT[i,k]=(ifelse(0 > (x[i]-pops[k])^3 , 0 ,(x[i]-pops[k])^3))
+	bigT=data.frame(TT)
+	for(i in 1:length(pops))
+	colnames(bigT)[i]=paste("X",as.character(pops[i]),sep="")
+
+	#model selection process
+
+	X=cbind(x,x.2,x.3,bigT)
+	r=regsubsets(X,y,,method="exhaustive",nbest=1,nmax=(K+4))
+	a=(summary(r)[c("bic")])
+	AA=data.frame(a)
+	V=cbind(1:length(t(AA)),AA)
+	aa=matrix(nrow=length(t(AA)),ncol=1)
+	for (i in 1:length(t(AA))) aa[i,1]=ifelse(min(V[,2])==V[i,2],V[i,1],0)
+	aaa=data.frame(aa)
+	aaaa=colSums(aaa)
+	c=coef(r,id=aaaa)
+	cc=data.frame(c)
+	C=matrix(nrow=2,ncol=length(c))
+	C[2,]=cc[,1]
+	C[1,]=names(c)
+
+	#matrix of spline knots
+
+	t=matrix(nrow=length(C[1,]),ncol=length(names(bigT)))
+	for (i in 1:length(C[1,])) 
+	for (j in 1:length(names(bigT))) 
+	t[i,j]=(ifelse(C[1,i]==(names(bigT)[j]),(names(bigT)[j]),NA))
+	tt=na.omit(c(t))
+
+	P=matrix(nrow=K, ncol=1)
+	for (i in 1:length(tt)) P[i]=ifelse(tt[i]=="NA", NA, paste("bigT",tt[i], sep="$"))
+	PP=na.omit(as.character(P))
+
+	#matrix of x,x^2, and x^3 terms 
+
+	qq=matrix(nrow=3,ncol=3)
+	for (i in 2:4) qq[(i-1),1]=ifelse(C[1,i]=="x","x",NA)
+	for (i in 2:4) qq[(i-1),2]=ifelse(C[1,i]=="x.2","x.2",NA)
+	for (i in 2:4) qq[(i-1),3]=ifelse(C[1,i]=="x.3","x.3",NA)
+	QQ=na.omit(as.character(qq))
+	Q=c(QQ)
+
+	EE=c(Q,PP)
+	E=paste(EE,collapse="+")
+
+	#paste of terms to be estimated by lm()
+
+	FF=as.formula(paste("y~ ", E))
+	f=lm(FF)
+
+	#plot of x and y axis and fitted values with vertical lines at knots 
+
+	plot(x,y,xlab="your_x",ylab="your_y")
+	lines(x,f$fitted.values, col="red")
+
+	lineyears=matrix(nrow=length(PP),ncol=length(colSums(W)))
+	for(i in 1:length(PP))
+	for(j in 1:length(colSums(W)))
+	lineyears[i,j]=ifelse(data.frame(strsplit(PP[i],"X"))[2,1]==colSums(W)[j],colSums(W)[j],0)
+
+	spline.knots=colSums(lineyears)[which(colSums(lineyears)!=0)]
+
+	par(new=TRUE)
+	for (i in 1:length(spline.knots)) 
+	abline(v=spline.knots[i], lty=2)
+}

R/lm.Kpart.R

@@ -0,0 +1,80 @@
+lm.Kpart <-
+function(data,K)
+{
+
+	#library(leaps)
+
+	#variable assignment
+
+	x=data[,2]
+	x.2=x^2
+	x.3=x^3
+	y=data[,1]
+	L=floor(length(y)/K)
+
+	#Min/Max algorithm to find absolute maximum deviate in the kth partition
+
+	m=matrix(nrow=K,ncol=1)
+	for (k in 1:K) m[k]=abs(max(y[(L*(k-1)+1):(k*L)])-mean(y[(L*(k-1)+1):(k*L)]))
+
+	W=matrix(nrow=L, ncol=K)
+	for (l in 1:L) for (k in 1:K) W[l,k]=ifelse(m[k] != abs(y[(L*(k-1)+l)]-mean(y[(L*(k-1)+1):(k*L)])), 0, x[(L*(k-1)+l)])
+	pops=colSums(W)
+
+	#matrix of potential spline knots
+
+	TT=matrix(nrow=length(x), ncol=K)
+	for (i in 1:length(x)) for (k in 1:K) TT[i,k]=(ifelse(0 > (x[i]-pops[k])^3 , 0 ,(x[i]-pops[k])^3))
+	bigT=data.frame(TT)
+	for(i in 1:length(pops))
+	colnames(bigT)[i]=paste("X",as.character(pops[i]),sep="")
+
+	#model selection process
+
+	X=cbind(x,x.2,x.3,bigT)
+	r=regsubsets(X,y,,method="exhaustive",nbest=1,nmax=(K+4))
+	a=(summary(r)[c("bic")])
+	AA=data.frame(a)
+	V=cbind(1:length(t(AA)),AA)
+	aa=matrix(nrow=length(t(AA)),ncol=1)
+	for (i in 1:length(t(AA))) aa[i,1]=ifelse(min(V[,2])==V[i,2],V[i,1],0)
+	aaa=data.frame(aa)
+	aaaa=colSums(aaa)
+	c=coef(r,id=aaaa)
+	cc=data.frame(c)
+	C=matrix(nrow=2,ncol=length(c))
+	C[2,]=cc[,1]
+	C[1,]=names(c)
+
+	#matrix of spline knots
+
+	t=matrix(nrow=length(C[1,]),ncol=length(names(bigT)))
+	for (i in 1:length(C[1,])) 
+	for (j in 1:length(names(bigT))) 
+	t[i,j]=(ifelse(C[1,i]==(names(bigT)[j]),(names(bigT)[j]),NA))
+	tt=na.omit(c(t))
+
+	P=matrix(nrow=K, ncol=1)
+	for (i in 1:length(tt)) P[i]=ifelse(tt[i]=="NA", NA, paste("bigT",tt[i], sep="$"))
+	PP=na.omit(as.character(P))
+
+	#matrix of x,x^2, and x^3 terms 
+
+	qq=matrix(nrow=3,ncol=3)
+	for (i in 2:4) qq[(i-1),1]=ifelse(C[1,i]=="x","x",NA)
+	for (i in 2:4) qq[(i-1),2]=ifelse(C[1,i]=="x.2","x.2",NA)
+	for (i in 2:4) qq[(i-1),3]=ifelse(C[1,i]=="x.3","x.3",NA)
+	QQ=na.omit(as.character(qq))
+	Q=c(QQ)
+
+	EE=c(Q,PP)
+	E=paste(EE,collapse="+")
+
+	#paste of terms to be estimated by lm()
+
+	FF=as.formula(paste("y~ ", E))
+	f=lm(FF)
+	#returns linear model which summary(lm(.)) can be used
+
+	return(f)
+}