####################################################
##
## R program to read in the melanoma data,
## obtain the deterministic peaks,
## and then fit the multivariate linear mixed effect 
## models using lme, as well as a standard 
## two-stage t test
##
## Browne, Dryden, Handley, Mian, Schadendorf
##                2009
##
######################################################


print("Which section of the spectrum? (1-10) or (0) for the full spectrum")
ans2<-readline()
chunk<- as.integer(ans2)


#print("Do you want to read in all the data? (y/n)")
ans1<-"y"

if (ans1=="y"){

x<-scan(file="melanoma.csv",what="",sep=",")

### 205 observations (101 stage 1, then 104 stage 4

y<-x[207:length(x)]
y<-as.double(y)
z<-matrix(y,206,18856)				   # dim(z)   = 206 18856
z<-t(z)						   # dim(z)   = 18856 206
w<-z[,2:206]					   # dim(w)   = 18856 205
mz<-z[,1]

p<-dim(w)[1]
n<-dim(w)[2]

ww<-w[c(1:p)[mz>=2000],]

###################################
#normalize total ion count

ref<-sum(ww)/205

wwsave<-ww

for (i in 1:205){
ww[,i]<-ww[,i]/sum(ww[,i])*ref
}
######################################

## drop m/z values < 2000 as they are too noisy

mzchop<-mz[mz>=2000]
pww<-dim(ww)[1]
image(-log(ww+3),xlab=" ",ylab=" ",axes=FALSE)


mark<-rep(0,times=5)
for (ii in 1:5){
mark[ii]<-min(c(1:pww)[mzchop>5000*ii])
}
mn<--0.1
mx<-1.1
for (i in 1:5){
lines(c(mark[i]/pww,mark[i]/pww),c(mn,mx),lty=2)
}



  image(mzchop,1:205,-log(ww+3),ylab="Sample number",xlab="m/z")
 lines(c(0,30000),c(101.5,101.5))
 text(28000,10,"Stage 1")
 text(28000,190,"Stage 4")


mzchopsave<-mzchop
wwsave<-ww




#######################################################################
dnorm2<- function( xx , mu, xi, delta, omega, pi){
dens<- pi*dnorm( xx , mu, sqrt(xi)*mu) + (1-pi)*dnorm( xx, mu*delta, omega*sqrt(xi)*mu) 
dens
}


detpeakfit2<-function(xi,delta,omega, pi,npeaks){

peakloc<- rep(0,times=npeaks)
peakhgt<- rep(0,times=npeaks)


if (npeaks==0){

y<-apply(ww,1,mean)

yhat<- mzchop*0

}


if (npeaks>0){

sumsq<-rep(0,times=npeaks)
loglike<-sumsq

y<-apply(ww,1,mean)
#plot(mzchop,y,type="l")
iloc<- c(1:pww)[y==max(y)]
mu<- mzchop[iloc]
yhat<- dnorm2( mzchop, mu, xi, delta, omega, pi)*y[iloc]/dnorm2( mu, mu, xi, delta, omega, pi) 

#lines(mzchop,yhat,col=2)


sumsq[1]<- sum((y-yhat)**(2)) 

if (chunk==0){
loglike[1] <- arima(  y-yhat , c(2,0,0) )$loglik
}
peakloc[1]<-iloc
peakhgt[1]<-y[iloc]


i<-1
print("npeaks")
print(npeaks)


if (npeaks>1){



for (i in 2:npeaks){
cat(c(" ",i))
resid <- (y - yhat)/(1+yhat)
peakloc[i] <- c(1:pww)[resid==max(resid)]
peakhgt[i] <- y[ peakloc[i] ]

Fmat<-matrix(0,i,i)

for ( j1 in 1:i ){
for ( j2 in 1:i){

mu1<- mzchop[peakloc[j1]]
mu2<- mzchop[peakloc[j2]]

Fmat[j1,j2] <- dnorm2( mu1 , mu2 , xi, delta, omega , pi)

}
}

cc <- solve(Fmat)%*%peakhgt[1:i]

sum<-0
for (j in 1:i){
mu<-mzchop[peakloc[j]]
sum <- sum + cc[j]*dnorm2( mzchop , mu , xi, delta, omega, pi)
}
yhat<-sum
sumsq[i]<- sum((y-yhat)**(2)) 
if (chunk==0){
loglike[i] <- arima(  y-yhat , c(2,0,0) )$loglik
}
#print(loglike[i])
#plot(acf(y-yhat))
#plot(pacf(y-yhat))
}
}
}


plot(mzchop,y,type="l",ylim=c(0,max(y)),ylab="intensity",xlab="m/z")
lines(mzchop,yhat,col=2)
points(mzchop[peakloc[1:i]],rep(0,times=i),pch=20)
#print("GOF")
GOF<-sum((y-yhat)**(2)) 
#print( sum((y-yhat)**(2)) )
out<-list(peakloc=0,y=0,yhat=0,GOF=0,ss=0)
out$peakloc<-peakloc
out$y<-y
out$yhat<-yhat
out$GOF<-GOF
out$ss<-sumsq
if (chunk==0){
out$loglike <- loglike
}
out
}

##===================================================================
}
##################################################################


