# R/gqTHOR.R In scottwalmsley/THOR: Determines a spectral count threshold cutoff.

#### Defines functions gq.THOR

```#' gq.THOR
#' Runs the Goldfeld-QUandt test for heteroskedasticity in a sliding window test.
#'
#' @param THOR.data
#' @param point
#' @param fraction
#' @param STEP
#' @param alpha
#' @param PLOT
#'
#' @return
#' @export
#'
#' @examples
gq.THOR = function (THOR.data, point = 0.5, fraction = 0.375,STEP = 1,alpha=0.05,PLOT=T){

require("lmtest")
nam = array();
RVAL = list();
if(PLOT==T){
pdf(file="M_THOR.pdf",onefile=TRUE,width=8,height=3)
}

for(N in 1:length(THOR.data\$r.M)){
cat(paste("N=",N,"\n"));
pvals =  array();
GQ = array();
SpC = array();

r.M = THOR.data\$r.M[[N]];
M = THOR.data\$M[[N]];
M

n = dim(r.M)[1];
nr = n - (n %% STEP);
nr = nr - 100;
id = seq(1,nr,by=STEP);

j=1;i=1;
for(i in id) {

fit = lm(Ref~y-1,data=r.M[i:n,]);
gq = gqtest(fit,alternative="greater");
GQ[j] = gq\$statistic;
pvals[j] = gq\$p.value;
SpC[j] = exp(M\$Ref[i]);
j=j+1;

}

csum = cumsum(pvals)/sum(pvals);
cor = cor(M\$Ref,M\$y,method="pearson");

CO.id = min(id[csum >= alpha]);
rM.CO = r.M\$Ref[CO.id];
M.CO =  exp(M\$Ref[CO.id]);

if(PLOT==T){

############ estimate regression,
############ and get 95% prediction interval.

w = r.M\$Ref

lm.fit = lm(y ~ Ref, data=r.M, weights=w)
lm.pred = predict(lm.fit, data=r.M, interval="predict", level=0.95, weights=w)
lm.pred = as.data.frame(lm.pred)

#dev.new(width=8,height=3.5);
par(mar=c(4,4,1,1))

plot(y ~ Ref,xlim=c(0,8),ylim=c(-3,3), data=r.M,cex=0.5,xlab="Ref",ylab=paste("R",N,sep=""))
abline(v=rM.CO, col=2)
lines(lm.pred\$fit ~ r.M\$Ref, col=3, lwd=1)
lines(lm.pred\$lwr ~ r.M\$Ref, col=2, lty=2)
lines(lm.pred\$upr ~ r.M\$Ref, col=2, lty=2)

### identify the points
is.outlier = rep(NA, n)

for(i in 1:n) {
if(r.M\$y[i] >= lm.pred\$lwr[i]  &  r.M\$y[i] <= lm.pred\$upr[i]) {
is.outlier[i] = FALSE
}
else{is.outlier[i] = TRUE}
}
points(y ~ Ref, cex=0.5,data=r.M[is.outlier, ], col=4);
}

RVAL[[N]] <- list(Ref=exp(M\$Ref),SpC = exp(M\$y),p.value = pvals,csum = csum,GQ = GQ,count = SpC,rm.co=rM.CO,m.co=M.CO,cor=cor);
#names(RVAL[[N]]) = paste("R",N,sep="");
}

if(PLOT==T){
dev.off();
}

#RVAL = list(Ref=M\$Ref,Result=RVAL);
class(RVAL) = "THOR.result";
#names(RVAL) = nam;
return(RVAL)

}
```
scottwalmsley/THOR documentation built on May 11, 2018, 12:04 a.m.