### R code from vignette source 'vignettes/flowQB/inst/doc/AdvancedflowQBNIH3.Rnw'
### Encoding: UTF-8

###################################################
### code chunk number 1: tab1
###################################################
rm(list=ls(all=TRUE))
library("flowQB")
File= system.file("extdata","NIH3.fcs",package="flowQB")
library(flowCore)
Data=read.FCS(File)


###################################################
### code chunk number 2: tab1
###################################################
Markers=c("FSC-A","SSC-A","B710-A","B515-A","G780-A","G710-A","G660-A","G610-A" 
,"G560-A", "R780-A","R710-A","R660-A","V800-A","V705-A","V655-A"
,"V605-A","V585-A","V565-A","V545-A","V450-A")
lgcl <- estimateLogicle(Data, channels = Markers)
TransformedData <- transform(Data, lgcl)


###################################################
### code chunk number 3: tab1
###################################################
TransformedMFI=data.frame(exprs(TransformedData))


###################################################
### code chunk number 4: tab1
###################################################
OriginalMFI=data.frame(exprs(Data))


###################################################
### code chunk number 5: tab1
###################################################
a1=3.3
b1=3.6
a2=2.5
b2=2.85
SingletsIndices=which(TransformedMFI[,1]> a1 & TransformedMFI[,1] < b1 
 & TransformedMFI[,3]> a2 & TransformedMFI[,3] < b2)


###################################################
### code chunk number 6: tab1
###################################################
nClusters=7


###################################################
### code chunk number 7: tab1
###################################################
MarkersForAssessment=c(3,4,5,6,7)
TSSINGLETS=TransformedMFI[SingletsIndices,MarkersForAssessment]
OSSINGLETS=data.frame(OriginalMFI[SingletsIndices,MarkersForAssessment])



###################################################
### code chunk number 8: tab1
###################################################
gps=gPS(TSSINGLETS[,c(1,3)],OSSINGLETS[,c(1,3)],7)



###################################################
### code chunk number 9: tab1
###################################################
rps=rPS(TSSINGLETS[,c(1,3)],OSSINGLETS[,c(1,3)],7)