if (chunk > 0){


## the splits have been obtained previously using a C program

splits<-c(737,942,659,982,1175,659,804,  415+349+1023+985+328+898 , 
          317+493+1108+368, 658+1051)
csplits<-cumsum(splits)
ind<-c((csplits[chunk]-splits[chunk]+1):csplits[chunk])


##combine original splits (normalized)
## npeaks using total 200 peaks using well-tuned good-looking parameters
peakvec<-c(34,24,17,20,29,15,18, 22,16,5)
 

        
npeak<-peakvec[chunk]
print(peakvec)
}

if (chunk==0){
ind<-c(1:length(mzchopsave))
npeak<-500 
}

mzchop<- mzchopsave[ind]
ww<-wwsave[ind,]
pww<-dim(ww)[1]


### well-tuned good-looking parameters (from approx median of subset splitting auto-fits)
delta<-1.006
omega<-2
pi<-0.5
xi<-0.8e-05


ans<-detpeakfit2(xi,delta,omega,pi,npeak)

func<-function(x){
ans<-detpeakfit2(abs(x[1]/1000000),1+x[2]/10000,x[3]/10,x[4]/10,npeak)
print(x)
ans$GOF
}

#print("Do you want to autofit the parameters?")
#ans3<-readline()

ans3<-"n"


if (ans3=="y"){

#iterative fit to find better parameters
par(mfrow=c(1,1))
autofit<-nlm( func, p=c(10,40,20,5) )
x<-autofit$estimate
xi<-abs(x[1]/1000000)
delta<-1+x[2]/10000
omega<-x[3]/10
pi<-x[4]/10

print("parameters")
print(c(xi,delta,omega,pi))

}
 

par(mfrow=c(2,2))


plot(mzchop,ww[,1],type="n",ylim=c(0,max(ww)),ylab="intensity",xlab="m/z")
for (i in 1:101){
lines( mzchop, ww[,i], col=2)
}
for (i in 102:205){
lines( mzchop, ww[,i], lty=2, col=3)
}

y1<-apply(ww[,1:101],1,mean)
y2<-apply(ww[,102:205],1,mean)

plot(mzchop,y1,type="n",ylim=c(0,max(ww)),ylab="intensity",xlab="m/z")
lines(mzchop,y1,col=2)
lines(mzchop,y2,col=3)
s1<-sqrt(diag(var(t(ww[,1:101]))))
s2<-sqrt(diag(var(t(ww[,102:205]))))

plot(mzchop,y1,type="n",ylim=c(-max(c(y1,y2)),max(c(y1,y2))),ylab="CI",xlab="m/z")
lines(mzchop,y1-y2+2*(s1/sqrt(101)+s2/sqrt(104)),col=4)
lines(mzchop,y1-y2-2*(s1/sqrt(101)+s2/sqrt(104)),col=4)
lines(mzchop,mzchop*0,col=5)

ans<-detpeakfit2(xi,delta,omega,pi,npeak)


####################################################################

y<-ans$y
yhat<-ans$yhat



nsplit<-length(ind)


mu<-sort(ans$peakloc)

