Nothing
`elrm` <-
function(formula, interest, r=4, iter=1000, dataset, burnIn=0, alpha=0.05)
{
if(iter %% 1 != 0)
{
stop("'iter' must be an integer");
}
if(iter < 1000)
{
stop("'iter' must be atleast 1,000");
}
if(burnIn %% 1 != 0)
{
stop("'burnIn' must be an integer");
}
if(iter < burnIn + 1000)
{
stop("'burnIn' must be atleast 1,000 iterations smaller than 'iter'");
}
info = getDesignMatrix(formula,interest,dataset=dataset);
design.matrix = info$design.matrix;
yCol = info$yCol;
mCol = info$mCol;
zCols = info$zCols;
wCols = info$wCols;
if((r %% 2) != 0)
{
stop("'r' must be an even number");
}
if(r > length(as.vector(design.matrix[,yCol])))
{
msg = paste("'r' must be smaller than or equal to",length(as.vector(design.matrix[,yCol])),"(length of the response vector)");
stop(msg);
}
out.filename = paste(tempdir(),"//elrmout",round(runif(1,0,10000),0),".txt",sep="");
tempdata.filename = paste(tempdir(),"//elrmtempdata",round(runif(1,0,10000),0),".txt",sep="");
write.table(design.matrix,tempdata.filename,quote=FALSE,row.names=FALSE,col.names=FALSE);
Z.matrix = matrix(nrow = nrow(design.matrix), ncol = length(zCols));
Z.matrix = as.matrix(design.matrix[,zCols]);
observed.response = as.matrix(design.matrix[,yCol]);
S.matrix = NULL;
message("Generating the Markov chain ...");
t1 = format(Sys.time(),"");
sample.size = max(floor(iter*0.05),1);
last.sample = observed.response;
k = 0;
while(k < iter)
{
progress((k/iter)*100);
Sys.sleep(0.10);
design.matrix[,yCol]=last.sample;
write.table(design.matrix,tempdata.filename,quote=FALSE,row.names=FALSE,col.names=FALSE);
if(k == 0)
{
.C("MCMC", as.integer(yCol), as.integer(mCol), as.integer(wCols), as.integer(length(wCols)), as.integer(zCols), as.integer(length(zCols)), as.integer(r), gsub("\\\\", "/", as.character(tempdata.filename)), gsub("\\\\", "/", as.character(out.filename)), as.integer(min(sample.size,iter-k)), as.integer(1), PACKAGE="elrm");
}
else
{
.C("MCMC", as.integer(yCol), as.integer(mCol), as.integer(wCols), as.integer(length(wCols)), as.integer(zCols), as.integer(length(zCols)), as.integer(r), gsub("\\\\", "/", as.character(tempdata.filename)), gsub("\\\\", "/", as.character(out.filename)), as.integer(min(sample.size,iter-k)), as.integer(0), PACKAGE="elrm");
}
mc = matrix(scan(out.filename,what=numeric(),skip=0,n=min(sample.size,iter-k)*nrow(dataset),sep="\t",quiet=T),nrow=min(sample.size,iter-k),ncol=nrow(dataset),byrow=T);
S.temp = mc %*% Z.matrix
k = k + sample.size;
S.matrix = rbind(S.matrix,S.temp);
last.sample = mc[nrow(mc),];
}
S.matrix = round(S.matrix,6);
progress((k/iter)*100);
message('\n');
Sys.sleep(0.10);
unlink(tempdata.filename);
unlink(out.filename);
t2 = format(Sys.time(),"");
dif1 = difftime(t2,t1,units="auto");
message(paste("Generation of the Markov Chain required", round(dif1,4), attr(dif1,"units")));
message("Conducting inference ...");
t3 = format(Sys.time(),"");
S.observed = round(as.vector(t(Z.matrix)%*%observed.response,mode='numeric'),6);
names(S.observed) = names(design.matrix)[zCols];
S.begin = as.matrix(S.matrix[0:burnIn,]);
if(burnIn == 1)
{
S.begin = as.matrix(t(S.begin));
}
S.matrix = as.matrix(S.matrix[(burnIn+1):iter,]);
marginal = getMarginalInf(S.matrix, S.observed, alpha);
if(length(zCols) > 1)
{
joint = getJointInf(S.matrix, S.observed);
t4 = format(Sys.time(),"");
dif2 = difftime(t4,t3,units="auto");
message(paste("Inference required", round(dif2,4), attr(dif2,"units")));
coeffs = as.vector(c(NA,marginal$estimate),mode='numeric');
names(coeffs) = c("joint",names(marginal$estimate));
pvals = as.vector(c(joint$pvalue,marginal$pvalue),mode='numeric');
pvals.se = as.vector(c(joint$pvalue.se,marginal$pvalue.se),mode='numeric');
names(pvals) = names(coeffs);
names(pvals.se) = names(coeffs);
marginal$distribution[["joint"]] = joint$distribution;
sizes = c(joint$mc.size,marginal$mc.size);
names(sizes) = c("joint",names(marginal$mc.size));
S.matrix = data.frame(rbind(S.begin,S.matrix));
colnames(S.matrix) = names(design.matrix)[zCols];
## S3 definition
object <- list(coeffs=coeffs,coeffs.ci=marginal$estimate.ci,p.values=pvals,p.values.se=pvals.se,mc=mcmc(S.matrix),mc.size=sizes,obs.suff.stat=S.observed,distribution=marginal$distribution,call.history=list(match.call()),dataset=dataset,last=last.sample,r=r,ci.level=(1-alpha)*100);
class(object) <- "elrm";
}
else
{
t4 = format(Sys.time(),"");
dif2 = difftime(t4,t3,units="auto");
message(paste("Inference required", round(dif2,4), attr(dif2,"units")));
S.matrix = data.frame(rbind(S.begin,S.matrix));
colnames(S.matrix) = names(design.matrix)[zCols];
## S3 definition
object <- list(coeffs=marginal$estimate,coeffs.ci=marginal$estimate.ci,p.values=marginal$pvalue,p.values.se=marginal$pvalue.se,mc=mcmc(S.matrix),mc.size=marginal$mc.size,obs.suff.stat=S.observed,distribution=marginal$distribution,call.history=list(match.call()),dataset=dataset,last=last.sample,r=r,ci.level=(1-alpha)*100);
class(object) <- "elrm";
}
return(object);
}
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.