Nothing
## This file contains the functions used for analyzing causal heterogeneity
## using support vector machines under a double LASSO constraint.
## The functions consist of
## scale.func: demeans and standardizes variables.
## scale.func.2: gives variables a standard deviation of one
## FindIt: a wrapper function for SVM.func
## maketwoway<-rescales covariates and makes all twoway interactions.
## SVM.func: the workhorse function. Given values of y, control, treatment, and lambdas,
## it returns betas, the loss (GCV), fitted values, and size of the margin.
## find.lambda: Conducts the alternating line search.
## A large value of lambda.c is selected (10), and alternating
## line searches of decreasing radius are conducted.
## lars: An adaptation of the LARS algorithm of Efron, et al. 2004.
## Due to the large numbers being considered in the fitting,
## the code is adapted to not crash. The original code is commented out.
## The remainder of the code is untouched.
##################################################################
##################################################################
##################################################################
scale.func<-function(x) {
x2<-x-mean(x)
x2/sd(x2)}
scale.func.2<-function(x) {x2<-x
x2*sd(x2)}
#####################################################
##change here.
#' FindIt for Estimating Heterogeneous Treatment Effects
#'
#' \code{FindIt} returns a model with the most predictive treatment-treatment
#' interactions or treatment-covariate interactions.
#'
#' Implements the alternating line search algorithm for estimating the tuning
#' parameters, as described in Imai and Ratkovic (2013).
#'
#' @param model.treat A formula that specifies outcome and treatment variables.
#' @param model.main An optional formula that specifies pre-treatment
#' covariates to be adjusted.
#' @param model.int A formula specifying pre-treatment covariates to be
#' interacted with treatment assignments when \code{treat.type="single"}.
#' @param data An optional data frame, list or environment (or object coercible
#' by 'as.data.frame' to a data frame) containing the variables in the model.
#' If not found in 'data', the variables are taken from 'environment(formula)',
#' typically the environment from which 'FindIt' is called.
#' @param type "binary" for a binary outcome variable, which needs to be
#' \code{integer} class; "continuous" for a continuous outcome variable.
#' @param treat.type "single" for interactions between a single treatment
#' variable, which needs to be \code{integer} class, and multiple pre-treatment
#' covariates specified with \code{model.int}; "multiple" is used when
#' treatment-treatment interactions are of interest and \code{treat} is a
#' matrix of multiple treatments.
#' @param nway An argument passed to \code{makeallway} when
#' \code{treat.type="multiple"}. \code{FindIt} generates treatment-treatment
#' interactions up to the order specified with this argument. In general, it is
#' recommended to use the number of factorial treatments. The current version
#' covers up to four way interactions.
#' @param search.lambdas Whether to search for the tuning parameters for the
#' LASSO constraints. If \code{FALSE}, \code{lambdas} must be supplied.
#' @param lambdas Tuning parameters to be given to \code{FindIt}; only used if
#' \code{search.lambdas=FALSE}.
#' @param make.twoway If \code{make.twoway=TRUE}, all possible two-way
#' interactions for the pre-treatment covariates specified in \code{model.main}
#' and \code{model.int} are generated within \code{FindIt}. The default is set
#' to be \code{TRUE}.
#' @param make.allway If \code{make.allway=TRUE}, all possible
#' treatment-treatment interactions for multiple treatments are generated when
#' \code{treat.type="multiple"}. Interactions of the order up to the value of
#' \code{nway} is computed.
#' @param wts An optional set of scaling weights. The default is 1.
#' @param scale.c A set of weights for recaling the pre-treatment covariates;
#' only used if \code{make.twoway=FALSE}. \code{maketwoway} is useful for
#' generating these.
#' @param scale.int A set of weights for recaling the covariates to be
#' interacted with treatment variables ; only used if \code{make.twoway=FALSE}.
#' \code{maketwoway} is useful for generating these.
#' @param fit.glmnet Whether to fit using the coordinate descent method in
#' glmnet (TRUE) or the regularization path method of LARS (FALSE).
#' @param make.reference Whether to make a reference matrix to check which
#' columns are dropped when \code{makeallway=TRUE}.
#' @param reference.main If \code{make.allway=FALSE} and researchers generate a
#' matrix of all possible interactions between factorial treatments,
#' \code{reference} from \code{makeallway} function is better to be passed to
#' \code{FindIt} through this argument.
#' @param threshold An argument passed to \code{makeallway} when
#' \code{treat.type="multiple"}. Threshold to drop correlated columns when
#' \code{makeallway} is used.
#' @return \item{coefs}{A named vector of scaled coefficients}
#' \item{coefs.orig}{A vector of coefficients on the original scale, if scale.c
#' and scale.t was used} \item{fit}{Fitted values on an SVM scale}
#' \item{names.out}{Names of the coefficients} \item{y}{A vector of observed
#' outcomes} \item{X.c}{A matrix of pre-treatment covariates to be adjusted}
#' \item{X.t}{A matrix of treatments and treatment-treatment interactions, or
#' treatment-covariate interactions} \item{GCV}{GCV statistic at the minimum}
#' \item{ATE}{When \code{treat.type="single"}, the estimated ATE. When
#' \code{treat.type="multiple"}, the estimated treatment effect of each unique
#' treatment combination} \item{lambdas}{Tuning parameters used for the fit}
#' \item{reference}{When \code{treat.type="multiple"}, after making all
#' interaction terms, columns with no variation or columns perfectly correlated
#' with one of other columns are automatically dropped. \code{reference} shows
#' which columns are kept and dropped.}
#' @import lars
#' @import glmnet
#' @author Naoki Egami, Marc Ratkovic and Kosuke Imai.
#' @references Imai, Kosuke and Marc Ratkovic. 2013. ``Estimating Treatment
#' Effect Heterogeneity in Randomized Program Evaluation.'' Annals of Applied
#' Statistics, Vol.7, No.1(March), pp. 443-470.
#' \url{http://imai.fas.harvard.edu/research/files/svm.pdf}
#'
#' Egami, Naoki and Kosuke Imai. 2019. Causal Interaction in
#' Factorial Experiments: Application to Conjoint Analysis, Journal of the American Statistical Association.
#' \url{http://imai.fas.harvard.edu/research/files/int.pdf}
#' @examples
#'
#'
#' ###################################################
#' ## Example 1: Treatment-Covariate Interaction
#' ###################################################
#' data(LaLonde)
#'
#' ## The model includes a treatment variable,
#' ## nine covariates to be interacted with the treatment variable,
#' ## and the same nine covariates to be adjusted.
#'
#' \dontrun{
#'
#' ## Run to find the LASSO parameters
#' F1 <-FindIt(model.treat= outcome ~ treat,
#' model.main= ~ age+educ+black+hisp+white+
#' marr+nodegr+log.re75+u75,
#' model.int= ~ age+educ+black+hisp+white+
#' marr+nodegr+log.re75+u75,
#' data = LaLonde,
#' type="binary",
#' treat.type="single")
#' }
#'
#' ## Fit with uncovered lambda parameters.
#' F1 <-FindIt(model.treat= outcome ~ treat,
#' model.main= ~ age+educ+black+hisp+white+
#' marr+nodegr+log.re75+u75,
#' model.int= ~ age+educ+black+hisp+white+
#' marr+nodegr+log.re75+u75,
#' data = LaLonde,
#' type="binary",
#' treat.type="single",
#' search.lambdas=FALSE,
#' lambdas = c(-3.8760,-4.0025) )
#'
#' summary(F1)
#'
#' ## Returns all the estimated treatment effects.
#' pred1 <- predict(F1)
#' ## Top10
#' head(pred1$data, n=10)
#' ## Bottom 10
#' tail(pred1$data ,n=10)
#'
#' ## Visualize all the estimated treatment effects.
#' \dontrun{
#' plot(pred1)
#' }
#'
#' ###################################################
#' ## Example 2: Treatment-Treatment Interaction
#' ###################################################
#'
#' \dontrun{
#' data(GerberGreen)
#'
#' ## The model includes four factorial treatments and
#' ## all two, three, four-way interactions between them.
#' ## Four pre-treatment covariates are adjusted.
#'
#' ## Run to search for lambdas.
#' F2<- FindIt(model.treat= voted98 ~ persngrp+phnscrpt+mailings+appeal,
#' nway=4,
#' model.main= ~ age+majorpty+vote96.1+vote96.0,
#' data = GerberGreen,
#' type="binary",
#' treat.type="multiple")
#'
#' ## Fit, given selected lambdas.
#' F2<- FindIt(model.treat= voted98 ~ persngrp+phnscrpt+mailings+appeal,
#' nway=4,
#' model.main= ~ age+majorpty+vote96.1+vote96.0,
#' data = GerberGreen,
#' type="binary",
#' treat.type="multiple",
#' search.lambdas=FALSE,
#' lambdas=c(-15.000,-6.237))
#'
#' ## Returns coefficient estimates.
#' summary(F2)
#'
#' ## Returns predicted values for unique treatment combinations.
#' pred2 <- predict(F2,unique=TRUE)
#' ## Top 10
#' head(pred2$data, n=10)
#' ## Bottom 10
#' tail(pred2$data, n=10)
#'
#' ## Visualize predicted values for each treatment combination.
#' plot(pred2)
#' }
#'
#' @export
FindIt <- function(model.treat, model.main,model.int,data=NULL,
type="binary",
treat.type="multiple",
nway,
search.lambdas=TRUE,lambdas=NULL,
make.twoway=TRUE,make.allway=TRUE,
wts=1,scale.c=1,scale.int=1,
fit.glmnet=TRUE,
make.reference=TRUE,
reference.main=NULL,
threshold=0.999999){
## if(make.allway==FALSE & is.null(reference.main)){
## warning("Need to specify the corresponding reference of treatment matrix")
## break
## }
unique <- FALSE
## Check nway
if(missing(nway) & treat.type=="multiple") {
stop("Need to specify nway")
}
if(missing(nway) & treat.type=="single") {
nway <- 2
}
if(sum(wts<0)>0) {
stop("Only non-negative weights allowed")
}
## Check whether the model has the main effect covariates or not.
if(!missing(model.main)|!missing(model.int)){
main <- TRUE
}
if(missing(model.int) & !missing(model.main)){
model.int <- model.main
main <- TRUE
}
if(!missing(model.int) & missing(model.main)){
model.main <- model.int
main <- TRUE
}
if(missing(model.main) & missing(model.int)){
main <- FALSE
}
## Extract the data
treat.frame <- model.frame(model.treat,data, na.action=NULL)
if(main){
main.frame <- model.frame(model.main, data, na.action=NULL)
int.frame <- model.frame(model.int, data, na.action=NULL)
}
## NA : List-wise deletion.
data.which <- c()
if(main){
for(i in 1:nrow(treat.frame)){
if(any(is.na(main.frame[i,])) | any(is.na(int.frame[i,])) |
any(is.na(treat.frame[i,])) )
{data.which[i] <- 0}
else{data.which[i] <- 1}
}
}
if(main==FALSE){
for(i in 1:nrow(treat.frame)){
if(any(is.na(treat.frame[i,])))
{data.which[i] <- 0}
else{data.which[i] <- 1}
}
}
treat.frame <- treat.frame[data.which==1,]
treat.orig <- treat.frame[,-1]
treat <- treat.frame[,-1]
treat <- as.matrix(treat)
y <- y.orig <- treat.frame[,1]
if(type=="binary"){
y.orig <- y
y <-(2*y-1)
}
if(main){
X.c <- X.c.orig <- main.frame[data.which==1,]
X.int <- X.int.orig <- int.frame[data.which==1,]
X.c <- X.c.orig <- as.data.frame(X.c)
X.int <- X.int.orig <- as.data.frame(X.int)
if(ncol(X.c)>1){
if(any(apply(X.c,2,FUN=function(x) length(unique(x)) == 1))){
warning("There is a column with no variation
in the main effect matrix.")
}
}
if(ncol(X.int)>1){
if(any(apply(X.int,2,FUN=function(x) length(unique(x)) == 1))){
warning("There is a column with no variation
in the interaction matrix.")
}
}
if(treat.type=="single" & main==FALSE){
warning("Need to specify the covariates to interact with treatment.")
}
}
## Make multiple-treatment matrix
if(treat.type == "multiple"){
if(make.allway==TRUE){
Allway <- makeallway(treat,threshold,make.reference=make.reference,nway=nway)
X.t <- Allway$FinalData
reference.main <- Allway$reference
}
if(make.allway==FALSE){
X.t <- as.matrix(treat)
## reference.main <- reference.main
reference.main <- "No reference"
}
}
## print("X.t")
## print(head(X.t))
## print(head(X.c))
## print(head(X.int))
## Maketwoway
if(main){
X.c.old <- X.c
X.int.old <- X.int
if(make.twoway==TRUE){
if(ncol(X.c.old)>1){
XC <- maketwoway(X.c,center=TRUE,wts=wts)
X.c <- XC$X
scale.c <- XC$scale.X
XC.deleted <- XC$D
frame.meanC <- XC$frame.mean
}
if(ncol(X.c.old)==1){
form <- ~ .
X.c.old <- as.data.frame(X.c.old)
X.c1 <- model.matrix(form,data=X.c.old)
X.c1 <- as.data.frame(X.c1)
X.c <- X.c1[,-1]
if(is.vector(X.c)){
scale.c <- sd(X.c)
frame.meanC <- mean(X.c)
X.c <- (X.c-mean(wts^.5*X.c))/sd(X.c)
X.c <- as.data.frame(X.c)
colnames(X.c) <- colnames(X.c1)[-1]
}else{
scale.c <- apply(X.c,2,sd)
frame.meanC <- apply(X.c,2,function(x) mean(wts^.5*x))
X.c <- apply(X.c,2,FUN=function(x) (x-mean(wts^.5*x))/sd(x))
}
}
if(ncol(X.int.old)>1){
XInt <- maketwoway(X.int,center=TRUE,wts=wts)
X.int <- XInt$X
scale.int <- XInt$scale.X
XInt.deleted <- XInt$D
frame.meanInt <- XInt$frame.mean
}
if(ncol(X.int.old)==1){
form <- ~ .
X.int.old<- as.data.frame(X.int.old)
X.int1 <- model.matrix(form,data=X.int.old)
X.int1 <- as.data.frame(X.int1)
X.int <- X.int1[,-1]
if(is.vector(X.int)){
scale.int <- sd(X.int)
frame.meanInt <- mean(X.int)
X.int <- (X.int-mean(wts^.5*X.int))/sd(X.int)
X.int <- as.data.frame(X.int)
colnames(X.int) <- colnames(X.int1)[-1]
}
else{
scale.int <- apply(X.int,2,sd)
frame.meanInt <- apply(X.int,2,function(x) mean(wts^.5*x))
X.int <- apply(X.int,2,FUN=function(x) (x-mean(wts^.5*x))/sd(x))
}
}
}
## Remove all elements with no variance and add an intercept.
if(is.matrix(X.c)){
scale.c <- scale.c[apply(X.c,2,sd)>0]
X.c<-X.c.1 <- X.c[,apply(X.c,2,sd)>0]
## X.c<-unique(X.c,MARGIN=2)
scale.c <- scale.c[colnames(X.c.1)==colnames(X.c)]
frame.meanC <- frame.meanC[colnames(X.c.1)==colnames(X.c)]
}
if(is.vector(X.c)){
X.c<-cbind(1,X.c)
colnames(X.c)[1]<-"Intercept"
}
else{
if(sd(X.c[,1])>0){
X.c<-cbind(1,X.c)
colnames(X.c)[1]<-"Intercept"
}
}
name.c <- colnames(X.c)
## This X.c is used in the SVM.
## Remove all elements with no variance and add an intercept.
if(is.matrix(X.int)){
scale.int <- scale.int[apply(X.int,2,sd)>0]
frame.meanInt <- frame.meanInt[apply(X.int,2,sd)>0]
X.int<- X.int.1<- X.int[,apply(X.int,2,sd)>0]
scale.int <- scale.int[colnames(X.int.1)==colnames(X.int)]
frame.meanInt <- frame.meanInt[colnames(X.int.1)==colnames(X.int)]
}
if(is.vector(X.int)){
X.int<-cbind(1,X.int)
colnames(X.int)[1]<-"Intercept"
}
else{
if(sd(X.int[,1])>0){
X.int<-cbind(1,X.int)
colnames(X.int)[1]<-"Intercept"
}
}
name.int <- colnames(X.int)
## This X.int is used in the SVM
## Make X.t
## I changed X.c to X.int
if(treat.type == "single") {
treat.temp <- as.numeric(treat>0)
X.t<-cbind(treat*1e4,treat.temp*X.int[,-1])
X.t[treat!=0,-1]<-apply(X.t[treat!=0,-1],2, FUN=function(x) x-mean(x))
## This X.t is used.
colnames(X.t)<-c("treat",paste("treat",colnames(X.int)[-1],sep=":"))
}
}
###################################
## Create the FindIt.fit function.
###################################
FindIt.fit<-function(y,X.c,X.t,X.int,treat,treat.type,
search.lambdas,lambdas,
wts,scale.c,scale.int,
reference.main,
fit.glmnet,
type,
X.c.orig,X.int.orig,treat.orig,
threshold,
unique,
main){
## Give wts mean 1.
wts<-wts/mean(wts)
y <- y*wts^.5
X.t <- X.t*wts^.5
name.t <- colnames(X.t)
if(main){
X.c <- X.c*wts^.5
## X.int <- X.int*wts^.5
name.c <- colnames(X.c)
name.int <- colnames(X.int)
}
## treat<-treat*wts^.5
X.t <- as.matrix(X.t)
if(main){
X.c <- as.matrix(X.c)
X.int <- as.matrix(X.int)
}
if(search.lambdas==TRUE){
if(main){
lambdas <-search.lambda(y,X.c,X.t,
fit.glmnet=fit.glmnet,
main=main,
type=type)
}
if(main==FALSE){
lambdas <- search.lambda.nocov(y=y,X.t=X.t,
fit.glmnet=fit.glmnet,
main=main,
type=type)
}
}
if(main){
A<-SVM.func(y,X.c[,-1],X.t,
lambdas[1],lambdas[2],
fit.glmnet=fit.glmnet,
type=type,
main=main)
}
if(main==FALSE){
A<-SVM.func(y=y,X.c[,-1],X.t=X.t,
lambda.c=1,lambda.t=lambdas,
fit.glmnet=fit.glmnet,
type=type,
main=main)
}
## When we input X.c doesnt have intercept.
if(treat.type=="single"){
## Keep the names of coefficients.
reference.coef <- as.data.frame(A$beta[-c(1:(dim(X.c)[2]))])
rownames(reference.coef) <- c("treat",colnames(X.int[,-1]))
if(type=="binary"){
reference.coef <- reference.coef/2
}
## Generate outcome treatment and no treatment.
## I changed the code.
preds.treat<-cbind(X.c,cbind(1e4,X.int[,-1]))%*%A$beta
preds.control<-cbind(X.c,0*cbind(1e4,X.int[,-1]))%*%A$beta
## preds <- cbind(X.c,X.t) %*% A$beta
## diff <- mean(y.inter)- mean(preds)
## cbind(1e4,X.int[,-1]) is X.t when treat =1.
## preds.treat<-cbind(X.c,cbind(1e4,X.c[,-1]))%*%A$beta
## in here, X.c includes the intercept.
## 1e4 means treated group
## preds.control<-cbind(X.c,0*cbind(1e4,X.c[,-1]))%*%A$beta
if(type=="binary"){
preds.treat <- sign(preds.treat)*pmin(abs(preds.treat),1)
preds.control <- sign(preds.control)*pmin(abs(preds.control),1)
ATE <- mean(preds.treat-preds.control)/2
preds.diff<- (preds.treat-preds.control)/2
}
if(type=="continuous"){
preds.diff<-preds.treat-preds.control
ATE<-mean(preds.treat-preds.control)
}
scale.out<-c(1,scale.c,1e-4,scale.int)
Treatment.version <- NULL
}
if(treat.type == "multiple"){
## make the reference matrix.
if(main){
## preds <- cbind(X.c,X.t) %*% A$beta
## diff <- mean(y.inter) - mean(preds)
scale.out<-c(1,scale.c,rep(1,dim(X.t)[2]))
coef.orig.ref <- A$beta[-c(1:ncol(X.c))]/scale.out[-c(1:ncol(X.c))]
}
if(main==FALSE){
## preds <- cbind(1,X.t) %*% A$beta
## diff <- mean(y.inter) - mean(preds)
scale.out <- c(1,rep(1,ncol(X.t)))
## I need one more for intercept.
coef.orig.ref <- A$beta[-1]/scale.out[-1]
}
if(type=="binary"){
reference.coef <- as.data.frame(signif(coef.orig.ref/2,6))
}
if(type=="continuous"){
reference.coef <- as.data.frame(signif(coef.orig.ref,6))
}
rownames(reference.coef) <- colnames(X.t)
###############################
## create the Treatment.version
###############################
coefs <- A$beta
if(type=="binary"){
if(main){
preds.treat <- cbind(X.c,X.t) %*% coefs
preds.control <- X.c %*% coefs[c(1:ncol(X.c))]
preds.treat <- sign(preds.treat)*pmin(abs(preds.treat),1)
preds.control <- sign(preds.control)*pmin(abs(preds.control),1)
preds.diff <- (preds.treat-preds.control)/2
}
if(main==FALSE){
preds.treat <- X.t %*% coefs[-1]
preds.treat <- sign(preds.treat)*pmin(abs(preds.treat),1)
preds.diff <- preds.treat/2
}
ATE <- mean(preds.diff)
}
if(type=="continuous"){
if(main){
preds.diff <- X.t %*% coefs[-c(1:(dim(X.c)[2]))]
}
if(main==FALSE){
preds.diff <- X.t %*% coefs[-1]
}
ATE <- mean(preds.diff)
}
if(unique==FALSE){
##preds.control <- X.c %*% coefs[1:ncol(X.c)]
##preds <- cbind(X.c,X.t)%*%coefs
##preds.diff <- preds - preds.control
preds.diff <- preds.diff
Treatment.version <- cbind(preds.diff,treat.orig)
colnames(Treatment.version) <- c("Treatment.effect",
colnames(treat.orig))
Treatment.version <- as.data.frame(Treatment.version)
}
if(unique==TRUE){
X.t.u <- X.t
preds.diff.u <- as.data.frame(preds.diff)
rownames(preds.diff.u) <- rownames(X.t.u) <- seq(1:nrow(X.t.u))
treat.unique <- unique(treat.orig,MARGIN=1)
X.t.unique <- unique(X.t.u, MARGIN=1)
X.t.unique2 <- unique(X.t, MARGIN=1)
preds.diff <- preds.diff.u[rownames(preds.diff.u) %in%
rownames(X.t.unique),]
preds.diff <- as.data.frame(preds.diff)
##name back
rownames(X.t.unique) <- rownames(X.t.unique2)
rownames(preds.diff) <- rownames(X.t.unique2)
Treatment.version <- cbind(preds.diff,treat.unique)
colnames(Treatment.version) <- c("Treatment.effect",
colnames(treat.orig))
Treatment.version <- as.data.frame(Treatment.version)
Treatment.version<-
Treatment.version[order(Treatment.version$Treatment.effect,
decreasing=TRUE),]
}
}
if(main){
names.out <- c(colnames(X.c),colnames(X.t))
}
else{
names.out <- c("Intercept",colnames(X.t))
}
if(main){
output<-(list("coefs"=A$beta,"coefs.orig"=A$beta/scale.out,
"fit"=A$fit,
"scale.out"=scale.out,
"scale.c" = scale.c,
"scale.int"=scale.int,
"frame.meanC" = frame.meanC,
"frame.meanInt" = frame.meanInt,
"names.out"=names.out,
"y"=y,"X.c"=X.c,"X.t"=X.t,"X.int"=X.int,
"GCV"=exp(A$loss),
"name.c" =name.c,
"name.int" = name.int,
"name.t" =name.t,
"ATE" = ATE,
"Treatment.version" = Treatment.version,
"lambdas"=lambdas,
"reference.main"=reference.main,
"reference.coef"=reference.coef)
)
}
if(main==FALSE){
output<-(list("coefs"=coefs,"coefs.orig"=A$beta/scale.out,
"fit"=A$fit,
"scale.out"=scale.out,
"names.out"=names.out,
"y"=y,"X.t"=X.t,
"GCV"=exp(A$loss),
"name.t" =name.t,
"ATE" = ATE,
"Treatment.version" = Treatment.version,
"lambdas"=lambdas,
"reference.main"=reference.main,
"reference.coef"=reference.coef)
)
}
class(output)<-c("FindIt","list")
output
}
###########################################
## The FindIt.fit function is finished
##########################################
if(main==FALSE){
X.c <- X.int <- scale.c <- scale.int <- frame.meanC <- frame.meanInt <- NULL
X.c.orig <- X.int.orig <- NULL
}
Fit.output <- FindIt.fit(y=y,X.c=X.c,X.t=X.t,X.int=X.int,treat=treat,
treat.type=treat.type,
search.lambdas = search.lambdas,
lambdas = lambdas,
wts = wts,
scale.c = scale.c,
scale.int = scale.int,
reference.main = reference.main,
fit.glmnet = fit.glmnet,
type = type,
X.c.orig = X.c.orig,
X.int.orig = X.int.orig,
treat.orig = treat.orig,
threshold = threshold,
unique = unique,
main=main)
#########################################################
#### Make the reference matrix when treat.type is multiple.
#########################################################
if(make.reference==TRUE & make.allway==TRUE){
if(treat.type=="multiple"){
reference.main <- Fit.output$reference.main
reference.coef <- Fit.output$reference.coef
coef <- c()
for(i in 1:nrow(reference.main)){
if(rownames(reference.main)[i] %in% rownames(reference.coef)){
coef[i] <-
reference.coef[rownames(reference.coef)==
rownames(reference.main)[i],1]
}
else{coef[i] <- NA
}
}
reference.Final <- as.data.frame(cbind(reference.main,coef))
colnames(reference.Final) <- c(colnames(reference.main),
"Coefficients"
)
}
}else{
reference.Final <- "No reference"
}
if(treat.type=="single"){
reference.Final <- "No reference"
}
if(main){
output<-(list("coefs"=Fit.output$coefs,"coefs.orig"=Fit.output$coefs.orig,
"fit"=Fit.output$fit,
"scale.out"=Fit.output$scale.out,
"scale.c" = Fit.output$scale.c,
"scale.int" = Fit.output$scale.int,
"frame.meanC" = frame.meanC,
"frame.meanInt" = frame.meanInt,
"names.out"=Fit.output$names.out,
"y"=Fit.output$y,"X.c"=Fit.output$X.c,
"X.t"=Fit.output$X.t, "X.int"=Fit.output$X.int,
"treat" = Fit.output$treat,
"GCV"=Fit.output$GCV,
"y.orig"=y.orig, "X.c.orig" = X.c.orig,
"X.int.orig" = X.int.orig, "treat.orig" = treat.orig,
"name.c" = Fit.output$name.c,
"name.int" = Fit.output$name.int,
"name.t" = Fit.output$name.t,
"treat.type" = treat.type,
"model.main"=model.main,
"model.int" = model.int,
"model.treat"=model.treat,
"ATE"=Fit.output$ATE,
"data"=data,
"Treatment.version" = Fit.output$Treatment.version,
"lambdas"=Fit.output$lambdas,
"threshold" = threshold,
"make.twoway" = make.twoway,
"make.allway" = make.allway,
"model.main"=model.main,
"model.int"=model.int, "model.treat"=model.treat,
"reference"= reference.Final,"type" = type,
"main" =main,
"nway" =nway
))
}
if(main==FALSE){
output<-(list("coefs"=Fit.output$coefs,"coefs.orig"=Fit.output$coefs.orig,
"fit"=Fit.output$fit,
"scale.out"=Fit.output$scale.out,
"names.out"=Fit.output$names.out,
"y"=Fit.output$y,
"X.t"=Fit.output$X.t,
"treat" = Fit.output$treat,
"GCV"=Fit.output$GCV,
"y.orig"=y.orig,
"treat.orig" = treat.orig,
"name.t" = Fit.output$name.t,
"treat.type" = treat.type,
"model.treat"=model.treat,
"ATE"=Fit.output$ATE,
"data"=data,
"Treatment.version" = Fit.output$Treatment.version,
"lambdas"=Fit.output$lambdas,
"threshold" = threshold,
"make.twoway" = make.twoway,
"make.allway" = make.allway,
"reference"= reference.Final,
"type" = type,
"main" =main,
"nway" =nway
))
}
class(output)<-c("FindIt","list")
output
}
#####################################################
SVM.func<-function(y,X.c,X.t,lambda.c,lambda.t,wts=1,fit.glmnet,type,main){
### From the construction, this X.c doesnt include the intercept.
n<-length(y)
if(main){
X<-X2<-cbind(X.c,X.t)
}
else{
X <- X2 <- X.t
}
##Originally written to handle a string of lambda.c and lambda.t
##loss.run creates a matrix of losses, across values of lambda.
## loss.run<-matrix(0,nrow=length(lambda.c),ncol=length(lambda.t))
## min.params<-c(lambda.c[1],lambda.t[1])
## for(i in 1:length(lambda.c)){
## for(j in 1:length(lambda.t)){
##Scale the X matrix by lambdas.
X3<-X2
if(main){
if(is.matrix(X.c)){
X3[,1:dim(X.c)[2]]<-1/exp(lambda.c)*X2[,1:dim(X.c)[2]]
X3[,-c(1:dim(X.c)[2])]<-1/exp(lambda.t)*X2[,-c(1:dim(X.c)[2])]
}
if(is.vector(X.c)){
X3[,1]<-1/exp(lambda.c)*X2[,1]
X3[,-1]<-1/exp(lambda.t)*X2[,-1]
}
}
if(main==FALSE){
X3 <- 1/exp(lambda.t)*X2
}
##X3.c <- 1/exp(lambda.c[i])*X.c
##X3.t <- 1/exp(lambda.t[j])*X.t
##X3 <- cbind(X.c,X.t)
##Declare original values.
X.new<-X3
which.use<-rep(T,n)
beta.new<-beta.curr<-rep(0,dim(X2)[2]+1)
## this +1 is for the intercept.
## The loop to fit the LASSO.
## Center on the margin, fit the lasso, update coefficients,
##then drop all outside the margin (y*yhat>1)
if(type=="binary"){
for(beta.run in 1:100){
X.new.2<-apply(X.new[which.use,],MARGIN=2,FUN=function(x) x-mean(x))
X.new.2<-apply(X.new.2,MARGIN=2,
FUN=function(x)
if(length(unique(x))==1) rnorm(sum(which.use))
else x)
if(fit.glmnet==TRUE){
glmnet.1<-glmnet(X.new.2,y[which.use]-mean(y[which.use]),
family="gaussian",
lambda=c(5,4,3,2,seq(1.5,1,-.1)),
standardize=F)
beta.new[-1]<-as.vector(glmnet.1$beta[,10])
}else{
##This was the old L1 optimizer.
lasso1<-lars(X.new.2,(y[which.use]-mean(y[which.use])),
max.steps=dim(X)[2]+2,
normalize=F,
type="lasso",
eps=0)
beta.new[-1]<-as.vector(predict(lasso1,s=1,type="coef",
mode="lambda")$coef)
if(log(sum(which.use))/2*sum(beta.new[-1]!=0)>.9*sum(which.use)){
beta.new[-1]<-
as.vector(predict(lasso1,
s=min(c(floor(.9*sum(which.use)*2/
log(n)),dim(X.new)[2]*.8)),
type="coef",mode="step")$coef)
}
}
beta.curr<-.5*beta.new+.5*beta.curr
beta.new[1]<-mean(y[which.use])-mean(X.new[which.use,]%*%beta.new[-1])
## This is for the intercept. This corrects the demean tricks.
run.diff<-(mean((beta.new[-1]-beta.curr[-1])^2)/sum(beta.new[-1]^2+1e-10))
if(run.diff<1e-6) break
which.use<-(y*cbind(1,X.new)%*%beta.new)<=y^2
}
##Find fitted values.
fits<-pmin(abs(cbind(1,X.new)%*%beta.new),y^2)*sign((cbind(1,X.new)%*%beta.new))
fits2<-cbind(1,X.new)%*%beta.new
fits2[sign(fits2)==sign(y)]<-pmin(abs(fits2[sign(fits2)==sign(y)]),
y[sign(fits2)==sign(y)]^2)*sign(fits2[sign(fits2)==sign(y)])
}
if(type=="continuous"){
X.new.2<-apply(X.new,MARGIN=2,FUN=function(x) x-mean(x))
X.new.2<-apply(X.new.2,MARGIN=2,
FUN=function(x) if(length(unique(x))==1) rnorm(sum(which.use)) else x)
if(fit.glmnet==TRUE){
glmnet.1<-glmnet(X.new.2,y-mean(y),
family="gaussian",
lambda=c(5,4,3,2,seq(1.5,1,-.1)),
standardize=F)
beta.new[-1]<-as.vector(glmnet.1$beta[,10])
}else{
##This was the old L1 optimizer.
lasso1<-lars(X.new.2,(y-mean(y)),
max.steps=dim(X)[2]+2,normalize=F,type="lasso",eps=0)
beta.new[-1]<-as.vector(predict(lasso1,s=1,type="coef",mode="lambda")$coef)
if(log(length(y))/2*sum(beta.new[-1]!=0)>.9*length(y))
beta.new[-1]<-
as.vector(predict(
lasso1,
s=min(c(floor(.9*length(y)*2/log(n)),dim(X.new)[2]*.8)),
type="coef",mode="step")$coef)
}
beta.curr<-.5*beta.new+.5*beta.curr
beta.new[1]<-mean(y)-mean(X.new%*%beta.new[-1])
## This is for the intercept
##run.diff<-(mean((beta.new[-1]-beta.curr[-1])^2)/sum(beta.new[-1]^2+1e-10))
##if(run.diff<1e-6) break
##which.use<-(y*cbind(1,X.new)%*%beta.new)<=y^2
##}
## We dont need loop. Because this is not SVM.
##Find fitted values.
fits2 <- cbind(1,X.new)%*%beta.new
fits <- fits2
##fits<-pmin(abs(cbind(1,X.new)%*%beta.new),y^2)*sign((cbind(1,X.new)%*%beta.new))
##fits2<-cbind(1,X.new)%*%beta.new
##fits2[sign(fits2)==sign(y)]<-pmin(abs(fits2[sign(fits2)==sign(y)]),
## y[sign(fits2)==sign(y)]^2)*sign(fits2[sign(fits2)==sign(y)])
}
##Calculate degrees of freedom
edf<- 1+sum(beta.new[-1]!=0)
##*log(n)/2#sum(diag(hat.mat))#+1-sum(beta.new[-c(1:(dim(X.c)[2]+1))]!=0)
##Left in for estimates of standard deviation.
stdev<-NA
##GCV statistic
if(type=="binary"){
loss.run<-sum((y^2-y*fits2)^2)/(sum(which.use)-edf)^2*(mean(which.use))^2*n
}
if(type=="continuous"){
loss.run<-sum((y-fits2)^2)/(n*((1-edf/n)^2))
}
## Gather minimum loss function.
## if(i*j>1) if(loss.run[i,j]==min(loss.run[loss.run!=0]))
## min.params<-c(lambda.c[i],lambda.t[j])
## }
## }
##Scale betas back.
if(main){
if(is.matrix(X.c)){
beta.new[c(1:dim(X.c)[2]+1)]<-beta.new[c(1:dim(X.c)[2]+1)]/exp(lambda.c)
##for the main covariates and intercepts
beta.new[-c(1:(dim(X.c)[2]+1))]<-beta.new[-c(1:(dim(X.c)[2]+1))]/exp(lambda.t)
##for the interaction covariates
}
if(is.vector(X.c)){
beta.new[c(1:2)]<-beta.new[c(1:2)]/exp(lambda.c)
##for the main covariates and intercepts
beta.new[-c(1:2)]<-beta.new[-c(1:2)]/exp(lambda.t)
##for the interaction covariates
}
}
if(main==FALSE){
beta.new <- beta.new/exp(lambda.t)
}
beta.new<-as.vector(beta.new)
X<-as.matrix(X)
##Calculate intercept
if(type=="binary"){
beta.new[1]<-mean(y[which.use])-mean((X%*%beta.new[-1])[which.use])
}
if(type=="continuous"){
beta.new[1]<-mean(y)-mean((X%*%beta.new[-1]))
}
output<-list(##"lambdas"=min.params,
"beta"=beta.new,"fits"=fits,
"loss"=log(loss.run),##"marg"=mean(which.use),
"edf"=edf,"sd"=stdev)
invisible(output)
}
#####################################################
#####################################################
search.lambda<-function(y=y,X.c=X.c,X.t=X.t,fit.glmnet,main,type){
lambda.find<-function(lambda) as.numeric(SVM.func(y,X.c,X.t,
lambda[1],lambda[2],
fit.glmnet=fit.glmnet,
main=main,
type=type)$loss)
lambda.c.find<-function(lambda) as.numeric(SVM.func(y,X.c,X.t,
lambda,lambda.t,
fit.glmnet=fit.glmnet,
main=main,
type=type)$loss)
lambda.t.find<-function(lambda) as.numeric(SVM.func(y,X.c,X.t,
lambda.c,lambda,
fit.glmnet=fit.glmnet,
main=main,
type=type)$loss)
lambda.c.old<-lambda.t.old<-999
lambda.c.seek<-seq(-15,10,1)
lambda.t.seek<-seq(-15,10,1)
lambda.t<-25
lambda.c.out<-sapply(lambda.c.seek,lambda.c.find)
lambda.c<-lambda.c.seek[min(which(lambda.c.out==min(lambda.c.out)))]
lambda.t.out<-sapply(lambda.t.seek,lambda.t.find)
lambda.t<-lambda.t.seek[min(which(lambda.t.out==min(lambda.t.out)))]
print(c(lambda.c,lambda.t))
print(range(lambda.c.seek))
print(range(lambda.t.seek))
##lambda.c<- -15
##lambda.t<-5
lambda.c.old<-lambda.c
lambda.t.old<-lambda.t
for(lambda.loop in 1:5){
lambda.c.seek<-seq(lambda.c-5,lambda.c+5,1)
lambda.c.seek<-unique(pmax(lambda.c.seek,-15))
lambda.t.seek<-seq(lambda.t-5,lambda.t+5,1)
lambda.c.out<-sapply(lambda.c.seek,lambda.c.find)
lambda.c<-lambda.c.seek[min(which(lambda.c.out==min(lambda.c.out)))]
lambda.c<-max(-15,lambda.c)
lambda.t.out<-sapply(lambda.t.seek,lambda.t.find)
lambda.t<-lambda.t.seek[min(which(lambda.t.out==min(lambda.t.out)))]
print(c(lambda.c,lambda.t))
print(range(lambda.c.seek))
print(range(lambda.t.seek))
if(lambda.c.old==lambda.c&lambda.t.old==lambda.t) break
lambda.c.old<-lambda.c
lambda.t.old<-lambda.t
}
for(lambda.loop in 1:5){
lambda.c.seek<-seq(lambda.c-2.5,lambda.c+2.5,.5)
lambda.c.seek<-unique(pmax(lambda.c.seek,-15))
lambda.t.seek<-seq(lambda.t-2.5,lambda.t+2.5,.5)
lambda.c.out<-sapply(lambda.c.seek,lambda.c.find)
lambda.c<-lambda.c.seek[min(which(lambda.c.out==min(lambda.c.out)))]
lambda.c<-max(-15,lambda.c)
lambda.t.out<-sapply(lambda.t.seek,lambda.t.find)
lambda.t<-lambda.t.seek[min(which(lambda.t.out==min(lambda.t.out)))]
print(c(lambda.c,lambda.t))
print(range(lambda.c.seek))
print(range(lambda.t.seek))
if(lambda.c.old==lambda.c&lambda.t.old==lambda.t) break
lambda.c.old<-lambda.c
lambda.t.old<-lambda.t
}
for(lambda.loop in 1:5){
lambda.c.seek<-seq(lambda.c-1.5,lambda.c+1.5,.25)
lambda.c.seek<-unique(pmax(lambda.c.seek,-15))
lambda.t.seek<-seq(lambda.t-1.5,lambda.t+1.5,.25)
lambda.c.out<-sapply(lambda.c.seek,lambda.c.find)
lambda.c<-lambda.c.seek[min(which(lambda.c.out==min(lambda.c.out)))]
lambda.c<-max(-15,lambda.c)
lambda.t.out<-sapply(lambda.t.seek,lambda.t.find)
lambda.t<-lambda.t.seek[min(which(lambda.t.out==min(lambda.t.out)))]
print(c(lambda.c,lambda.t))
print(range(lambda.c.seek))
print(range(lambda.t.seek))
if(lambda.c.old==lambda.c&lambda.t.old==lambda.t) break
lambda.c.old<-lambda.c
lambda.t.old<-lambda.t
}
for(lambda.loop in 1:5){
lambda.c.seek<-seq(lambda.c-.5,lambda.c+.5,.1)
lambda.c.seek<-unique(pmax(lambda.c.seek,-15))
lambda.t.seek<-seq(lambda.t-.5,lambda.t+.5,.1)
lambda.c.out<-sapply(lambda.c.seek,lambda.c.find)
lambda.c<-lambda.c.seek[min(which(lambda.c.out==min(lambda.c.out)))]
lambda.c<-max(-15,lambda.c)
lambda.t.out<-sapply(lambda.t.seek,lambda.t.find)
lambda.t<-lambda.t.seek[min(which(lambda.t.out==min(lambda.t.out)))]
print(c(lambda.c,lambda.t))
print(range(lambda.c.seek))
print(range(lambda.t.seek))
if(lambda.c.old==lambda.c&lambda.t.old==lambda.t) break
lambda.c.old<-lambda.c
lambda.t.old<-lambda.t
}
for(lambda.loop in 1:5){
lambda.c.seek<-seq(lambda.c-.25,lambda.c+.25,.05)
lambda.c.seek<-unique(pmax(lambda.c.seek,-15))
lambda.t.seek<-seq(lambda.t-.25,lambda.t+.25,.05)
lambda.c.out<-sapply(lambda.c.seek,lambda.c.find)
lambda.c<-lambda.c.seek[min(which(lambda.c.out==min(lambda.c.out)))]
lambda.c<-max(-15,lambda.c)
lambda.t.out<-sapply(lambda.t.seek,lambda.t.find)
lambda.t<-lambda.t.seek[min(which(lambda.t.out==min(lambda.t.out)))]
print(c(lambda.c,lambda.t))
print(range(lambda.c.seek))
print(range(lambda.t.seek))
if(lambda.c.old==lambda.c&lambda.t.old==lambda.t) break
lambda.c.old<-lambda.c
lambda.t.old<-lambda.t
}
for(lambda.loop in 1:5){
lambda.c.seek<-seq(lambda.c-.15,lambda.c+.15,.025)
lambda.c.seek<-unique(pmax(lambda.c.seek,-15))
lambda.t.seek<-seq(lambda.t-.15,lambda.t+.15,.025)
lambda.c.out<-sapply(lambda.c.seek,lambda.c.find)
lambda.c<-lambda.c.seek[min(which(lambda.c.out==min(lambda.c.out)))]
lambda.c<-max(-15,lambda.c)
lambda.t.out<-sapply(lambda.t.seek,lambda.t.find)
lambda.t<-lambda.t.seek[min(which(lambda.t.out==min(lambda.t.out)))]
print(c(lambda.c,lambda.t))
print(range(lambda.c.seek))
print(range(lambda.t.seek))
if(lambda.c.old==lambda.c&lambda.t.old==lambda.t) break
lambda.c.old<-lambda.c
lambda.t.old<-lambda.t
}
for(lambda.loop in 1:5){
lambda.c.seek<-seq(lambda.c-.1,lambda.c+.1,.01)
lambda.c.seek<-unique(pmax(lambda.c.seek,-15))
lambda.t.seek<-seq(lambda.t-.1,lambda.t+.1,.01)
lambda.c.out<-sapply(lambda.c.seek,lambda.c.find)
lambda.c<-lambda.c.seek[min(which(lambda.c.out==min(lambda.c.out)))]
lambda.c<-max(-15,lambda.c)
lambda.t.out<-sapply(lambda.t.seek,lambda.t.find)
lambda.t<-lambda.t.seek[min(which(lambda.t.out==min(lambda.t.out)))]
print(c(lambda.c,lambda.t))
print(range(lambda.c.seek))
print(range(lambda.t.seek))
if(lambda.c.old==lambda.c&lambda.t.old==lambda.t) break
lambda.c.old<-lambda.c
lambda.t.old<-lambda.t
}
for(lambda.loop in 1:5){
lambda.c.seek<-seq(lambda.c-.05,lambda.c+.05,.005)
lambda.c.seek<-unique(pmax(lambda.c.seek,-15))
lambda.t.seek<-seq(lambda.t-.05,lambda.t+.05,.005)
lambda.c.out<-sapply(lambda.c.seek,lambda.c.find)
lambda.c<-lambda.c.seek[min(which(lambda.c.out==min(lambda.c.out)))]
lambda.c<-max(-15,lambda.c)
lambda.t.out<-sapply(lambda.t.seek,lambda.t.find)
lambda.t<-lambda.t.seek[min(which(lambda.t.out==min(lambda.t.out)))]
print(c(lambda.c,lambda.t))
print(range(lambda.c.seek))
print(range(lambda.t.seek))
if(lambda.c.old==lambda.c&lambda.t.old==lambda.t) break
lambda.c.old<-lambda.c
lambda.t.old<-lambda.t
}
for(lambda.loop in 1:5){
lambda.c.seek<-seq(lambda.c-.025,lambda.c+.025,.0025)
lambda.c.seek<-unique(pmax(lambda.c.seek,-15))
lambda.t.seek<-seq(lambda.t-.025,lambda.t+.025,.0025)
lambda.c.out<-sapply(lambda.c.seek,lambda.c.find)
lambda.c<-lambda.c.seek[min(which(lambda.c.out==min(lambda.c.out)))]
lambda.c<-max(-15,lambda.c)
lambda.t.out<-sapply(lambda.t.seek,lambda.t.find)
lambda.t<-lambda.t.seek[min(which(lambda.t.out==min(lambda.t.out)))]
print(c(lambda.c,lambda.t))
print(range(lambda.c.seek))
print(range(lambda.t.seek))
if(lambda.c.old==lambda.c&lambda.t.old==lambda.t) break
lambda.c.old<-lambda.c
lambda.t.old<-lambda.t
}
for(lambda.loop in 1:5){
lambda.c.seek<-seq(lambda.c-.01,lambda.c+.01,.001)
lambda.c.seek<-unique(pmax(lambda.c.seek,-15))
lambda.t.seek<-seq(lambda.t-.01,lambda.t+.01,.001)
lambda.c.out<-sapply(lambda.c.seek,lambda.c.find)
lambda.c<-lambda.c.seek[min(which(lambda.c.out==min(lambda.c.out)))]
lambda.c<-max(-15,lambda.c)
lambda.t.out<-sapply(lambda.t.seek,lambda.t.find)
lambda.t<-lambda.t.seek[min(which(lambda.t.out==min(lambda.t.out)))]
print(c(lambda.c,lambda.t))
print(range(lambda.c.seek))
print(range(lambda.t.seek))
if(lambda.c.old==lambda.c&lambda.t.old==lambda.t) break
lambda.c.old<-lambda.c
lambda.t.old<-lambda.t
}
for(lambda.loop in 1:5){
lambda.c.seek<-seq(lambda.c-.0025,lambda.c+.0025,.0005)
lambda.c.seek<-unique(pmax(lambda.c.seek,-15))
lambda.t.seek<-seq(lambda.t-.0025,lambda.t+.0025,.0005)
lambda.c.out<-sapply(lambda.c.seek,lambda.c.find)
lambda.c<-lambda.c.seek[min(which(lambda.c.out==min(lambda.c.out)))]
lambda.c<-max(-15,lambda.c)
lambda.t.out<-sapply(lambda.t.seek,lambda.t.find)
lambda.t<-lambda.t.seek[min(which(lambda.t.out==min(lambda.t.out)))]
print(c(lambda.c,lambda.t))
print(range(lambda.c.seek))
print(range(lambda.t.seek))
if(lambda.c.old==lambda.c&lambda.t.old==lambda.t) break
lambda.c.old<-lambda.c
lambda.t.old<-lambda.t
}
optim.lambda<-c(lambda.c,lambda.t)
invisible(optim.lambda)
}
#####################################################
search.lambda.nocov<-function(y=y,X.t=X.t,fit.glmnet,main,type){
lambda.t.find<-function(lambda) as.numeric(SVM.func(y,X.c=NULL,X.t,
lambda.c=1,lambda,
fit.glmnet=fit.glmnet,
main=FALSE,
type=type)$loss)
## This lambda.c will not be used.
## lambda.c.old<-lambda.t.old<-999
lambda.t.old<-999
##############################
## The Initial Try
##############################
## lambda.c.seek<-seq(-15,10,1)
lambda.t.seek<-seq(-15,10,1)
## lambda.t<-25
## lambda.c.out<-sapply(lambda.c.seek,lambda.c.find)
## lambda.c<-lambda.c.seek[min(which(lambda.c.out==min(lambda.c.out)))]
lambda.t.out<-sapply(lambda.t.seek,lambda.t.find)
lambda.t<-lambda.t.seek[min(which(lambda.t.out==min(lambda.t.out)))]
## print(c(lambda.c,lambda.t))
## print(range(lambda.c.seek))
print(range(lambda.t.seek))
##lambda.c<- -15
##lambda.t<-5
## lambda.c.old<-lambda.c
lambda.t.old<-lambda.t
for(lambda.loop in 1:5){
## lambda.c.seek<-seq(lambda.c-5,lambda.c+5,1)
## lambda.c.seek<-unique(pmax(lambda.c.seek,-15))
lambda.t.seek<-seq(lambda.t-5,lambda.t+5,1)
## lambda.c.out<-sapply(lambda.c.seek,lambda.c.find)
## lambda.c<-lambda.c.seek[min(which(lambda.c.out==min(lambda.c.out)))]
## lambda.c<-max(-15,lambda.c)
lambda.t.out<-sapply(lambda.t.seek,lambda.t.find)
lambda.t<-lambda.t.seek[min(which(lambda.t.out==min(lambda.t.out)))]
## print(c(lambda.c,lambda.t))
## print(range(lambda.c.seek))
print(range(lambda.t.seek))
## if(lambda.c.old==lambda.c&lambda.t.old==lambda.t) break
if(lambda.t.old==lambda.t) break
## lambda.c.old<-lambda.c
lambda.t.old<-lambda.t
}
for(lambda.loop in 1:5){
## lambda.c.seek<-seq(lambda.c-2.5,lambda.c+2.5,.5)
## lambda.c.seek<-unique(pmax(lambda.c.seek,-15))
lambda.t.seek<-seq(lambda.t-2.5,lambda.t+2.5,.5)
## lambda.c.out<-sapply(lambda.c.seek,lambda.c.find)
## lambda.c<-lambda.c.seek[min(which(lambda.c.out==min(lambda.c.out)))]
## lambda.c<-max(-15,lambda.c)
lambda.t.out<-sapply(lambda.t.seek,lambda.t.find)
lambda.t<-lambda.t.seek[min(which(lambda.t.out==min(lambda.t.out)))]
## print(c(lambda.c,lambda.t))
## print(range(lambda.c.seek))
print(range(lambda.t.seek))
## if(lambda.c.old==lambda.c&lambda.t.old==lambda.t) break
if(lambda.t.old==lambda.t) break
## lambda.c.old<-lambda.c
lambda.t.old<-lambda.t
}
for(lambda.loop in 1:5){
## lambda.c.seek<-seq(lambda.c-1.5,lambda.c+1.5,.25)
## lambda.c.seek<-unique(pmax(lambda.c.seek,-15))
lambda.t.seek<-seq(lambda.t-1.5,lambda.t+1.5,.25)
## lambda.c.out<-sapply(lambda.c.seek,lambda.c.find)
## lambda.c<-lambda.c.seek[min(which(lambda.c.out==min(lambda.c.out)))]
## lambda.c<-max(-15,lambda.c)
lambda.t.out<-sapply(lambda.t.seek,lambda.t.find)
lambda.t<-lambda.t.seek[min(which(lambda.t.out==min(lambda.t.out)))]
## print(c(lambda.c,lambda.t))
## print(range(lambda.c.seek))
print(range(lambda.t.seek))
## if(lambda.c.old==lambda.c&lambda.t.old==lambda.t) break
if(lambda.t.old==lambda.t) break
## lambda.c.old<-lambda.c
lambda.t.old<-lambda.t
}
for(lambda.loop in 1:5){
## lambda.c.seek<-seq(lambda.c-.5,lambda.c+.5,.1)
## lambda.c.seek<-unique(pmax(lambda.c.seek,-15))
lambda.t.seek<-seq(lambda.t-.5,lambda.t+.5,.1)
## lambda.c.out<-sapply(lambda.c.seek,lambda.c.find)
## lambda.c<-lambda.c.seek[min(which(lambda.c.out==min(lambda.c.out)))]
## lambda.c<-max(-15,lambda.c)
lambda.t.out<-sapply(lambda.t.seek,lambda.t.find)
lambda.t<-lambda.t.seek[min(which(lambda.t.out==min(lambda.t.out)))]
## print(c(lambda.c,lambda.t))
## print(range(lambda.c.seek))
print(range(lambda.t.seek))
## if(lambda.c.old==lambda.c&lambda.t.old==lambda.t) break
if(lambda.t.old==lambda.t) break
## lambda.c.old<-lambda.c
lambda.t.old<-lambda.t
}
for(lambda.loop in 1:5){
## lambda.c.seek<-seq(lambda.c-.25,lambda.c+.25,.05)
## lambda.c.seek<-unique(pmax(lambda.c.seek,-15))
lambda.t.seek<-seq(lambda.t-.25,lambda.t+.25,.05)
## lambda.c.out<-sapply(lambda.c.seek,lambda.c.find)
## lambda.c<-lambda.c.seek[min(which(lambda.c.out==min(lambda.c.out)))]
## lambda.c<-max(-15,lambda.c)
lambda.t.out<-sapply(lambda.t.seek,lambda.t.find)
lambda.t<-lambda.t.seek[min(which(lambda.t.out==min(lambda.t.out)))]
## print(c(lambda.c,lambda.t))
## print(range(lambda.c.seek))
print(range(lambda.t.seek))
## if(lambda.c.old==lambda.c&lambda.t.old==lambda.t) break
if(lambda.t.old==lambda.t) break
## lambda.c.old<-lambda.c
lambda.t.old<-lambda.t
}
for(lambda.loop in 1:5){
## lambda.c.seek<-seq(lambda.c-.15,lambda.c+.15,.025)
## lambda.c.seek<-unique(pmax(lambda.c.seek,-15))
lambda.t.seek<-seq(lambda.t-.15,lambda.t+.15,.025)
## lambda.c.out<-sapply(lambda.c.seek,lambda.c.find)
## lambda.c<-lambda.c.seek[min(which(lambda.c.out==min(lambda.c.out)))]
## lambda.c<-max(-15,lambda.c)
lambda.t.out<-sapply(lambda.t.seek,lambda.t.find)
lambda.t<-lambda.t.seek[min(which(lambda.t.out==min(lambda.t.out)))]
## print(c(lambda.c,lambda.t))
## print(range(lambda.c.seek))
print(range(lambda.t.seek))
## if(lambda.c.old==lambda.c&lambda.t.old==lambda.t) break
if(lambda.t.old==lambda.t) break
## lambda.c.old<-lambda.c
lambda.t.old<-lambda.t
}
for(lambda.loop in 1:5){
## lambda.c.seek<-seq(lambda.c-.1,lambda.c+.1,.01)
## lambda.c.seek<-unique(pmax(lambda.c.seek,-15))
lambda.t.seek<-seq(lambda.t-.1,lambda.t+.1,.01)
## lambda.c.out<-sapply(lambda.c.seek,lambda.c.find)
## lambda.c<-lambda.c.seek[min(which(lambda.c.out==min(lambda.c.out)))]
## lambda.c<-max(-15,lambda.c)
lambda.t.out<-sapply(lambda.t.seek,lambda.t.find)
lambda.t<-lambda.t.seek[min(which(lambda.t.out==min(lambda.t.out)))]
## print(c(lambda.c,lambda.t))
## print(range(lambda.c.seek))
print(range(lambda.t.seek))
## if(lambda.c.old==lambda.c&lambda.t.old==lambda.t) break
if(lambda.t.old==lambda.t) break
## lambda.c.old<-lambda.c
lambda.t.old<-lambda.t
}
for(lambda.loop in 1:5){
## lambda.c.seek<-seq(lambda.c-.05,lambda.c+.05,.005)
## lambda.c.seek<-unique(pmax(lambda.c.seek,-15))
lambda.t.seek<-seq(lambda.t-.05,lambda.t+.05,.005)
## lambda.c.out<-sapply(lambda.c.seek,lambda.c.find)
## lambda.c<-lambda.c.seek[min(which(lambda.c.out==min(lambda.c.out)))]
## lambda.c<-max(-15,lambda.c)
lambda.t.out<-sapply(lambda.t.seek,lambda.t.find)
lambda.t<-lambda.t.seek[min(which(lambda.t.out==min(lambda.t.out)))]
## print(c(lambda.c,lambda.t))
## print(range(lambda.c.seek))
print(range(lambda.t.seek))
## if(lambda.c.old==lambda.c&lambda.t.old==lambda.t) break
if(lambda.t.old==lambda.t) break
## lambda.c.old<-lambda.c
lambda.t.old<-lambda.t
}
for(lambda.loop in 1:5){
## lambda.c.seek<-seq(lambda.c-.025,lambda.c+.025,.0025)
## lambda.c.seek<-unique(pmax(lambda.c.seek,-15))
lambda.t.seek<-seq(lambda.t-.025,lambda.t+.025,.0025)
## lambda.c.out<-sapply(lambda.c.seek,lambda.c.find)
## lambda.c<-lambda.c.seek[min(which(lambda.c.out==min(lambda.c.out)))]
## lambda.c<-max(-15,lambda.c)
lambda.t.out<-sapply(lambda.t.seek,lambda.t.find)
lambda.t<-lambda.t.seek[min(which(lambda.t.out==min(lambda.t.out)))]
## print(c(lambda.c,lambda.t))
## print(range(lambda.c.seek))
print(range(lambda.t.seek))
## if(lambda.c.old==lambda.c&lambda.t.old==lambda.t) break
if(lambda.t.old==lambda.t) break
## lambda.c.old<-lambda.c
lambda.t.old<-lambda.t
}
for(lambda.loop in 1:5){
## lambda.c.seek<-seq(lambda.c-.01,lambda.c+.01,.001)
## lambda.c.seek<-unique(pmax(lambda.c.seek,-15))
lambda.t.seek<-seq(lambda.t-.01,lambda.t+.01,.001)
## lambda.c.out<-sapply(lambda.c.seek,lambda.c.find)
## lambda.c<-lambda.c.seek[min(which(lambda.c.out==min(lambda.c.out)))]
## lambda.c<-max(-15,lambda.c)
lambda.t.out<-sapply(lambda.t.seek,lambda.t.find)
lambda.t<-lambda.t.seek[min(which(lambda.t.out==min(lambda.t.out)))]
## print(c(lambda.c,lambda.t))
## print(range(lambda.c.seek))
print(range(lambda.t.seek))
## if(lambda.c.old==lambda.c&lambda.t.old==lambda.t) break
if(lambda.t.old==lambda.t) break
## lambda.c.old<-lambda.c
lambda.t.old<-lambda.t
}
for(lambda.loop in 1:5){
## lambda.c.seek<-seq(lambda.c-.0025,lambda.c+.0025,.0005)
## lambda.c.seek<-unique(pmax(lambda.c.seek,-15))
lambda.t.seek<-seq(lambda.t-.0025,lambda.t+.0025,.0005)
## lambda.c.out<-sapply(lambda.c.seek,lambda.c.find)
## lambda.c<-lambda.c.seek[min(which(lambda.c.out==min(lambda.c.out)))]
## lambda.c<-max(-15,lambda.c)
lambda.t.out<-sapply(lambda.t.seek,lambda.t.find)
lambda.t<-lambda.t.seek[min(which(lambda.t.out==min(lambda.t.out)))]
## print(c(lambda.c,lambda.t))
## print(range(lambda.c.seek))
print(range(lambda.t.seek))
## if(lambda.c.old==lambda.c&lambda.t.old==lambda.t) break
if(lambda.t.old==lambda.t) break
## lambda.c.old<-lambda.c
lambda.t.old<-lambda.t
}
## optim.lambda<-c(lambda.c,lambda.t)
optim.lambda<- lambda.t
invisible(optim.lambda)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.