print("no of peaks here")
print(length(mu))


wsub<-ww
response<-c(wsub)
peak <-factor( rep(1:nsplit,times=205) )
patient <- factor( rep( c(1:205), each=nsplit ) )
group <- factor ( c( rep("I",times=nsplit*101), rep("IV",times=nsplit*104) ) )

# basis vectors

npeak<-length(mu)
basismat<-matrix(0,nsplit,npeak)
Zmat1<-matrix(0,205*nsplit,npeak)
Zmat2<-matrix(0,205*nsplit,npeak)
for (i in 1:npeak){
basismat[,i]<-dnorm2( mzchop , mzchop[mu[i]] , xi, delta, omega, pi)
Zmat1[,i]<-c(rep(basismat[,i],times=101),rep(basismat[,i],times=104))
Zmat2[,i]<-c(rep(0,times=101*nsplit),rep(basismat[,i],times=104))
}


mel<-data.frame( response, Zmat1, Zmat2, patient, group)
print("fitting fixed effects model")
outlm<- lm( response ~ 0+Zmat1+Zmat2, data=mel)

library(lme4)

print("Fitting linear mixed effects model")

############
#********************
# need to adjust the number of random effects
# to give the correct no of peaks for independent 
#  random effects
##**************************************************

if (npeak==35){
out4<- lmer( response ~ 0+Zmat1[,1:npeak]+Zmat2[,1:npeak] + 
    (0+Zmat1[,1] | patient ) 
  + (0+Zmat1[,2] | patient ) 
  + (0+Zmat1[,3] | patient ) 
  + (0+Zmat1[,4] | patient ) 
  + (0+Zmat1[,5] | patient ) 
  + (0+Zmat1[,6] | patient ) 
  + (0+Zmat1[,7] | patient ) 
  + (0+Zmat1[,8] | patient ) 
  + (0+Zmat1[,9] | patient ) 
  + (0+Zmat1[,10] | patient ) 
  + (0+Zmat1[,11] | patient ) 
  + (0+Zmat1[,12] | patient ) 
  + (0+Zmat1[,13] | patient ) 
  + (0+Zmat1[,14] | patient ) 
  + (0+Zmat1[,15] | patient ) 
  + (0+Zmat1[,16] | patient ) 
  + (0+Zmat1[,17] | patient ) 
  + (0+Zmat1[,18] | patient ) 
  + (0+Zmat1[,19] | patient ) 
  + (0+Zmat1[,20] | patient ) 
  + (0+Zmat1[,21] | patient ) 
  + (0+Zmat1[,22] | patient ) 
  + (0+Zmat1[,23] | patient ) 
  + (0+Zmat1[,24] | patient ) 
  + (0+Zmat1[,25] | patient ) 
  + (0+Zmat1[,26] | patient ) 
  + (0+Zmat1[,27] | patient ) 
  + (0+Zmat1[,28] | patient ) 
  + (0+Zmat1[,29] | patient ) 
  + (0+Zmat1[,30] | patient )  
 + (0+Zmat1[,31] | patient ) 
  + (0+Zmat1[,32] | patient ) 
  + (0+Zmat1[,33] | patient ) 
  + (0+Zmat1[,34] | patient ) 
  + (0+Zmat1[,35] | patient ) 
 , data = mel)
}