###################################################
### code chunk number 10: tab1
###################################################
library(xtable)
PSdataB515=data.frame(NE=gps[,1]
,gMeans=gps[,2]
,rMeans=rps[,2]
,gSDs=gps[,3]
,rSDs=rps[,3]
,gRDs=gps[,4]
,rRDs=rps[,4])
row.names(PSdataB515)=paste("Peak",1:nClusters,sep="")
print(xtable(PSdataB515, caption = 
"Extracted Gaussian (g), Robust (r) peak statistics and associated RDs. 
NE is the number of events.", label = "tab:one"))


###################################################
### code chunk number 11: tab1
###################################################
mrps=MultilevelrPS(TSSINGLETS[,2:4],OSSINGLETS[,2:4],7)
mgps=MultilevelgPS(TSSINGLETS[,2:4],OSSINGLETS[,2:4],7)
library(xtable)
print(xtable(mrps, caption = "Multi level Kmeans data with Robust (r) statistics  for  peaks. 
NE is the number of events.", label = "tab:mr"))
print(xtable(mgps, caption = 
"Multi level Kmeans data with Gaussian (g) statistics  for  peaks. 
NE is the number of events.", label = "tab:mg"))


###################################################
### code chunk number 12: tab1
###################################################
pi=1
pf=6


###################################################
### code chunk number 13: tab1
###################################################
rPSdataB515=data.frame(NE=rps[,1],rMeans=rps[,2]
 ,rSDs=rps[,3])
CVs=rPSdataB515[,3]/rPSdataB515[,2]


###################################################
### code chunk number 14: tab1
###################################################
p=1
Q2D=MFI2MESF(rPSdataB515[,2:3],p,0)


###################################################
### code chunk number 15: tab1
###################################################
NE=rPSdataB515[pi:pf,1]
S=Q2D[pi:pf,2]
c0=0
c1=1
c2=0
WWs=0.5*(NE-1)/((1+c2*S^2+c1*S+c0)^2)
for( iter in 1:6)
{
QQBw=qrWEIGHTEDMESF(Q2D,pi,pf,WWs)
COFS=c(c0,as.numeric(QQBw)[5],as.numeric(QQBw)[8],c1,as.numeric(QQBw)[6],as.numeric(QQBw)[9],100*c2,as.numeric(QQBw)[7],as.numeric(QQBw)[10])
if(iter==2) COEFFS=COFS
if(iter>2) COEFFS=rbind(COEFFS,COFS)
c1=1/as.numeric(QQBw)[1]
c2=as.numeric(QQBw)[4]
c0=as.numeric(QQBw)[2]/as.numeric(QQBw)[1]
WWs=0.5*(NE-1)/((1+c2*S^2+c1*S+c0)^2)
}


###################################################
### code chunk number 16: tab1
###################################################
print(WWs)



###################################################
### code chunk number 17: tab1
###################################################
COEFFS=data.frame(COEFFS)
names(COEFFS)= c( "c0" , "Pc0" ,"StdErrorc0", "c1" , "Pc1","StdErrorc1","c2percent" , "Pc2","StdErrorc2" )
row.names(COEFFS)=paste("Iter",1:5,sep="")
print(xtable(COEFFS, caption = "Computed coefficients $(c_{0},c_{1},c_{2})$ 
using weighted quadratic regression where $Pc_{i}$ is the pvalue associated to $c_{i}$ and $StdErrorc_{i}$ is Std. Error of $c_{i}$.", label = "tab:3"))


###################################################
### code chunk number 18: tab1
###################################################
Q=1/c1
B=c0/c1
CVintrinsic=sqrt(c2)
print(paste("Detection efficiency-Q=",round(Q,4)))
print(paste("Background illumination-B=",round(B,4)))
print(paste("Intrinsic CV of the beads-CVintrinsic=",round(CVintrinsic,4)))


###################################################
### code chunk number 19: tab1
###################################################
delta=c1^2-c0*c2*4
R1=-c1-sqrt(delta)
R1=R1/(2*c2)
R2=-c1+sqrt(delta)
R2=R2/(2*c2)
print(paste("Discriminant Delta=",round(delta,4)))
print(paste("Root1=",round(R1,4)))
print(paste("Root2=",round(R2,4)))


###################################################
### code chunk number 20: tab1
###################################################
mrPSdataB515=data.frame(NE=mrps[,1],rMeans=mrps[,3]
 ,rSDs=mrps[,6])
CVs=mrPSdataB515[,3]/mrPSdataB515[,2]


###################################################
### code chunk number 21: tab1
###################################################
p=1
Q2D=MFI2MESF(mrPSdataB515[,2:3],p,0)


###################################################
### code chunk number 22: tab1
###################################################
NE=mrPSdataB515[pi:pf,1]
S=Q2D[pi:pf,2]
c0=0
c1=1
c2=0
WWs=0.5*(NE-1)/((1+c2*S^2+c1*S+c0)^2)
for( iter in 1:6)
{
QQBw=qrWEIGHTEDMESF(Q2D,pi,pf,WWs)
COFS=c(c0,as.numeric(QQBw)[5],as.numeric(QQBw)[8],c1,as.numeric(QQBw)[6],as.numeric(QQBw)[9],100*c2,as.numeric(QQBw)[7],as.numeric(QQBw)[10])
if(iter==2) COEFFS=COFS
if(iter>2) COEFFS=rbind(COEFFS,COFS)
c1=1/as.numeric(QQBw)[1]
c2=as.numeric(QQBw)[4]
c0=as.numeric(QQBw)[2]/as.numeric(QQBw)[1]
WWs=0.5*(NE-1)/((1+c2*S^2+c1*S+c0)^2)
}


###################################################
### code chunk number 23: tab1
###################################################
print(WWs)



###################################################
### code chunk number 24: tab1
###################################################
COEFFS=data.frame(COEFFS)
names(COEFFS)= c( "c0" , "Pc0" ,"StdErrorc0", "c1" , "Pc1","StdErrorc1","c2percent" , "Pc2","StdErrorc2" )
row.names(COEFFS)=paste("Iter",1:5,sep="")
print(xtable(COEFFS, caption = "Multi-level Kmeans-Robust: Computed coefficients $(c_{0},c_{1},c_{2})$ 
using weighted quadratic regression where $Pc_{i}$ is the pvalue associated to $c_{i}$ and $StdErrorc_{i}$ is Std. Error of $c_{i}$.", label = "tab:mrq"))


###################################################
### code chunk number 25: tab1
###################################################
Q=1/c1
B=c0/c1
CVintrinsic=sqrt(c2)
print(paste("Detection efficiency-Q=",round(Q,4)))
print(paste("Background illumination-B=",round(B,4)))
print(paste("Intrinsic CV of the beads-CVintrinsic=",round(CVintrinsic,4)))


###################################################
### code chunk number 26: tab1
###################################################
delta=c1^2-c0*c2*4
R1=-c1-sqrt(delta)
R1=R1/(2*c2)
R2=-c1+sqrt(delta)
R2=R2/(2*c2)
print(paste("Discriminant Delta=",round(delta,4)))
print(paste("Root1=",round(R1,4)))
print(paste("Root2=",round(R2,4)))


###################################################
### code chunk number 27: tab1
###################################################
mgPSdataB515=data.frame(NE=mgps[,1],rMeans=mgps[,3]
 ,rSDs=mgps[,6])
CVs=mgPSdataB515[,3]/mgPSdataB515[,2]


###################################################
### code chunk number 28: tab1
###################################################
p=1
Q2D=MFI2MESF(mgPSdataB515[,2:3],p,0)


###################################################
### code chunk number 29: tab1
###################################################
NE=mgPSdataB515[pi:pf,1]
S=Q2D[pi:pf,2]
c0=0
c1=1
c2=0
WWs=0.5*(NE-1)/((1+c2*S^2+c1*S+c0)^2)
for( iter in 1:6)
{
QQBw=qrWEIGHTEDMESF(Q2D,pi,pf,WWs)
COFS=c(c0,as.numeric(QQBw)[5],as.numeric(QQBw)[8],c1,as.numeric(QQBw)[6],as.numeric(QQBw)[9],100*c2,as.numeric(QQBw)[7],as.numeric(QQBw)[10])
if(iter==2) COEFFS=COFS
if(iter>2) COEFFS=rbind(COEFFS,COFS)
c1=1/as.numeric(QQBw)[1]
c2=as.numeric(QQBw)[4]
c0=as.numeric(QQBw)[2]/as.numeric(QQBw)[1]
WWs=0.5*(NE-1)/((1+c2*S^2+c1*S+c0)^2)
}


###################################################
### code chunk number 30: tab1
###################################################
print(WWs)



###################################################
### code chunk number 31: tab1
###################################################
COEFFS=data.frame(COEFFS)
names(COEFFS)= c( "c0" , "Pc0" ,"StdErrorc0", "c1" , "Pc1","StdErrorc1","c2percent" , "Pc2","StdErrorc2" )
row.names(COEFFS)=paste("Iter",1:5,sep="")
print(xtable(COEFFS, caption = "Multi-level Kmeans-Gaussian: Computed coefficients $(c_{0},c_{1},c_{2})$ 
using weighted quadratic regression where $Pc_{i}$ is the pvalue associated to $c_{i}$ and $StdErrorc_{i}$ is Std. Error of $c_{i}$.", label = "tab:mgq"))


###################################################
### code chunk number 32: tab1
###################################################
Q=1/c1
B=c0/c1
CVintrinsic=sqrt(c2)
print(paste("Detection efficiency-Q=",round(Q,4)))
print(paste("Background illumination-B=",round(B,4)))
print(paste("Intrinsic CV of the beads-CVintrinsic=",round(CVintrinsic,4)))


###################################################
### code chunk number 33: tab1
###################################################
delta=c1^2-c0*c2*4
R1=-c1-sqrt(delta)
R1=R1/(2*c2)
R2=-c1+sqrt(delta)
R2=R2/(2*c2)
print(paste("Discriminant Delta=",round(delta,4)))
print(paste("Root1=",round(R1,4)))
print(paste("Root2=",round(R2,4)))