`flrti` <- function(Y,X,sigma,deriv=2,weight=1,int=T,wrap=F,plot=F){ p <- ncol(X) beta.zero <- 0 extra.con <- NULL if (length(weight)>1) deriv <- (1:4)[weight[-1]>0] else{ weight <- c(weight,0,0,0,0) weight[deriv+1] <- 1} ### different weights when there is intercept ### if(int){ X <- cbind(1,X) weights<-c(0,rep(weight[1],p),0,rep(weight[1],p)) } else weights <- rep(weight[1],2*p) for (i in 2:5) if (weight[i]>0) weights <- c(weights,rep(weight[i],2*(p-i+1))) A <- makeA(p,deriv,int) colscale <- apply(X,2,function(x){sqrt(sum(x^2))}) newX <- t(t(X)/colscale) if (wrap) if (int) extra.con <- c(0,1/colscale[2],rep(0,p-2),-1/colscale[p+1],0,-1/colscale[2],rep(0,p-2),1/colscale[p+1],rep(0,2*nrow(A))) else extra.con <- c(1/colscale[1],rep(0,p-2),-1/colscale[p],-1/colscale[1],rep(0,p-2),1/colscale[p],rep(0,2*nrow(A))) lin.fit <- linearoptim(Y,newX,sigma,sqrt(2*p),weights,A%*%diag(1/colscale),extra.con) beta <- lin.fit$eta/colscale fit <- X%*%beta residuals <- Y-fit Rsq <- 1- sum(residuals^2)/sum((Y-mean(Y))^2) if (int){ beta.zero <- beta[1] beta <- beta[-1]} beta <- p*beta if (plot) plot(seq(0,1,length=p),beta,type='l',ylab="FLRTI Estimate", xlab="X axis") list(beta=beta,gamma=lin.fit$gamma,sol=lin.fit$sol, lin.fit=lin.fit$lpout,code=lin.fit$code,sigma=sigma,weight=weight[1], deriv=deriv,int=int,beta.zero=beta.zero,Rsq=Rsq,fitted.values=fit,residuals=residuals) } `flrti.boot` <- function(Y=NULL,X=NULL,obj=NULL,B=100,CI=.95,bootfit=NULL){ if (is.null(bootfit)){ n <- length(Y) sigma <- obj$sigma weight <- obj$weight int <- obj$int wrap <- obj$wrap deriv <- obj$deriv bootbeta <- matrix(0,B,ncol(X)) for (b in 1:B){ status <- 1 while (status!=0){ s <- sample(n,n,replace=T) bootY <- Y[s] bootX <- X[s,] bootfit <- flrti(bootY,bootX,sigma,deriv,weight,int,wrap) status <- bootfit$code} bootbeta[b,] <- bootfit$beta }} else{ obj <- bootfit$obj bootbeta <- bootfit$bootbeta} uci <- apply(bootbeta,2,function(x)quantile(x,1-(1-CI)/2)) lci <- apply(bootbeta,2,function(x)quantile(x,(1-CI)/2)) p <- length(uci) plot(seq(0,1,length=p),obj$beta,ylim=c(min(lci),max(uci)),type='l',lwd=3, xlab="X axis", ylab="Beta Confidence Intervals") lines(seq(0,1,length=p),lci,col=2) lines(seq(0,1,length=p),uci,col=2) list(obj=obj,bootbeta=bootbeta,uci=uci,lci=lci)} `flrti.cv` <- function(Y,X,sigma,cv=2,deriv=2,weight=1,s=NULL,int=T,wrap=F){ n <- length(Y) if (is.null(s)) s <- sample(n,replace=T) siglen <- length(sigma) wlen <- length(weight) error <- matrix(0,wlen,siglen) allerrors <- array(0,c(cv,wlen,siglen)) for (j in 1:siglen) for (k in 1:wlen){ for (i in 1:cv){ test.index <- s[(1+(i-1)*floor(n/cv)):(i*floor(n/cv))] if (i==cv) test.index <- s[(1+(i-1)*floor(n/cv)):n] testY <- Y[test.index] testX <- X[test.index,] trainY <- Y[-test.index] trainX <- X[-test.index,] trainfit <-flrti(trainY,trainX,sigma[j],deriv,weight[k],int=int,wrap=wrap) if (trainfit$code==0) allerrors[i,k,j]<-sum((testY-predict.flrti(trainfit,testX))^2)*n/length(testY) error[k,j] <- error[k,j]+allerrors[i,k,j]*length(testY)/n } code.errors <- sum(allerrors[,k,j]==0) if (code.errors