if (npeak==34){
out4<- lmer( response ~ 0+Zmat1[,1:npeak]+Zmat2[,1:npeak] + 
    (0+Zmat1[,1] | patient ) 
  + (0+Zmat1[,2] | patient ) 
  + (0+Zmat1[,3] | patient ) 
  + (0+Zmat1[,4] | patient ) 
  + (0+Zmat1[,5] | patient ) 
  + (0+Zmat1[,6] | patient ) 
  + (0+Zmat1[,7] | patient ) 
  + (0+Zmat1[,8] | patient ) 
  + (0+Zmat1[,9] | patient ) 
  + (0+Zmat1[,10] | patient ) 
  + (0+Zmat1[,11] | patient ) 
  + (0+Zmat1[,12] | patient ) 
  + (0+Zmat1[,13] | patient ) 
  + (0+Zmat1[,14] | patient ) 
  + (0+Zmat1[,15] | patient ) 
  + (0+Zmat1[,16] | patient ) 
  + (0+Zmat1[,17] | patient ) 
  + (0+Zmat1[,18] | patient ) 
  + (0+Zmat1[,19] | patient ) 
  + (0+Zmat1[,20] | patient ) 
  + (0+Zmat1[,21] | patient ) 
  + (0+Zmat1[,22] | patient ) 
  + (0+Zmat1[,23] | patient ) 
  + (0+Zmat1[,24] | patient ) 
  + (0+Zmat1[,25] | patient ) 
  + (0+Zmat1[,26] | patient ) 
  + (0+Zmat1[,27] | patient ) 
  + (0+Zmat1[,28] | patient ) 
  + (0+Zmat1[,29] | patient ) 
  + (0+Zmat1[,30] | patient )  
 + (0+Zmat1[,31] | patient ) 
  + (0+Zmat1[,32] | patient ) 
  + (0+Zmat1[,33] | patient ) 
  + (0+Zmat1[,34] | patient ) 
 , data = mel)
}


if (npeak==30){
out4<- lmer( response ~ 0+Zmat1[,1:npeak]+Zmat2[,1:npeak] + 
    (0+Zmat1[,1] | patient ) 
  + (0+Zmat1[,2] | patient ) 
  + (0+Zmat1[,3] | patient ) 
  + (0+Zmat1[,4] | patient ) 
  + (0+Zmat1[,5] | patient ) 
  + (0+Zmat1[,6] | patient ) 
  + (0+Zmat1[,7] | patient ) 
  + (0+Zmat1[,8] | patient ) 
  + (0+Zmat1[,9] | patient ) 
  + (0+Zmat1[,10] | patient ) 
  + (0+Zmat1[,11] | patient ) 
  + (0+Zmat1[,12] | patient ) 
  + (0+Zmat1[,13] | patient ) 
  + (0+Zmat1[,14] | patient ) 
  + (0+Zmat1[,15] | patient ) 
  + (0+Zmat1[,16] | patient ) 
  + (0+Zmat1[,17] | patient ) 
  + (0+Zmat1[,18] | patient ) 
  + (0+Zmat1[,19] | patient ) 
  + (0+Zmat1[,20] | patient ) 
  + (0+Zmat1[,21] | patient ) 
  + (0+Zmat1[,22] | patient ) 
  + (0+Zmat1[,23] | patient ) 
  + (0+Zmat1[,24] | patient ) 
  + (0+Zmat1[,25] | patient ) 
  + (0+Zmat1[,26] | patient ) 
  + (0+Zmat1[,27] | patient ) 
  + (0+Zmat1[,28] | patient ) 
  + (0+Zmat1[,29] | patient ) 
  + (0+Zmat1[,30] | patient )  , data = mel)
}

if (npeak==24){
out4<- lmer( response ~ 0+Zmat1[,1:npeak]+Zmat2[,1:npeak] + 
    (0+Zmat1[,1] | patient ) 
  + (0+Zmat1[,2] | patient ) 
  + (0+Zmat1[,3] | patient ) 
  + (0+Zmat1[,4] | patient ) 
  + (0+Zmat1[,5] | patient ) 
  + (0+Zmat1[,6] | patient ) 
  + (0+Zmat1[,7] | patient ) 
  + (0+Zmat1[,8] | patient ) 
  + (0+Zmat1[,9] | patient ) 
  + (0+Zmat1[,10] | patient ) 
  + (0+Zmat1[,11] | patient ) 
  + (0+Zmat1[,12] | patient ) 
  + (0+Zmat1[,13] | patient ) 
  + (0+Zmat1[,14] | patient ) 
  + (0+Zmat1[,15] | patient ) 
  + (0+Zmat1[,16] | patient ) 
  + (0+Zmat1[,17] | patient ) 
  + (0+Zmat1[,18] | patient ) 
  + (0+Zmat1[,19] | patient ) 
  + (0+Zmat1[,20] | patient ) 
  + (0+Zmat1[,21] | patient ) 
  + (0+Zmat1[,22] | patient ) 
  + (0+Zmat1[,23] | patient ) 
  + (0+Zmat1[,24] | patient ) 
 , data = mel)
}

if (npeak==29){
out4<- lmer( response ~ 0+Zmat1[,1:npeak]+Zmat2[,1:npeak] + 
    (0+Zmat1[,1] | patient ) 
  + (0+Zmat1[,2] | patient ) 
  + (0+Zmat1[,3] | patient ) 
  + (0+Zmat1[,4] | patient ) 
  + (0+Zmat1[,5] | patient ) 
  + (0+Zmat1[,6] | patient ) 
  + (0+Zmat1[,7] | patient ) 
  + (0+Zmat1[,8] | patient ) 
  + (0+Zmat1[,9] | patient ) 
  + (0+Zmat1[,10] | patient ) 
  + (0+Zmat1[,11] | patient ) 
  + (0+Zmat1[,12] | patient ) 
  + (0+Zmat1[,13] | patient ) 
  + (0+Zmat1[,14] | patient ) 
  + (0+Zmat1[,15] | patient ) 
  + (0+Zmat1[,16] | patient ) 
  + (0+Zmat1[,17] | patient ) 
  + (0+Zmat1[,18] | patient ) 
  + (0+Zmat1[,19] | patient ) 
  + (0+Zmat1[,20] | patient ) 
  + (0+Zmat1[,21] | patient ) 
  + (0+Zmat1[,22] | patient ) 
  + (0+Zmat1[,23] | patient ) 
  + (0+Zmat1[,24] | patient ) 
  + (0+Zmat1[,25] | patient ) 
  + (0+Zmat1[,26] | patient ) 
  + (0+Zmat1[,27] | patient ) 
  + (0+Zmat1[,28] | patient ) 
  + (0+Zmat1[,29] | patient )   , data = mel)
}
if (npeak == 20){
out4<- lmer( response ~ 0+Zmat1[,1:npeak]+Zmat2[,1:npeak] + 
    (0+Zmat1[,1] | patient ) 
  + (0+Zmat1[,2] | patient ) 
  + (0+Zmat1[,3] | patient ) 
  + (0+Zmat1[,4] | patient ) 
  + (0+Zmat1[,5] | patient ) 
  + (0+Zmat1[,6] | patient ) 
  + (0+Zmat1[,7] | patient ) 
  + (0+Zmat1[,8] | patient ) 
  + (0+Zmat1[,9] | patient ) 
  + (0+Zmat1[,10] | patient ) 
  + (0+Zmat1[,11] | patient ) 
  + (0+Zmat1[,12] | patient ) 
  + (0+Zmat1[,13] | patient ) 
  + (0+Zmat1[,14] | patient ) 
  + (0+Zmat1[,15] | patient ) 
  + (0+Zmat1[,16] | patient ) 
  + (0+Zmat1[,17] | patient ) 
  + (0+Zmat1[,18] | patient ) 
  + (0+Zmat1[,19] | patient ) 
  + (0+Zmat1[,20] | patient )  , data = mel)
}
if (npeak == 10){
out4<- lmer( response ~ 0+Zmat1[,1:npeak]+Zmat2[,1:npeak] + 
    (0+Zmat1[,1] | patient ) 
  + (0+Zmat1[,2] | patient ) 
  + (0+Zmat1[,3] | patient ) 
  + (0+Zmat1[,4] | patient ) 
  + (0+Zmat1[,5] | patient ) 
  + (0+Zmat1[,6] | patient ) 
  + (0+Zmat1[,7] | patient ) 
  + (0+Zmat1[,8] | patient ) 
  + (0+Zmat1[,9] | patient ) 
  + (0+Zmat1[,10] | patient ) 
 , data = mel)
}
if (npeak == 5){
out4<- lmer( response ~ 0+Zmat1[,1:npeak]+Zmat2[,1:npeak] + 
    (0+Zmat1[,1] | patient ) 
  + (0+Zmat1[,2] | patient ) 
  + (0+Zmat1[,3] | patient ) 
  + (0+Zmat1[,4] | patient ) 
  + (0+Zmat1[,5] | patient ) 
 , data = mel)
}
if (npeak == 7){
out4<- lmer( response ~ 0+Zmat1[,1:npeak]+Zmat2[,1:npeak] + 
    (0+Zmat1[,1] | patient ) 
  + (0+Zmat1[,2] | patient ) 
  + (0+Zmat1[,3] | patient ) 
  + (0+Zmat1[,4] | patient ) 
  + (0+Zmat1[,5] | patient ) 
  + (0+Zmat1[,6] | patient ) 
  + (0+Zmat1[,7] | patient ) 
 , data = mel)
}

if (npeak==17){
out4<- lmer( response ~ 0+Zmat1[,1:npeak]+Zmat2[,1:npeak] + 
    (0+Zmat1[,1] | patient ) 
  + (0+Zmat1[,2] | patient ) 
  + (0+Zmat1[,3] | patient ) 
  + (0+Zmat1[,4] | patient ) 
  + (0+Zmat1[,5] | patient ) 
  + (0+Zmat1[,6] | patient ) 
  + (0+Zmat1[,7] | patient ) 
  + (0+Zmat1[,8] | patient ) 
  + (0+Zmat1[,9] | patient ) 
  + (0+Zmat1[,10] | patient ) 
  + (0+Zmat1[,11] | patient ) 
  + (0+Zmat1[,12] | patient ) 
  + (0+Zmat1[,13] | patient ) 
  + (0+Zmat1[,14] | patient ) 
  + (0+Zmat1[,15] | patient ) 
  + (0+Zmat1[,16] | patient ) 
  + (0+Zmat1[,17] | patient ) 
 , data = mel)
}

if (npeak==16){
out4<- lmer( response ~ 0+Zmat1[,1:npeak]+Zmat2[,1:npeak] + 
    (0+Zmat1[,1] | patient ) 
  + (0+Zmat1[,2] | patient ) 
  + (0+Zmat1[,3] | patient ) 
  + (0+Zmat1[,4] | patient ) 
  + (0+Zmat1[,5] | patient ) 
  + (0+Zmat1[,6] | patient ) 
  + (0+Zmat1[,7] | patient ) 
  + (0+Zmat1[,8] | patient ) 
  + (0+Zmat1[,9] | patient ) 
  + (0+Zmat1[,10] | patient ) 
  + (0+Zmat1[,11] | patient ) 
  + (0+Zmat1[,12] | patient ) 
  + (0+Zmat1[,13] | patient ) 
  + (0+Zmat1[,14] | patient ) 
  + (0+Zmat1[,15] | patient ) 
  + (0+Zmat1[,16] | patient ) 
 , data = mel)
}
if (npeak==28){
out4<- lmer( response ~ 0+Zmat1[,1:npeak]+Zmat2[,1:npeak] + 
    (0+Zmat1[,1] | patient ) 
  + (0+Zmat1[,2] | patient ) 
  + (0+Zmat1[,3] | patient ) 
  + (0+Zmat1[,4] | patient ) 
  + (0+Zmat1[,5] | patient ) 
  + (0+Zmat1[,6] | patient ) 
  + (0+Zmat1[,7] | patient ) 
  + (0+Zmat1[,8] | patient ) 
  + (0+Zmat1[,9] | patient ) 
  + (0+Zmat1[,10] | patient ) 
  + (0+Zmat1[,11] | patient ) 
  + (0+Zmat1[,12] | patient ) 
  + (0+Zmat1[,13] | patient ) 
  + (0+Zmat1[,14] | patient ) 
  + (0+Zmat1[,15] | patient ) 
  + (0+Zmat1[,16] | patient ) 
  + (0+Zmat1[,17] | patient ) 
  + (0+Zmat1[,18] | patient ) 
  + (0+Zmat1[,19] | patient ) 
  + (0+Zmat1[,20] | patient ) 
  + (0+Zmat1[,21] | patient ) 
  + (0+Zmat1[,22] | patient ) 
  + (0+Zmat1[,23] | patient ) 
  + (0+Zmat1[,24] | patient ) 
  + (0+Zmat1[,25] | patient ) 
  + (0+Zmat1[,26] | patient ) 
  + (0+Zmat1[,27] | patient ) 
  + (0+Zmat1[,28] | patient ) 
   , data = mel)
}
if (npeak==16){
out4<- lmer( response ~ 0+Zmat1[,1:npeak]+Zmat2[,1:npeak] + 
    (0+Zmat1[,1] | patient ) 
  + (0+Zmat1[,2] | patient ) 
  + (0+Zmat1[,3] | patient ) 
  + (0+Zmat1[,4] | patient ) 
  + (0+Zmat1[,5] | patient ) 
  + (0+Zmat1[,6] | patient ) 
  + (0+Zmat1[,7] | patient ) 
  + (0+Zmat1[,8] | patient ) 
  + (0+Zmat1[,9] | patient ) 
  + (0+Zmat1[,10] | patient ) 
  + (0+Zmat1[,11] | patient ) 
  + (0+Zmat1[,12] | patient ) 
  + (0+Zmat1[,13] | patient ) 
  + (0+Zmat1[,14] | patient ) 
  + (0+Zmat1[,15] | patient ) 
  + (0+Zmat1[,16] | patient ) 
 , data = mel)
}
if (npeak==18){
out4<- lmer( response ~ 0+Zmat1[,1:npeak]+Zmat2[,1:npeak] + 
    (0+Zmat1[,1] | patient ) 
  + (0+Zmat1[,2] | patient ) 
  + (0+Zmat1[,3] | patient ) 
  + (0+Zmat1[,4] | patient ) 
  + (0+Zmat1[,5] | patient ) 
  + (0+Zmat1[,6] | patient ) 
  + (0+Zmat1[,7] | patient ) 
  + (0+Zmat1[,8] | patient ) 
  + (0+Zmat1[,9] | patient ) 
  + (0+Zmat1[,10] | patient ) 
  + (0+Zmat1[,11] | patient ) 
  + (0+Zmat1[,12] | patient ) 
  + (0+Zmat1[,13] | patient ) 
  + (0+Zmat1[,14] | patient ) 
  + (0+Zmat1[,15] | patient ) 
  + (0+Zmat1[,16] | patient ) 
  + (0+Zmat1[,17] | patient ) 
  + (0+Zmat1[,18] | patient ) 
 , data = mel)
}

if (npeak==15){
out4<- lmer( response ~ 0+Zmat1[,1:npeak]+Zmat2[,1:npeak] + 
    (0+Zmat1[,1] | patient ) 
  + (0+Zmat1[,2] | patient ) 
  + (0+Zmat1[,3] | patient ) 
  + (0+Zmat1[,4] | patient ) 
  + (0+Zmat1[,5] | patient ) 
  + (0+Zmat1[,6] | patient ) 
  + (0+Zmat1[,7] | patient ) 
  + (0+Zmat1[,8] | patient ) 
  + (0+Zmat1[,9] | patient ) 
  + (0+Zmat1[,10] | patient ) 
  + (0+Zmat1[,11] | patient ) 
  + (0+Zmat1[,12] | patient ) 
  + (0+Zmat1[,13] | patient ) 
  + (0+Zmat1[,14] | patient ) 
  + (0+Zmat1[,15] | patient ) 
 , data = mel)
}



if (npeak==6){
out4<- lmer( response ~ 0+Zmat1[,1:npeak]+Zmat2[,1:npeak] + 
    (0+Zmat1[,1] | patient ) 
  + (0+Zmat1[,2] | patient ) 
  + (0+Zmat1[,3] | patient ) 
  + (0+Zmat1[,4] | patient ) 
  + (0+Zmat1[,5] | patient ) 
  + (0+Zmat1[,6] | patient ) 
 , data = mel)
}


if (npeak==2){
out4<- lmer( response ~ 0+Zmat1[,1:npeak]+Zmat2[,1:npeak] + 
    (0+Zmat1[,1] | patient ) 
  + (0+Zmat1[,2] | patient ) 
   , data = mel)
}
if (npeak==9){
out4<- lmer( response ~ 0+Zmat1[,1:npeak]+Zmat2[,1:npeak] + 
    (0+Zmat1[,1] | patient ) 
  + (0+Zmat1[,2] | patient ) 
  + (0+Zmat1[,3] | patient ) 
  + (0+Zmat1[,4] | patient ) 
  + (0+Zmat1[,5] | patient ) 
  + (0+Zmat1[,6] | patient ) 
  + (0+Zmat1[,7] | patient ) 
  + (0+Zmat1[,8] | patient ) 
  + (0+Zmat1[,9] | patient ) 
  , data = mel)
}
if (npeak==22){
out4<- lmer( response ~ 0+Zmat1[,1:npeak]+Zmat2[,1:npeak] + 
    (0+Zmat1[,1] | patient ) 
  + (0+Zmat1[,2] | patient ) 
  + (0+Zmat1[,3] | patient ) 
  + (0+Zmat1[,4] | patient ) 
  + (0+Zmat1[,5] | patient ) 
  + (0+Zmat1[,6] | patient ) 
  + (0+Zmat1[,7] | patient ) 
  + (0+Zmat1[,8] | patient ) 
  + (0+Zmat1[,9] | patient ) 
  + (0+Zmat1[,10] | patient ) 
  + (0+Zmat1[,11] | patient ) 
  + (0+Zmat1[,12] | patient ) 
  + (0+Zmat1[,13] | patient ) 
  + (0+Zmat1[,14] | patient ) 
  + (0+Zmat1[,15] | patient ) 
  + (0+Zmat1[,16] | patient ) 
  + (0+Zmat1[,17] | patient ) 
  + (0+Zmat1[,18] | patient ) 
  + (0+Zmat1[,19] | patient ) 
  + (0+Zmat1[,20] | patient ) 
  + (0+Zmat1[,21] | patient ) 
  + (0+Zmat1[,22] | patient ) 
   , data = mel)
}
if (npeak==4){
out4<- lmer( response ~ 0+Zmat1[,1:npeak]+Zmat2[,1:npeak] + 
    (0+Zmat1[,1] | patient ) 
  + (0+Zmat1[,2] | patient ) 
  + (0+Zmat1[,3] | patient ) 
  + (0+Zmat1[,4] | patient ) 
  , data = mel)
}
if (npeak==15){
out4<- lmer( response ~ 0+Zmat1[,1:npeak]+Zmat2[,1:npeak] + 
    (0+Zmat1[,1] | patient ) 
  + (0+Zmat1[,2] | patient ) 
  + (0+Zmat1[,3] | patient ) 
  + (0+Zmat1[,4] | patient ) 
  + (0+Zmat1[,5] | patient ) 
  + (0+Zmat1[,6] | patient ) 
  + (0+Zmat1[,7] | patient ) 
  + (0+Zmat1[,8] | patient ) 
  + (0+Zmat1[,9] | patient ) 
  + (0+Zmat1[,10] | patient ) 
  + (0+Zmat1[,11] | patient ) 
  + (0+Zmat1[,12] | patient ) 
  + (0+Zmat1[,13] | patient ) 
  + (0+Zmat1[,14] | patient ) 
  + (0+Zmat1[,15] | patient ) 
   , data = mel)
}
if (npeak==1){
out4<- lmer( response ~ 0+Zmat1[,1:npeak]+Zmat2[,1:npeak] + 
    (0+Zmat1[,1] | patient ) 
   , data = mel)
}



print(summary(outlm))
print(summary(out4))


beta<-fixef(out4)
Vcov <- vcov(out4, useScale = FALSE)
sebeta <- sqrt(diag(Vcov))
tstat<-c(beta/sebeta)[(npeak+1):(2*npeak)]

# from FDR=0.01 we have threshold 3.65023 with normalized data

points(mzchop[sort(ans$peakloc)][abs(tstat)>3.65023],mzchop[sort(ans$peakloc)][abs(tstat)>3.65023]*0,col=5,cex=2)

dev.print(file="chunk-results.ps")



print("Our mixed effects model results")
print(cbind(mzchop[mu],as.double(tstat)))

print("Compare with: standard t-stat at peaks")


# standard t - test
standardt<-  -(y1-y2) /  sqrt(s1**2/101 + s2**2/104) 

print(cbind(mzchop[mu], standardt[mu] ) )
plot(standardt[mu],as.double(tstat))
lines(c(-33,33),c(-33,33))