SCVA
Single-Case Visual Analysis

Package index
Search the SCVA package
Vignettes
Functions
14
Source code
6
Man pages
8
Home
/
CRAN
/
SCVA
/
R/graph.TREND.R

R/graph.TREND.R
In SCVA: Single-Case Visual Analysis

Defines functions graph.TREND

Documented in graph.TREND 

graph.TREND <-
function(design,TREND,CL,tr,data=read.table(file.choose(new=FALSE)),xlab=NULL,ylab=NULL,ylim=NULL,legendxy=NULL,labels=c("A","B","A","B")){
  
  x<-1:nrow(data)
  MT<-nrow(data)
  
  if(design=="AB"|design=="CRD"|design=="ATD"|design=="RBD"|design=="Custom"){
    A<-data[,2][data[,1]=="A"]
    B<-data[,2][data[,1]=="B"]
    MTA<-x[data[,1]=="A"]
    MTB<-x[data[,1]=="B"]
    if(TREND=="VLP"){
      if(CL=="mean"){
        CLA<-mean(A,na.rm=TRUE)
        CLB<-mean(B,na.rm=TRUE)
      }
      if(CL=="median"){
        CLA<-median(A,na.rm=TRUE)
        CLB<-median(B,na.rm=TRUE)
      }
      if(CL=="bmed"){
        aa<-sort(A)
        bb<-sort(B)
        if(length(aa)<5){
          CLA<-median(A,na.rm=TRUE)
        }
        if(length(aa)==5|length(aa)==7|length(aa)==9|length(aa)==11){
          CLA<-(aa[ceiling(length(aa)/2)-1]+aa[ceiling(length(aa)/2)]+aa[ceiling(length(aa)/2)+1])/3
        }
        if(length(aa)>=13&length(aa)%%2==1){
          CLA<-(aa[ceiling(length(aa)/2)-2]+aa[ceiling(length(aa)/2)-1]+aa[ceiling(length(aa)/2)]+aa[ceiling(length(aa)/2)+1] +aa[ceiling(length(aa)/2)+2])/5
        }
        if(length(aa)==6|length(aa)==8|length(aa)==10|length(aa)==12){
          CLA<-1/6*aa[length(aa)/2-1]+1/3*aa[length(aa)/2]+1/3*aa[length(aa)/2+1]+1/6*aa[length(aa)/2+2]
        }
        if(length(aa)>13&length(aa)%%2==0){
          CLA<-1/10*aa[length(aa)/2-2]+1/5*aa[length(aa)/2-1]+1/5*aa[length(aa)/2]+1/5*aa[length(aa)/2+1]+1/5*aa[length(aa)/2+2]+1/10*aa[length(aa)/2+3]
        }
        if(length(bb)<5){
          CLB<-median(B,na.rm=TRUE)
        }
        if(length(bb)==5|length(bb)==7|length(bb)==9|length(bb)==11){
          CLB<-(bb[ceiling(length(bb)/2)-1]+bb[ceiling(length(bb)/2)]+bb[ceiling(length(bb)/2)+1])/3
        }
        if(length(bb)>=13&length(bb)%%2==1){
          CLB<-(bb[ceiling(length(bb)/2)-2]+bb[ceiling(length(bb)/2)-1]+bb[ceiling(length(bb)/2)]+bb[ceiling(length(bb)/2)+1]+bb[ceiling(length(bb)/2)+2])/5
        }
        if(length(bb)==6|length(bb)==8|length(bb)==10|length(bb)==12){
          CLB<-1/6*bb[length(bb)/2-1]+1/3*bb[length(bb)/2]+1/3*bb[length(bb)/2+1]+1/6*bb[length(bb)/2+2]
        }
        if(length(bb)>13&length(bb)%%2==0){
          CLB<-1/10*bb[length(bb)/2-2]+1/5*bb[length(bb)/2-1]+1/5*bb[length(bb)/2]+1/5*bb[length(bb)/2+1]+1/5*bb[length(bb)/2+2]+1/10*bb[length(bb)/2+3]
        }
      } 
      if(CL=="trimmean"){
        CLA<-mean(A,trim=tr,na.rm=TRUE)	
        CLB<-mean(B,trim=tr,na.rm=TRUE)
      }  
      if(CL=="mest"){
        hpsi<-function(x,bend=1.28){
          hpsi<-ifelse(abs(x)<=bend,x,bend*sign(x))
          hpsi
        }
        mest<-function(x,bend=1.28,na.rm=TRUE){
          if(na.rm)x<-x[!is.na(x)]
          if(length(x)==0) return(NA)
          if(mad(x)==0)stop("MAD=0. The M-estimator cannot be computed.")
          y<-(x-median(x))/mad(x)
          A<-sum(hpsi(y,bend))
          B<-length(x[abs(y)<=bend])
          mest<-median(x)+mad(x)*A/B
          repeat{
            y<-(x-mest)/mad(x)
            A<-sum(hpsi(y,bend))
            B<-length(x[abs(y)<=bend])
            newmest<-mest+mad(x)*A/B
            if(abs(newmest-mest) <.0001)break
            mest<-newmest
          }
          mest
        }
        CLA<-mest(A,bend=tr)	
        CLB<-mest(B,bend=tr)
      } 
      residualsA<-A-CLA
      residualsB<-B-CLB
      residuals<-c(residualsA,residualsB)
      if (is.null(xlab)){
        xlab <- "Measurement Times"
      }
      if (is.null(ylab)){
        ylab <- "Residuals"
      }
      plot(x,residuals,xlab=xlab,ylab=ylab,ylim=ylim,type="n")
      lines(c(1,length(residualsA)),c(0,0))
      for(it in 1:length(residualsA)){
        lines(c(it,it),c(residualsA[it],0))
      }
      lines(c(length(residualsA)+1,length(residuals)),c(0,0))
      for(it in 1:length(residualsB)){
        lines(c(length(residualsA)+it,length(residualsA)+it),c(residualsB[it],0))
      }
      mtext(labels[1],side=3,at=(sum(data[,1]=="A")+1)/2)
      mtext(labels[2],side=3,at=(sum(data[,1]=="A")+(sum(data[,1]=="B")+1)/2))
    }
    if(TREND=="RTL"|TREND=="SM"|TREND=="LSR"|TREND=="RM3"|TREND=="RM5"|TREND=="RM42"){
      if (is.null(xlab)){
        xlab <- "Measurement Times"
      }
      if (is.null(ylab)){
        ylab <- "Scores"
      }
      if(design=="AB"){
        plot(x,data[,2],xlab=xlab,ylab=ylab,ylim=ylim,pch=16)
        lines(c(sum(data[,1]=="A")+0.5,sum(data[,1]=="A")+0.5),c(min(data[,2],ylim[1],na.rm=TRUE)-5,max(data[,2],ylim[2],na.rm=TRUE)+5),lty=2)
        mtext(labels[1],side=3,at=(sum(data[,1]=="A")+1)/2)
        mtext(labels[2],side=3,at=(sum(data[,1]=="A")+(sum(data[,1]=="B")+1)/2))
      }
      if(design=="CRD"|design=="ATD"|design=="RBD"|design=="Custom"){
        plot(x,data[,2],type="n",xlab=xlab,ylab=ylab,ylim=ylim)
        points(x[data[,1]=="A"],data[,2][data[,1]=="A"],pch=1)
        points(x[data[,1]=="B"],data[,2][data[,1]=="B"],pch=16)
      }
      if(TREND=="LSR"){
        interceptLSRa<-coefficients(lm(A~MTA,na.action=na.omit))[1]
        slopeLSRa<-coefficients(lm(A~MTA,na.action=na.omit))[2]
        interceptLSRb<-coefficients(lm(B~MTB,na.action=na.omit))[1]
        slopeLSRb<-coefficients(lm(B~MTB,na.action=na.omit))[2]
        if(design=="AB"){
          lines(c(1,length(A)),c(interceptLSRa+slopeLSRa,interceptLSRa+slopeLSRa*length(A)),lty=2)
          lines(c(length(A)+1,MT),c(interceptLSRb+slopeLSRb*(length(A)+1),interceptLSRb+slopeLSRb*MT),lty=2)
        }
        if(design=="CRD"|design=="ATD"|design=="RBD"|design=="Custom"){
          lines(c(1,MT),c(interceptLSRa+slopeLSRa,interceptLSRa+slopeLSRa*MT),lty=1)
          lines(c(1,MT),c(interceptLSRb+slopeLSRb,interceptLSRb+slopeLSRb*MT),lty=16)
          trend_labels<-c(labels[1:2],paste("trend line (LSR)",labels[1:2]))
          if(is.null(legendxy))
            legend(locator(1),lty=c(0,0,3,6),pch=c(1,16,46,46),legend=trend_labels,cex=0.8)
          else
            legend(legendxy[1],y=legendxy[2],lty=c(0,0,3,6),pch=c(1,16,46,46),legend=trend_labels,cex=0.8)
        }
      }
      if(TREND=="SM"){
        if(length(A)%%2==0){
          a1<-A[1:(length(A)/2)]
          a2<-A[(length(A)/2+1):length(A)]
          MTa1<-MTA[1:(length(A)/2)]
          MTa2<-MTA[(length(A)/2+1):length(A)]
          medVALUEa1<-median(a1,na.rm=TRUE)
          medVALUEa2<-median(a2,na.rm=TRUE)
          medTIMEa1<-median(MTa1,na.rm=TRUE)
          medTIMEa2<-median(MTa2,na.rm=TRUE)
          timeA<-c(medTIMEa1,medTIMEa2)
          valuesA<-c(medVALUEa1,medVALUEa2)
          interceptA<-coefficients(lm(valuesA~timeA))[1]
          slopeA<-coefficients(lm(valuesA~timeA))[2]
        }
        if(length(A)%%2==1){
          a11<-A[1:ceiling(length(A)/2)]
          a21<-A[(ceiling(length(A)/2)+1):length(A)]
          a12<-A[1:floor(length(A)/2)]
          a22<-A[(floor(length(A)/2)+1):length(A)]
          MTa11<-MTA[1:ceiling(length(A)/2)]
          MTa21<-MTA[(ceiling(length(A)/2)+1):length(A)]
          MTa12<-MTA[1:floor(length(A)/2)]
          MTa22<-MTA[(floor(length(A)/2)+1):length(A)]
          medVALUEa11<-median(a11,na.rm=TRUE)
          medVALUEa21<-median(a21,na.rm=TRUE)
          medVALUEa12<-median(a12,na.rm=TRUE)
          medVALUEa22<-median(a22,na.rm=TRUE)
          medTIMEa11<-median(MTa11,na.rm=TRUE)
          medTIMEa21<-median(MTa21,na.rm=TRUE)
          medTIMEa12<-median(MTa12,na.rm=TRUE)
          medTIMEa22<-median(MTa22,na.rm=TRUE)
          timeA1<-c(medTIMEa11,medTIMEa21)
          valuesA1<-c(medVALUEa11,medVALUEa21)
          interceptA1<-coefficients(lm(valuesA1~timeA1))[1]
          slopeA1<-coefficients(lm(valuesA1~timeA1))[2]
          timeA2<-c(medTIMEa12,medTIMEa22)
          valuesA2<-c(medVALUEa12,medVALUEa22)
          interceptA2<-coefficients(lm(valuesA2~timeA2))[1]
          slopeA2<-coefficients(lm(valuesA2~timeA2))[2]
        }
        if(length(B)%%2==0){
          b1<-B[1:(length(B)/2)]
          b2<-B[(length(B)/2+1):length(B)]
          MTb1<-MTB[1:(length(B)/2)]
          MTb2<-MTB[(length(B)/2+1):length(B)]
          medVALUEb1<-median(b1,na.rm=TRUE)
          medVALUEb2<-median(b2,na.rm=TRUE)
          medTIMEb1<-median(MTb1,na.rm=TRUE)
          medTIMEb2<-median(MTb2,na.rm=TRUE)
          timeB<-c(medTIMEb1,medTIMEb2)
          valuesB<-c(medVALUEb1,medVALUEb2)
          interceptB<-coefficients(lm(valuesB~timeB))[1]
          slopeB<-coefficients(lm(valuesB~timeB))[2]
        }
        if(length(B)%%2==1){
          b11<-B[1:ceiling(length(B)/2)]
          b21<-B[(ceiling(length(B)/2)+1):length(B)]
          b12<-B[1:floor(length(B)/2)]
          b22<-B[(floor(length(B)/2)+1):length(B)]
          MTb11<-MTB[1:ceiling(length(B)/2)]
          MTb21<-MTB[(ceiling(length(B)/2)+1):length(B)]
          MTb12<-MTB[1:floor(length(B)/2)]
          MTb22<-MTB[(floor(length(B)/2)+1):length(B)]
          medVALUEb11<-median(b11,na.rm=TRUE)
          medVALUEb21<-median(b21,na.rm=TRUE)
          medVALUEb12<-median(b12,na.rm=TRUE)
          medVALUEb22<-median(b22,na.rm=TRUE)
          medTIMEb11<-median(MTb11,na.rm=TRUE)
          medTIMEb21<-median(MTb21,na.rm=TRUE)
          medTIMEb12<-median(MTb12,na.rm=TRUE)
          medTIMEb22<-median(MTb22,na.rm=TRUE)
          timeB1<-c(medTIMEb11,medTIMEb21)
          valuesB1<-c(medVALUEb11,medVALUEb21)
          interceptB1<-coefficients(lm(valuesB1~timeB1))[1]
          slopeB1<-coefficients(lm(valuesB1~timeB1))[2]
          timeB2<-c(medTIMEb12,medTIMEb22)
          valuesB2<-c(medVALUEb12,medVALUEb22)
          interceptB2<-coefficients(lm(valuesB2~timeB2))[1]
          slopeB2<-coefficients(lm(valuesB2~timeB2))[2]
        }
        if(design=="CRD"|design=="ATD"|design=="RBD"|design=="Custom"){
          if(length(A)%%2==0){
            lines(c(1,MT),c(interceptA+slopeA,interceptA+slopeA*MT),lty=3)
          }
          if(length(A)%%2==1){
            lines(c(1,MT),c(interceptA1+slopeA1,interceptA1+slopeA1*MT),lty=3)
            lines(c(1,MT),c(interceptA2+slopeA2,interceptA2+slopeA2*MT),lty=3)
          }
          if(length(B)%%2==0){
            lines(c(1,MT),c(interceptB+slopeB,interceptB+slopeB*MT),lty=6)
          }
          if(length(B)%%2==1){
            lines(c(1,MT),c(interceptB1+slopeB1,interceptB1+slopeB1*MT),lty=6)
            lines(c(1,MT),c(interceptB2+slopeB2,interceptB2+slopeB2*MT),lty=6)
          }
          trend_labels<-c(labels[1:2],paste("trend line (split-middle)",labels[1:2]))
          if(is.null(legendxy))
            legend(locator(1),lty=c(0,0,3,6),pch=c(1,16,46,46),legend=trend_labels,cex=0.8)
          else
            legend(legendxy[1],y=legendxy[2],lty=c(0,0,3,6),pch=c(1,16,46,46),legend=trend_labels,cex=0.8)
        }
        if(design=="AB"){
          if(length(A)%%2==0){
            lines(c(1,length(A)),c(interceptA+slopeA,interceptA+slopeA*length(A)),lty=3)
          }
          if(length(A)%%2==1){
            lines(c(1,length(A)),c(interceptA1+slopeA1,interceptA1+slopeA1*length(A)),lty=3)
            lines(c(1,length(A)),c(interceptA2+slopeA2,interceptA2+slopeA2*length(A)),lty=3)
          }
          if(length(B)%%2==0){
            lines(c(length(A)+1,MT),c(interceptB+slopeB*(length(A)+1),interceptB+slopeB*MT),lty=3)
          }
          if(length(B)%%2==1){
            lines(c(length(A)+1,MT),c(interceptB1+slopeB1,interceptB1+slopeB1*MT),lty=3)
            lines(c(length(A)+1,MT),c(interceptB2+slopeB2*(length(A)+1),interceptB2+slopeB2*MT),lty=3)
          }
        }
      }
      if(TREND=="RTL"){
        if(length(A)%%3==0){
          a1<-A[1:(length(A)/3)]
          a2<-A[(length(A)/3+1):(length(A)/3*2)]
          a3<-A[(length(A)/3*2+1):(length(A))]
          MTa1<-MTA[1:(length(A)/3)]
          MTa2<-MTA[(length(A)/3+1):(length(A)/3*2)]
          MTa3<-MTA[(length(A)/3*2+1):(length(A))]
        }
        if(length(A)%%3==1){
          a1<-A[1:floor(length(A)/3)]
          a2<-A[(floor(length(A)/3)+1):(floor(length(A)/3)*2+1)]
          a3<-A[(floor(length(A)/3)*2+2):length(A)]
          MTa1<-MTA[1:floor(length(A)/3)]
          MTa2<-MTA[(floor(length(A)/3)+1):(floor(length(A)/3)*2+1)]
          MTa3<-MTA[(floor(length(A)/3)*2+2):length(A)]
        }
        if(length(A)%%3==2){
          a1<-A[1:(ceiling(length(A)/3))]
          a2<-A[(ceiling(length(A)/3)+1):(ceiling(length(A)/3)+floor(length(A)/3))]
          a3<-A[(ceiling(length(A)/3)+floor(length(A)/3)+1):length(A)]
          MTa1<-MTA[1:(ceiling(length(A)/3))]
          MTa2<-MTA[(ceiling(length(A)/3)+1):(ceiling(length(A)/3)+floor(length(A)/3))]
          MTa3<-MTA[(ceiling(length(A)/3)+floor(length(A)/3)+1):length(A)]
        }
        if(length(B)%%3==0){
          b1<-B[1:(length(B)/3)]
          b2<-B[(length(B)/3+1):(length(B)/3*2)]
          b3<-B[(length(B)/3*2+1):(length(B))]
          MTb1<-MTB[1:(length(B)/3)]
          MTb2<-MTB[(length(B)/3+1):(length(B)/3*2)]
          MTb3<-MTB[(length(B)/3*2+1):(length(B))]
        }
        if(length(B)%%3==1){
          b1<-B[1:floor(length(B)/3)]
          b2<-B[(floor(length(B)/3)+1):(floor(length(B)/3)*2+1)]
          b3<-B[(floor(length(B)/3)*2+2):length(B)]
          MTb1<-MTB[1:floor(length(B)/3)]
          MTb2<-MTB[(floor(length(B)/3)+1):(floor(length(B)/3)*2+1)]
          MTb3<-MTB[(floor(length(B)/3)*2+2):length(B)]
        }
        if(length(B)%%3==2){
          b1<-B[1:(ceiling(length(B)/3))]
          b2<-B[(ceiling(length(B)/3)+1):(ceiling(length(B)/3)+floor(length(B)/3))]
          b3<-B[(ceiling(length(B)/3)+floor(length(B)/3)+1):length(B)]
          MTb1<-MTB[1:(ceiling(length(B)/3))]
          MTb2<-MTB[(ceiling(length(B)/3)+1):(ceiling(length(B)/3)+floor(length(B)/3))]
          MTb3<-MTB[(ceiling(length(B)/3)+floor(length(B)/3)+1):length(B)]
        }
        medTIMEa1<-median(MTa1,na.rm=TRUE)
        medTIMEa2<-median(MTa2,na.rm=TRUE)
        medTIMEa3<-median(MTa3,na.rm=TRUE)
        medVALUEa1<-median(a1,na.rm=TRUE)
        medVALUEa2<-median(a2,na.rm=TRUE)
        medVALUEa3<-median(a3,na.rm=TRUE)
        medTIMEb1<-median(MTb1,na.rm=TRUE)
        medTIMEb2<-median(MTb2,na.rm=TRUE)
        medTIMEb3<-median(MTb3,na.rm=TRUE)
        medVALUEb1<-median(b1,na.rm=TRUE)
        medVALUEb2<-median(b2,na.rm=TRUE)
        medVALUEb3<-median(b3,na.rm=TRUE)
        slopeA<-(medVALUEa3-medVALUEa1)/(medTIMEa3-medTIMEa1)
        interceptA<-(1/3)*((medVALUEa1+medVALUEa2+medVALUEa3)-slopeA*(medTIMEa1+medTIMEa2+medTIMEa3))
        slopeB<-(medVALUEb3-medVALUEb1)/(medTIMEb3-medTIMEb1)
        interceptB<-(1/3)*((medVALUEb1+medVALUEb2+medVALUEb3)-slopeB*(medTIMEb1+medTIMEb2+medTIMEb3))
        if(design=="CRD"|design=="ATD"|design=="RBD"|design=="Custom"){
          lines(c(1,MT),c(interceptA+slopeA,interceptA+slopeA*MT),lty=3)
          lines(c(1,MT),c(interceptB+slopeB,interceptB+slopeB*MT),lty=6)
        }
        if(design=="AB"){
          lines(c(1,length(A)),c(interceptA+slopeA,interceptA+slopeA*length(A)),lty=3)
          lines(c(length(A)+1,MT),c(interceptB+slopeB*(length(A)+1),interceptB+slopeB*MT),lty=6)
        }
        points(medTIMEa1,medVALUEa1,cex=1.2)
        points(medTIMEa2,medVALUEa2,cex=1.2)
        points(medTIMEa3,medVALUEa3,cex=1.2)
        lines(c(medTIMEa1,medTIMEa2),c(medVALUEa1,medVALUEa2),lty=3)
        lines(c(medTIMEa2,medTIMEa3),c(medVALUEa2,medVALUEa3),lty=3)
        points(medTIMEb1,medVALUEb1,cex=1.2)
        points(medTIMEb2,medVALUEb2,cex=1.2)
        points(medTIMEb3,medVALUEb3,cex=1.2)
        lines(c(medTIMEb1,medTIMEb2),c(medVALUEb1,medVALUEb2),lty=6)
        lines(c(medTIMEb2,medTIMEb3),c(medVALUEb2,medVALUEb3),lty=6)
        if(design=="CRD"|design=="ATD"|design=="RBD"|design=="Custom"){
          trend_labels<-c(labels[1:2],paste("resistant trend line",labels[1:2]))
          if(is.null(legendxy))
            legend(locator(1),lty=c(0,0,3,6),pch=c(1,16,46,46),legend=trend_labels,cex=0.8)
          else
            legend(legendxy[1],y=legendxy[2],lty=c(0,0,3,6),pch=c(1,16,46,46),legend=trend_labels,cex=0.8)
        }
      }
      if(TREND=="RM3"){
        RM3a<-numeric()
        for(it in 1:(length(A)-2)){
          RM3a<-c(RM3a,median(A[it:(it+2)],na.rm=TRUE))
        }
        times3A<-numeric()
        for(it in 1:(length(A)-2)){
          times3A<-c(times3A,median(MTA[it:(it+2)],na.rm=TRUE))
        }
        for(it in 1:(length(A)-2)){
          points(times3A[it],RM3a[it],pch=3)
          lines(c(times3A[it],times3A[it+1]),c(RM3a[it],RM3a[it+1]))
        }
        RM3b<-numeric()
        for(it in 1:(length(B)-2)){
          RM3b<-c(RM3b,median(B[it:(it+2)],na.rm=TRUE))
        }
        times3B<-numeric()
        for(it in (1:(length(B)-2))){
          times3B<-c(times3B,median(MTB[it:(it+2)],na.rm=TRUE))
        }
        for(it in 1:(length(B)-2)){
          points(times3B[it],RM3b[it],pch=4)
          lines(c(times3B[it],times3B[it+1]),c(RM3b[it],RM3b[it+1]),lty=2)
        }
        if(design=="CRD"|design=="ATD"|design=="RBD"|design=="Custom"){
          trend_labels<-c(labels[1:2],paste("running medians (3)",labels[1:2]))
          if(is.null(legendxy))
            legend(locator(1),lty=c(0,0,1,2),pch=c(1,16,3,4),legend=trend_labels,cex=0.8)
          else
            legend(legendxy[1],y=legendxy[2],lty=c(0,0,1,2),pch=c(1,16,3,4),legend=trend_labels,cex=0.8)
        }
      }
      if(TREND=="RM5"){
        RM5a<-numeric()
        for(it in 1:(length(A)-4)){
          RM5a<-c(RM5a,median(A[it:(it+4)],na.rm=TRUE))
        }
        times5A<-numeric()
        for(it in 1:(length(A)-4)){
          times5A<-c(times5A,median(MTA[it:(it+4)],na.rm=TRUE))
        }
        for(it in 1:(length(A)-4)){
          points(times5A[it],RM5a[it],pch=3)
          lines(c(times5A[it],times5A[it+1]),c(RM5a[it],RM5a[it+1]))
        }
        RM5b<-numeric()
        for(it in 1:(length(B)-4)){
          RM5b<-c(RM5b,median(B[it:(it+4)],na.rm=TRUE))
        }
        times5B<-numeric()
        for(it in (1:(length(B)-4))){
          times5B<-c(times5B,median(MTB[it:(it+4)],na.rm=TRUE))
        }
        for(it in 1:(length(B)-4)){
          points(times5B[it],RM5b[it],pch=4)
          lines(c(times5B[it],times5B[it+1]),c(RM5b[it],RM5b[it+1]),lty=2)
        }
        if(design=="CRD"|design=="ATD"|design=="RBD"|design=="Custom"){
          trend_labels<-c(labels[1:2],paste("running medians (5)",labels[1:2]))
          if(is.null(legendxy))
            legend(locator(1),lty=c(0,0,1,2),pch=c(1,16,3,4),legend=trend_labels,cex=0.8)
          else
            legend(legendxy[1],y=legendxy[2],lty=c(0,0,1,2),pch=c(1,16,3,4),legend=trend_labels,cex=0.8)
        }
      }
      if(TREND=="RM42"){
        RM4a<-numeric()
        for(it in 1:(length(A)-3)){
          RM4a<-c(RM4a,median(A[it:(it+3)],na.rm=TRUE))
        }
        RM42a<-numeric()
        for(it in 1:(length(RM4a)-1)){
          RM42a<-c(RM42a,mean(RM4a[it:(it+1)],na.rm=TRUE))
        }
        times42A<-numeric()
        for(it in 1:(length(A)-4)){
          times42A<-c(times42A,median(MTA[it:(it+4)],na.rm=TRUE))
        }
        for(it in 1:(length(A)-4)){
          points(times42A[it],RM42a[it],pch=3)
          lines(c(times42A[it],times42A[it+1]),c(RM42a[it],RM42a[it+1]))
        }
        RM4b<-numeric()
        for(it in 1:(length(B)-3)){
          RM4b<-c(RM4b,median(B[it:(it+3)],na.rm=TRUE))
        }
        RM42b<-numeric()
        for(it in 1:(length(RM4b)-1)){
          RM42b<-c(RM42b,mean(RM4b[it:(it+1)],na.rm=TRUE))
        }
        times42B<-numeric()
        for(it in 1:(length(B)-4)){
          times42B<-c(times42B,median(MTB[it:(it+4)],na.rm=TRUE))
        }
        for(it in 1:(length(B)-4)){
          points(times42B[it],RM42b[it],pch=4)
          lines(c(times42B[it],times42B[it+1]),c(RM42b[it],RM42b[it+1]),lty=2)
        }
        if(design=="CRD"|design=="ATD"|design=="RBD"|design=="Custom"){
          trend_labels<-c(labels[1:2],paste("running medians (42)",labels[1:2]))
          if(is.null(legendxy))
            legend(locator(1),lty=c(0,0,1,2),pch=c(1,16,3,4),legend=trend_labels,cex=0.8)
          else
            legend(legendxy[1],y=legendxy[2],lty=c(0,0,1,2),pch=c(1,16,3,4),legend=trend_labels,cex=0.8)
        }
      }
    }
  }
  
  if(design=="ABA"|design=="ABAB"){
    A1<-data[,2][data[,1]=="A1"]
    B1<-data[,2][data[,1]=="B1"]
    A2<-data[,2][data[,1]=="A2"]
    B2<-data[,2][data[,1]=="B2"]
    if(TREND=="VLP"){
      if(CL=="mean"){
        CLA1<-mean(A1,na.rm=TRUE)
        CLB1<-mean(B1,na.rm=TRUE)
        CLA2<-mean(A2,na.rm=TRUE)
        CLB2<-mean(B2,na.rm=TRUE)
      }
      if(CL=="median"){
        CLA1<-median(A1,na.rm=TRUE)
        CLB1<-median(B1,na.rm=TRUE)
        CLA2<-median(A2,na.rm=TRUE)
        CLB2<-median(B2,na.rm=TRUE)
      }
      if(CL=="bmed"){
        aa1<-sort(A1)
        bb1<-sort(B1)
        aa2<-sort(A2)
        bb2<-sort(B2)
        if(length(aa1)<5){
          CLA1<-median(data[,2][data[,1]=="A1"],na.rm=TRUE)
        }
        if(length(aa1)==5|length(aa1)==7|length(aa1)==9|length(aa1)==11){
          CLA1<-(aa1[ceiling(length(aa1)/2)-1]+aa1[ceiling(length(aa1)/2)]+aa1[ceiling(length(aa1)/2)+1])/3
        }
        if(length(aa1)>=13&length(aa1)%%2==1){
          CLA1<-(aa1[ceiling(length(aa1)/2)-2]+aa1[ceiling(length(aa1)/2)-1]+aa1[ceiling(length(aa1)/2)]+aa1[ceiling(length(aa1)/2)+1] +aa1[ceiling(length(aa1)/2)+2])/5
        }
        if(length(aa1)==6|length(aa1)==8|length(aa1)==10|length(aa1)==12){
          CLA1<-1/6*aa1[length(aa1)/2-1]+1/3*aa1[length(aa1)/2]+1/3*aa1[length(aa1)/2+1]+1/6*aa1[length(aa1)/2+2]
        }
        if(length(aa1)>13&length(aa1)%%2==0){
          CLA1<-1/10*aa1[length(aa1)/2-2]+1/5*aa1[length(aa1)/2-1]+1/5*aa1[length(aa1)/2]+1/5*aa1[length(aa1)/2+1]+1/5*aa1[length(aa1)/2+2]+1/10*aa1[length(aa1)/2+3]
        }
        if(length(bb1)<5){
          CLB1<-median(data[,2][data[,1]=="B1"],na.rm=TRUE)
        }
        if(length(bb1)==5|length(bb1)==7|length(bb1)==9|length(bb1)==11){
          CLB1<-(bb1[ceiling(length(bb1)/2)-1]+bb1[ceiling(length(bb1)/2)]+bb1[ceiling(length(bb1)/2)+1])/3
        }
        if(length(bb1)>=13&length(bb1)%%2==1){
          CLB1<-(bb1[ceiling(length(bb1)/2)-2]+bb1[ceiling(length(bb1)/2)-1]+bb1[ceiling(length(bb1)/2)]+bb1[ceiling(length(bb1)/2)+1]+bb1[ceiling(length(bb1)/2)+2])/5
        }
        if(length(bb1)==6|length(bb1)==8|length(bb1)==10|length(bb1)==12){
          CLB1<-1/6*bb1[length(bb1)/2-1]+1/3*bb1[length(bb1)/2]+1/3*bb1[length(bb1)/2+1]+1/6*bb1[length(bb1)/2+2]
        }
        if(length(bb1)>13&length(bb1)%%2==0){
          CLB1<-1/10*bb1[length(bb1)/2-2]+1/5*bb1[length(bb1)/2-1]+1/5*bb1[length(bb1)/2]+1/5*bb1[length(bb1)/2+1]+1/5*bb1[length(bb1)/2+2]+1/10*bb1[length(bb1)/2+3]
        }
        if(length(aa2)<5){
          CLA2<-median(data[,2][data[,1]=="A2"],na.rm=TRUE)
        }
        if(length(aa2)==5|length(aa2)==7|length(aa2)==9|length(aa2)==11){
          CLA2<-(aa2[ceiling(length(aa2)/2)-1]+aa2[ceiling(length(aa2)/2)]+aa2[ceiling(length(aa2)/2)+1])/3
        }
        if(length(aa2)>=13&length(aa2)%%2==1){
          CLA2<-(aa2[ceiling(length(aa2)/2)-2]+aa2[ceiling(length(aa2)/2)-1]+aa2[ceiling(length(aa2)/2)]+aa2[ceiling(length(aa2)/2)+1] +aa2[ceiling(length(aa2)/2)+2])/5
        }
        if(length(aa2)==6|length(aa2)==8|length(aa2)==10|length(aa2)==12){
          CLA2<-1/6*aa2[length(aa2)/2-1]+1/3*aa2[length(aa2)/2]+1/3*aa2[length(aa2)/2+1]+1/6*aa2[length(aa2)/2+2]
        }
        if(length(aa2)>13&length(aa2)%%2==0){
          CLA2<-1/10*aa2[length(aa2)/2-2]+1/5*aa2[length(aa2)/2-1]+1/5*aa2[length(aa2)/2]+1/5*aa2[length(aa2)/2+1]+1/5*aa2[length(aa2)/2+2]+1/10*aa2[length(aa2)/2+3]
        }
        if(length(bb2)<5){
          CLB2<-median(data[,2][data[,1]=="B2"])
        }
        if(length(bb2)==5|length(bb2)==7|length(bb2)==9|length(bb2)==11){
          CLB2<-(bb2[ceiling(length(bb2)/2)-1]+bb2[ceiling(length(bb2)/2)]+bb2[ceiling(length(bb2)/2)+1])/3
        }
        if(length(bb2)>=13&length(bb2)%%2==1){
          CLB2<-(bb2[ceiling(length(bb2)/2)-2]+bb2[ceiling(length(bb2)/2)-1]+bb2[ceiling(length(bb2)/2)]+bb2[ceiling(length(bb2)/2)+1]+bb2[ceiling(length(bb2)/2)+2])/5
        }
        if(length(bb2)==6|length(bb2)==8|length(bb2)==10|length(bb2)==12){
          CLB2<-1/6*bb2[length(bb2)/2-1]+1/3*bb2[length(bb2)/2]+1/3*bb2[length(bb2)/2+1]+1/6*bb2[length(bb2)/2+2]
        }
        if(length(bb2)>13&length(bb2)%%2==0){
          CLB2<-1/10*bb2[length(bb2)/2-2]+1/5*bb2[length(bb2)/2-1]+1/5*bb2[length(bb2)/2]+1/5*bb2[length(bb2)/2+1]+1/5*bb2[length(bb2)/2+2]+1/10*bb2[length(bb2)/2+3]
        }
      }
      if(CL=="trimmean"){
        CLA1<-mean(A1,trim=tr,na.rm=TRUE)
        CLB1<-mean(B1,trim=tr,na.rm=TRUE)
        CLA2<-mean(A2,trim=tr,na.rm=TRUE)
        CLB2<-mean(B2,trim=tr,na.rm=TRUE)
      }
      if(CL=="mest"){
        hpsi<-function(x,bend=1.28){
          hpsi<-ifelse(abs(x)<=bend,x,bend*sign(x))
          hpsi
        }
        mest<-function(x,bend=1.28,na.rm=TRUE){
          if(na.rm)x<-x[!is.na(x)]
          if(length(x)==0) return(NA)
          if(mad(x)==0)stop("MAD=0. The M-estimator cannot be computed.")
          y<-(x-median(x))/mad(x)
          A<-sum(hpsi(y,bend))
          B<-length(x[abs(y)<=bend])
          mest<-median(x)+mad(x)*A/B
          repeat{
            y<-(x-mest)/mad(x)
            A<-sum(hpsi(y,bend))
            B<-length(x[abs(y)<=bend])
            newmest<-mest+mad(x)*A/B
            if(abs(newmest-mest) <.0001)break
            mest<-newmest
          }
          mest
        }
        CLA1<-mest(A1,bend=tr)
        CLB1<-mest(B1,bend=tr)
        CLA2<-mest(A2,bend=tr)
        CLB2<-mest(B2,bend=tr)
      }
      residualsA1<-A1-CLA1
      residualsB1<-B1-CLB1
      residualsA2<-A2-CLA2
      residualsB2<-B2-CLB2
      residuals<-c(residualsA1,residualsB1,residualsA2,residualsB2)
      if (is.null(xlab)){
        xlab <- "Measurement Times"
      }
      if (is.null(ylab)){
        ylab <- "Residuals"
      }
      plot(x,residuals,xlab=xlab,ylab=ylab,ylim=ylim,type="n")
      lines(c(1,length(residualsA1)),c(0,0))
      for(it in 1:length(residualsA1)){
        lines(c(it,it),c(residualsA1[it],0))
      }
      lines(c(length(residualsA1)+1,length(residualsA1)+length(residualsB1)),c(0,0))
      for(it in 1:length(residualsB1)){
        lines(c(length(residualsA1)+it,length(residualsA1)+it),c(residualsB1[it],0))
      }
      lines(c(length(residualsA1)+length(residualsB1)+1,length(residualsA1)+length(residualsB1)+length(residualsA2)),c(0,0))
      for(it in 1:length(residualsA2)){
        lines(c(length(residualsA1)+length(residualsB1)+it,length(residualsA1)+length(residualsB1)+it),c(residualsA2[it],0))
      }
      mtext(labels[1],side=3,at=(sum(data[,1]=="A1")+1)/2)
      mtext(labels[2],side=3,at=(sum(data[,1]=="A1")+(sum(data[,1]=="B1")+1)/2))
      mtext(labels[3],side=3,at=(sum(data[,1]=="A1")+sum(data[,1]=="B1")+(sum(data[,1]=="A2")+1)/2))
      if(design=="ABAB"){
        lines(c(length(residualsA1)+length(residualsB1)+length(residualsA2)+1,MT),c(0,0))
        for(it in 1:length(residualsB2)){
          lines(c(length(residualsA1)+length(residualsB1)+length(residualsA2)+it,length(residualsA1)+length(residualsB1)+length(residualsA2)+it),c(residualsB2[it],0))
        }
        mtext(labels[4],side=3,at=(sum(data[,1]=="A1")+sum(data[,1]=="B1")+sum(data[,1]=="A2")+(sum(data[,1]=="B2")+1)/2))
      }
    }
    if(TREND=="RTL"|TREND=="SM"|TREND=="LSR"|TREND=="RM3"|TREND=="RM5"|TREND=="RM42"){
      if (is.null(xlab)){
        xlab <- "Measurement Times"
      }
      if (is.null(ylab)){
        ylab <- "Scores"
      }
      plot(x,data[,2],xlab=xlab,ylab=ylab,ylim=ylim,pch=16)
      lines(c(sum(data[,1]=="A1")+0.5,sum(data[,1]=="A1")+0.5),c(min(data[,2],ylim[1],na.rm=TRUE)-5,max(data[,2],ylim[2],na.rm=TRUE)+5),lty=2)
      lines(c(sum(data[,1]=="A1")+sum(data[,1]=="B1")+0.5,sum(data[,1]=="A1")+sum(data[,1]=="B1")+0.5),
            c(min(data[,2],ylim[1],na.rm=TRUE)-5,max(data[,2],ylim[2],na.rm=TRUE)+5),lty=2)
      mtext(labels[1],side=3,at=(sum(data[,1]=="A1")+1)/2)
      mtext(labels[2],side=3,at=(sum(data[,1]=="A1")+(sum(data[,1]=="B1")+1)/2))
      mtext(labels[3],side=3,at=(sum(data[,1]=="A1")+sum(data[,1]=="B1")+(sum(data[,1]=="A2")+1)/2))
      if(design=="ABAB"){
        lines(c(sum(data[,1]=="A1")+sum(data[,1]=="B1")+sum(data[,1]=="A2")+0.5,sum(data[,1]=="A1")+sum(data[,1]=="B1")+sum(data[,1]=="A2")+0.5),
              c(min(data[,2],ylim[1],na.rm=TRUE)-5,max(data[,2],ylim[2],na.rm=TRUE)+5),lty=2)
        mtext(labels[4],side=3,at=(sum(data[,1]=="A1")+sum(data[,1]=="B1")+sum(data[,1]=="A2")+(sum(data[,1]=="B2")+1)/2))
      }
      if(TREND=="LSR"){
        xa1<-1:length(A1)
        interceptLSRa1<-coefficients(lm(A1~xa1,na.action=na.omit))[1]
        slopeLSRa1<-coefficients(lm(A1~xa1,na.action=na.omit))[2]
        xb1<-1:length(B1)
        interceptLSRb1<-coefficients(lm(B1~xb1,na.action=na.omit))[1]
        slopeLSRb1<-coefficients(lm(B1~xb1,na.action=na.omit))[2]	
        xa2<-1:length(A2)
        interceptLSRa2<-coefficients(lm(A2~xa2,na.action=na.omit))[1]
        slopeLSRa2<-coefficients(lm(A2~xa2,na.action=na.omit))[2]
        lines(c(1,length(A1)),c(interceptLSRa1+slopeLSRa1,interceptLSRa1+slopeLSRa1*length(A1)),lty=2)
        lines(c(length(A1)+1,(length(A1)+length(B1))),c(interceptLSRb1+slopeLSRb1,interceptLSRb1+slopeLSRb1*length(B1)),lty=2)
        lines(c(length(A1)+length(B1)+1,(length(A1)+length(B1)+length(A2))),c(interceptLSRa2+slopeLSRa2,interceptLSRa2+slopeLSRa2*length(A2)),lty=2)
        if(design=="ABAB"){
          xb2<-1:length(B2)
          interceptLSRb2<-coefficients(lm(B2~xb2,na.action=na.omit))[1]
          slopeLSRb2<-coefficients(lm(B2~xb2,na.action=na.omit))[2]
          lines(c(length(A1)+length(B1)+length(A2)+1,MT),c(interceptLSRb2+slopeLSRb2,interceptLSRb2+slopeLSRb2*length(B2)),lty=2)
        }
      }
      if(TREND=="SM"){
        if(length(A1)%%2==0){
          a11<-A1[1:(length(A1)/2)]
          a12<-A1[(length(A1)/2+1):length(A1)]
          medTIMEa11<-median(x[1:length(a11)])
          medTIMEa12<-median(x[(length(a11)+1):length(A1)])
          medVALUEa11<-median(a11,na.rm=TRUE)
          medVALUEa12<-median(a12,na.rm=TRUE)
          timeA1<-c(medTIMEa11,medTIMEa12)
          valuesA1<-c(medVALUEa11,medVALUEa12)
          interceptA1<-coefficients(lm(valuesA1~timeA1,na.action=na.omit))[1]
          slopeA1<-coefficients(lm(valuesA1~timeA1,na.action=na.omit))[2]
          lines(c(1,length(A1)),c(interceptA1+slopeA1,interceptA1+slopeA1*length(A1)),lty=3)
        }
        if(length(A1)%%2==1){
          a111<-A1[1:ceiling(length(A1)/2)]
          a121<-A1[(ceiling(length(A1)/2)+1):length(A1)]
          a112<-A1[1:floor(length(A1)/2)]
          a122<-A1[(floor(length(A1)/2)+1):length(A1)]
          medTIMEa111<-median(x[1:length(a111)])
          medTIMEa121<-median(x[(length(a111)+1):length(A1)])
          medTIMEa112<-median(x[1:length(a112)])
          medTIMEa122<-median(x[(length(a112)+1):length(A1)])
          medVALUEa111<-median(a111,na.rm=TRUE)
          medVALUEa121<-median(a121,na.rm=TRUE)
          medVALUEa112<-median(a112,na.rm=TRUE)
          medVALUEa122<-median(a122,na.rm=TRUE)
          timeA11<-c(medTIMEa111,medTIMEa121)
          valuesA11<-c(medVALUEa111,medVALUEa121)
          interceptA11<-coefficients(lm(valuesA11~timeA11))[1]
          slopeA11<-coefficients(lm(valuesA11~timeA11))[2]
          lines(c(1,length(A1)),c(interceptA11+slopeA11,interceptA11+slopeA11*length(A1)),lty=3)
          timeA12<-c(medTIMEa112,medTIMEa122)
          valuesA12<-c(medVALUEa112,medVALUEa122)
          interceptA12<-coefficients(lm(valuesA12~timeA12))[1]
          slopeA12<-coefficients(lm(valuesA12~timeA12))[2]
          lines(c(1,length(A1)),c(interceptA12+slopeA12,interceptA12+slopeA12*length(A1)),lty=3)
        }
        if(length(B1)%%2==0){
          b11<-B1[1:(length(B1)/2)]
          b12<-B1[(length(B1)/2+1):length(B1)]
          medTIMEb11<-median(x[(length(A1)+1):(length(A1)+length(b11))])
          medTIMEb12<-median(x[(length(A1)+length(b11)+1):(length(A1)+length(B1))])
          medVALUEb11<-median(b11,na.rm=TRUE)
          medVALUEb12<-median(b12,na.rm=TRUE)
          timeB1<-c(medTIMEb11,medTIMEb12)
          valuesB1<-c(medVALUEb11,medVALUEb12)
          interceptB1<-coefficients(lm(valuesB1~timeB1))[1]
          slopeB1<-coefficients(lm(valuesB1~timeB1))[2]
          lines(c(length(A1)+1,length(A1)+length(B1)),c(interceptB1+slopeB1*(length(A1)+1),interceptB1+slopeB1*(length(A1)+length(B1))),lty=3)
        }
        if(length(B1)%%2==1){
          b111<-B1[1:ceiling(length(B1)/2)]
          b121<-B1[(ceiling(length(B1)/2)+1):length(B1)]
          b112<-B1[1:floor(length(B1)/2)]
          b122<-B1[(floor(length(B1)/2)+1):length(B1)]
          medTIMEb111<-median(x[(length(A1)+1):(length(A1)+length(b111))])
          medTIMEb121<-median(x[(length(A1)+length(b111)+1):(length(A1)+length(B1))])
          medTIMEb112<-median(x[(length(A1)+1):(length(A1)+length(b112))])
          medTIMEb122<-median(x[(length(A1)+length(b112)+1):(length(A1)+length(B1))])
          medVALUEb111<-median(b111,na.rm=TRUE)
          medVALUEb121<-median(b121,na.rm=TRUE)
          medVALUEb112<-median(b112,na.rm=TRUE)
          medVALUEb122<-median(b122,na.rm=TRUE)
          timeB11<-c(medTIMEb111,medTIMEb121)
          valuesB11<-c(medVALUEb111,medVALUEb121)
          interceptB11<-coefficients(lm(valuesB11~timeB11))[1]
          slopeB11<-coefficients(lm(valuesB11~timeB11))[2]
          lines(c(length(A1)+1,length(A1)+length(B1)),c(interceptB11+slopeB11*(length(A1)+1),interceptB11+slopeB11*(length(A1)+length(B1))),lty=3)
          timeB12<-c(medTIMEb112,medTIMEb122)
          valuesB12<-c(medVALUEb112,medVALUEb122)
          interceptB12<-coefficients(lm(valuesB12~timeB12))[1]
          slopeB12<-coefficients(lm(valuesB12~timeB12))[2]
          lines(c(length(A1)+1,length(A1)+length(B1)),c(interceptB12+slopeB12*(length(A1)+1),interceptB12+slopeB12*(length(A1)+length(B1))),lty=3)
        }
        if(length(A2)%%2==0){
          a21<-A2[1:(length(A2)/2)]
          a22<-A2[(length(A2)/2+1):length(A2)]
          medTIMEa21<-median(x[(length(A1)+length(B1)+1):(length(A1)+length(B1)+length(a21))])
          medTIMEa22<-median(x[(length(A1)+length(B1)+length(a21)+1):(length(A1)+length(B1)+length(A2))])
          medVALUEa21<-median(a21,na.rm=TRUE)
          medVALUEa22<-median(a22,na.rm=TRUE)
          timeA2<-c(medTIMEa21,medTIMEa22)
          valuesA2<-c(medVALUEa21,medVALUEa22)
          interceptA2<-coefficients(lm(valuesA2~timeA2))[1]
          slopeA2<-coefficients(lm(valuesA2~timeA2))[2]
          lines(c(length(A1)+length(B1)+1,length(A1)+length(B1)+length(A2)),c(interceptA2+slopeA2*(length(A1)+length(B1)+1),interceptA2+slopeA2*(length(A1)+length(B1)+length(A2))),lty=3)
        }
        if(length(A2)%%2==1){
          a211<-A2[1:ceiling(length(A2)/2)]
          a221<-A2[(ceiling(length(A2)/2)+1):length(A2)]
          a212<-A2[1:floor(length(A2)/2)]
          a222<-A2[(floor(length(A2)/2)+1):length(A2)]
          medTIMEa211<-median(x[(length(A1)+length(B1)+1):(length(A1)+length(B1)+length(a211))])
          medTIMEa221<-median(x[(length(A1)+length(B1)+length(a211)+1):(length(A1)+length(B1)+length(A2))])
          medTIMEa212<-median(x[(length(A1)+length(B1)+1):(length(A1)+length(B1)+length(a212))])
          medTIMEa222<-median(x[(length(A1)+length(B1)+length(a212)+1):(length(A1)+length(B1)+length(A2))])
          medVALUEa211<-median(a211,na.rm=TRUE)
          medVALUEa221<-median(a221,na.rm=TRUE)
          medVALUEa212<-median(a212,na.rm=TRUE)
          medVALUEa222<-median(a222,na.rm=TRUE)
          timeA21<-c(medTIMEa211,medTIMEa221)
          valuesA21<-c(medVALUEa211,medVALUEa221)
          interceptA21<-coefficients(lm(valuesA21~timeA21))[1]
          slopeA21<-coefficients(lm(valuesA21~timeA21))[2]
          lines(c(length(A1)+length(B1)+1,length(A1)+length(B1)+length(A2)),c(interceptA21+slopeA21*(length(A1)+length(B1)+1),interceptA21+slopeA21*(length(A1)+length(B1)+length(A2))),lty=3)
          timeA22<-c(medTIMEa212,medTIMEa222)
          valuesA22<-c(medVALUEa212,medVALUEa222)
          interceptA22<-coefficients(lm(valuesA22~timeA22))[1]
          slopeA22<-coefficients(lm(valuesA22~timeA22))[2]
          lines(c(length(A1)+length(B1)+1,length(A1)+length(B1)+length(A2)),c(interceptA22+slopeA22*(length(A1)+length(B1)+1),interceptA22+slopeA22*(length(A1)+length(B1)+length(A2))),lty=3)
        }
        if(design=="ABAB"){
          if(length(B2)%%2==0){
            b21<-B2[1:(length(B2)/2)]
            b22<-B2[(length(B2)/2+1):length(B2)]
            medTIMEb21<-median(x[(length(A1)+length(B1)+length(A2)+1):(length(A1)+length(B1)+length(A2)+length(b21))])
            medTIMEb22<-median(x[(length(A1)+length(B1)+length(A2)+length(b21)+1):MT])
            medVALUEb21<-median(b21,na.rm=TRUE)
            medVALUEb22<-median(b22,na.rm=TRUE)
            timeB2<-c(medTIMEb21,medTIMEb22)
            valuesB2<-c(medVALUEb21,medVALUEb22)
            interceptB2<-coefficients(lm(valuesB2~timeB2))[1]
            slopeB2<-coefficients(lm(valuesB2~timeB2))[2]
            lines(c(length(A1)+length(B1)+length(A2)+1,MT),c(interceptB2+slopeB2*(length(A1)+length(B1)+length(A2)+1),interceptB2+slopeB2*MT),lty=3)
          }
          if(length(B2)%%2==1){
            b211<-B2[1:ceiling(length(B2)/2)]
            b221<-B2[(ceiling(length(B2)/2)+1):length(B2)]
            b212<-B2[1:floor(length(B2)/2)]
            b222<-B2[(floor(length(B2)/2)+1):length(B2)]
            medTIMEb211<-median(x[(length(A1)+length(B1)+length(A2)+1):(length(A1)+length(B1)+length(A2)+length(b211))])
            medTIMEb221<-median(x[(length(A1)+length(B1)+length(A2)+length(b211)+1):MT])
            medTIMEb212<-median(x[(length(A1)+length(B1)+length(A2)+1):(length(A1)+length(B1)+length(A2)+length(b212))])
            medTIMEb222<-median(x[(length(A1)+length(B1)+length(A2)+length(b212)+1):MT])
            medVALUEb211<-median(b211,na.rm=TRUE)
            medVALUEb221<-median(b221,na.rm=TRUE)
            medVALUEb212<-median(b212,na.rm=TRUE)
            medVALUEb222<-median(b222,na.rm=TRUE)
            timeB21<-c(medTIMEb211,medTIMEb221)
            valuesB21<-c(medVALUEb211,medVALUEb221)
            interceptB21<-coefficients(lm(valuesB21~timeB21))[1]
            slopeB21<-coefficients(lm(valuesB21~timeB21))[2]
            lines(c(length(A1)+length(B1)+length(A2)+1,MT),c(interceptB21+slopeB21*(length(A1)+length(B1)+length(A2)+1),interceptB21+slopeB21*MT),lty=3)
            timeB22<-c(medTIMEb212,medTIMEb222)
            valuesB22<-c(medVALUEb212,medVALUEb222)
            interceptB22<-coefficients(lm(valuesB22~timeB22))[1]
            slopeB22<-coefficients(lm(valuesB22~timeB22))[2]
            lines(c(length(A1)+length(B1)+length(A2)+1,MT),c(interceptB22+slopeB22*(length(A1)+length(B1)+length(A2)+1),interceptB22+slopeB22*MT),lty=3)
          }
        }
      } 
      if(TREND=="RTL"){
        if(length(A1)%%3==0){
          a11<-A1[1:(length(A1)/3)]
          a12<-A1[(length(A1)/3+1):(length(A1)/3*2)]
          a13<-A1[(length(A1)/3*2+1):(length(A1))]
        }
        if(length(A1)%%3==1){
          a11<-A1[1:floor(length(A1)/3)]
          a12<-A1[(floor(length(A1)/3)+1):(floor(length(A1)/3)*2+1)]
          a13<-A1[(floor(length(A1)/3)*2+2):length(A1)]
        }
        if(length(A1)%%3==2){
          a11<-A1[1:(ceiling(length(A1)/3))]
          a12<-A1[(ceiling(length(A1)/3)+1):(ceiling(length(A1)/3)+floor(length(A1)/3))]
          a13<-A1[(ceiling(length(A1)/3)+floor(length(A1)/3)+1):length(A1)]
        }
        if(length(B1)%%3==0){
          b11<-B1[1:(length(B1)/3)]
          b12<-B1[(length(B1)/3+1):(length(B1)/3*2)]
          b13<-B1[(length(B1)/3*2+1):(length(B1))]
        }
        if(length(B1)%%3==1){
          b11<-B1[1:floor(length(B1)/3)]
          b12<-B1[(floor(length(B1)/3)+1):(floor(length(B1)/3)*2+1)]
          b13<-B1[(floor(length(B1)/3)*2+2):length(B1)]
        }
        if(length(B1)%%3==2){
          b11<-B1[1:(ceiling(length(B1)/3))]
          b12<-B1[(ceiling(length(B1)/3)+1):(ceiling(length(B1)/3)+floor(length(B1)/3))]
          b13<-B1[(ceiling(length(B1)/3)+floor(length(B1)/3)+1):length(B1)]
        }
        if(length(A2)%%3==0){
          a21<-A2[1:(length(A2)/3)]
          a22<-A2[(length(A2)/3+1):(length(A2)/3*2)]
          a23<-A2[(length(A2)/3*2+1):(length(A2))]
        }
        if(length(A2)%%3==1){
          a21<-A2[1:floor(length(A2)/3)]
          a22<-A2[(floor(length(A2)/3)+1):(floor(length(A2)/3)*2+1)]
          a23<-A2[(floor(length(A2)/3)*2+2):length(A2)]
        }
        if(length(A2)%%3==2){
          a21<-A2[1:(ceiling(length(A2)/3))]
          a22<-A2[(ceiling(length(A2)/3)+1):(ceiling(length(A2)/3)+floor(length(A2)/3))]
          a23<-A2[(ceiling(length(A2)/3)+floor(length(A2)/3)+1):length(A2)]
        }
        medTIMEa11<-median(x[1:length(a11)])
        medTIMEa12<-median(x[(length(a11)+1):(length(a11)+length(a12))])
        medTIMEa13<-median(x[(length(a11)+length(a12)+1):length(A1)])
        medVALUEa11<-median(a11,na.rm=TRUE)
        medVALUEa12<-median(a12,na.rm=TRUE)
        medVALUEa13<-median(a13,na.rm=TRUE)
        medTIMEb11<-median(x[(length(A1)+1):(length(A1)+length(b11))])
        medTIMEb12<-median(x[(length(A1)+length(b11)+1):(length(A1)+length(b11)+length(b12))])
        medTIMEb13<-median(x[(length(A1)+length(b11)+length(b12)+1):(length(A1)+length(B1))])
        medVALUEb11<-median(b11,na.rm=TRUE)
        medVALUEb12<-median(b12,na.rm=TRUE)
        medVALUEb13<-median(b13,na.rm=TRUE)
        medTIMEa21<-median(x[(length(A1)+length(B1)+1):(length(A1)+length(B1)+length(a21))])
        medTIMEa22<-median(x[(length(A1)+length(B1)+length(a21)+1):(length(A1)+length(B1)+length(a21)+length(a22))])
        medTIMEa23<-median(x[(length(A1)+length(B1)+length(a21)+length(a22)+1):(length(A1)+length(B1)+length(A2))])
        medVALUEa21<-median(a21,na.rm=TRUE)
        medVALUEa22<-median(a22,na.rm=TRUE)
        medVALUEa23<-median(a23,na.rm=TRUE)
        slopeA1<-(medVALUEa13-medVALUEa11)/(medTIMEa13-medTIMEa11)
        interceptA1<-(1/3)*((medVALUEa11+medVALUEa12+medVALUEa13)-slopeA1*(medTIMEa11+medTIMEa12+medTIMEa13))
        lines(c(1,length(A1)),c(interceptA1+slopeA1,interceptA1+slopeA1*length(A1)),lty=2)
        slopeB1<-(medVALUEb13-medVALUEb11)/(medTIMEb13-medTIMEb11)
        interceptB1<-(1/3)*((medVALUEb11+medVALUEb12+medVALUEb13)-slopeB1*(medTIMEb11+medTIMEb12+medTIMEb13))
        lines(c(length(A1)+1,(length(A1)+length(B1))),c(interceptB1+slopeB1*(length(A1)+1),interceptB1+slopeB1*(length(A1)+length(B1))),lty=2)
        slopeA2<-(medVALUEa23-medVALUEa21)/(medTIMEa23-medTIMEa21)
        interceptA2<-(1/3)*((medVALUEa21+medVALUEa22+medVALUEa23)-slopeA2*(medTIMEa21+medTIMEa22+medTIMEa23))
        lines(c((length(A1)+length(B1)+1),(length(A1)+length(B1)+length(A2))),c(interceptA2+slopeA2*(length(A1)+length(B1)+1),interceptA2+slopeA2*(length(A1)+length(B1)+length(A2))),lty=2)
        points(medTIMEa11,medVALUEa11,cex=1.2)
        points(medTIMEa12,medVALUEa12,cex=1.2)
        points(medTIMEa13,medVALUEa13,cex=1.2)
        lines(c(medTIMEa11,medTIMEa12),c(medVALUEa11,medVALUEa12),lty=3)
        lines(c(medTIMEa12,medTIMEa13),c(medVALUEa12,medVALUEa13),lty=3)
        points(medTIMEb11,medVALUEb11,cex=1.2)
        points(medTIMEb12,medVALUEb12,cex=1.2)
        points(medTIMEb13,medVALUEb13,cex=1.2)
        lines(c(medTIMEb11,medTIMEb12),c(medVALUEb11,medVALUEb12),lty=3)
        lines(c(medTIMEb12,medTIMEb13),c(medVALUEb12,medVALUEb13),lty=3)
        points(medTIMEa21,medVALUEa21,cex=1.2)
        points(medTIMEa22,medVALUEa22,cex=1.2)
        points(medTIMEa23,medVALUEa23,cex=1.2)
        lines(c(medTIMEa21,medTIMEa22),c(medVALUEa21,medVALUEa22),lty=3)
        lines(c(medTIMEa22,medTIMEa23),c(medVALUEa22,medVALUEa23),lty=3)
        if(design=="ABAB"){
          if(length(B2)%%3==0){
            b21<-B2[1:(length(B2)/3)]
            b22<-B2[(length(B2)/3+1):(length(B2)/3*2)]
            b23<-B2[(length(B2)/3*2+1):(length(B2))]
          }
          if(length(B2)%%3==1){
            b21<-B2[1:floor(length(B2)/3)]
            b22<-B2[(floor(length(B2)/3)+1):(floor(length(B2)/3)*2+1)]
            b23<-B2[(floor(length(B2)/3)*2+2):length(B2)]
          }
          if(length(B2)%%3==2){
            b21<-B2[1:(ceiling(length(B2)/3))]
            b22<-B2[(ceiling(length(B2)/3)+1):(ceiling(length(B2)/3)+floor(length(B2)/3))]
            b23<-B2[(ceiling(length(B2)/3)+floor(length(B2)/3)+1):length(B2)]
          }
          medTIMEb21<-median(x[(length(A1)+length(B1)+length(A2)+1):(length(A1)+length(B1)+length(A2)+length(b21))])
          medTIMEb22<-median(x[(length(A1)+length(B1)+length(A2)+length(b21)+1):(length(A1)+length(B1)+length(A2)+length(b21)+length(b22))])
          medTIMEb23<-median(x[(length(A1)+length(B1)+length(A2)+length(b21)+length(b22)+1):MT])
          medVALUEb21<-median(b21,na.rm=TRUE)
          medVALUEb22<-median(b22,na.rm=TRUE)
          medVALUEb23<-median(b23,na.rm=TRUE)
          slopeB2<-(medVALUEb23-medVALUEb21)/(medTIMEb23-medTIMEb21)
          interceptB2<-(1/3)*((medVALUEb21+medVALUEb22+medVALUEb23)-slopeB2*(medTIMEb21+medTIMEb22+medTIMEb23))
          lines(c((length(A1)++length(B1)+length(A2)+1),MT),c(interceptB2+slopeB2*(length(A1)+length(B1)+length(A2)+1),interceptB2+slopeB2*MT),lty=2)
          points(medTIMEb21,medVALUEb21,cex=1.2)
          points(medTIMEb22,medVALUEb22,cex=1.2)
          points(medTIMEb23,medVALUEb23,cex=1.2)
          lines(c(medTIMEb21,medTIMEb22),c(medVALUEb21,medVALUEb22),lty=3)
          lines(c(medTIMEb22,medTIMEb23),c(medVALUEb22,medVALUEb23),lty=3)
        }
      } 
      if(TREND=="RM3"){
        RM3A1<-numeric()
        for(it in 1:(length(A1)-2)){
          RM3A1<-c(RM3A1,median(A1[it:(it+2)],na.rm=TRUE))
        }
        times3A1<-numeric()
        for(it in 1:(length(A1)-2)){
          times3A1<-c(times3A1,median(x[it:(it+2)]))
        }
        for(it in 1:(length(A1)-2)){
          points(times3A1[it],RM3A1[it],pch=3)
          lines(c(times3A1[it],times3A1[it+1]),c(RM3A1[it],RM3A1[it+1]))
        }
        RM3B1<-numeric()
        for(it in 1:(length(B1)-2)){
          RM3B1<-c(RM3B1,median(B1[it:(it+2)],na.rm=TRUE))
        }
        times3B1<-numeric()
        for(it in (length(A1)+1):(length(A1)+length(B1)-2)){
          times3B1<-c(times3B1,median(x[it:(it+2)]))
        }
        for(it in 1:(length(B1)-2)){
          points(times3B1[it],RM3B1[it],pch=3)
          lines(c(times3B1[it],times3B1[it+1]),c(RM3B1[it],RM3B1[it+1]))
        }
        RM3A2<-numeric()
        for(it in 1:(length(A2)-2)){
          RM3A2<-c(RM3A2,median(A2[it:(it+2)],na.rm=TRUE))
        }
        times3A2<-numeric()
        for(it in (length(A1)+length(B1)+1):(length(A1)+length(B1)+length(A2)-2)){
          times3A2<-c(times3A2,median(x[it:(it+2)]))
        }
        for(it in 1:(length(A2)-2)){
          points(times3A2[it],RM3A2[it],pch=3)
          lines(c(times3A2[it],times3A2[it+1]),c(RM3A2[it],RM3A2[it+1]))
        }
        if(design=="ABAB"){
          RM3B2<-numeric()
          for(it in 1:(length(B2)-2)){
            RM3B2<-c(RM3B2,median(B2[it:(it+2)],na.rm=TRUE))
          }
          times3B2<-numeric()
          for(it in (length(A1)+length(B1)+length(A2)+1):(MT-2)){
            times3B2<-c(times3B2,median(x[it:(it+2)]))
          }
          for(it in 1:(length(B2)-2)){
            points(times3B2[it],RM3B2[it],pch=3)
            lines(c(times3B2[it],times3B2[it+1]),c(RM3B2[it],RM3B2[it+1]))
          }
        }
      }
      if(TREND=="RM5"){
        RM5A1<-numeric()
        for(it in 1:(length(A1)-4)){
          RM5A1<-c(RM5A1,median(A1[it:(it+4)],na.rm=TRUE))
        }
        times5A1<-numeric()
        for(it in 1:(length(A1)-4)){
          times5A1<-c(times5A1,median(x[it:(it+4)]))
        }
        for(it in 1:(length(A1)-4)){
          points(times5A1[it],RM5A1[it],pch=3)
          lines(c(times5A1[it],times5A1[it+1]),c(RM5A1[it],RM5A1[it+1]))
        }
        RM5B1<-numeric()
        for(it in 1:(length(B1)-4)){
          RM5B1<-c(RM5B1,median(B1[it:(it+4)],na.rm=TRUE))
        }
        times5B1<-numeric()
        for(it in (length(A1)+1):(length(A1)+length(B1)-4)){
          times5B1<-c(times5B1,median(x[it:(it+4)]))
        }
        for(it in 1:(length(B1)-4)){
          points(times5B1[it],RM5B1[it],pch=3)
          lines(c(times5B1[it],times5B1[it+1]),c(RM5B1[it],RM5B1[it+1]))
        }
        RM5A2<-numeric()
        for(it in 1:(length(A2)-4)){
          RM5A2<-c(RM5A2,median(A2[it:(it+4)],na.rm=TRUE))
        }
        times5A2<-numeric()
        for(it in (length(A1)+length(B1)+1):(length(A1)+length(B1)+length(A2)-4)){
          times5A2<-c(times5A2,median(x[it:(it+4)]))
        }
        for(it in 1:(length(A2)-4)){
          points(times5A2[it],RM5A2[it],pch=3)
          lines(c(times5A2[it],times5A2[it+1]),c(RM5A2[it],RM5A2[it+1]))
        }
        if(design=="ABAB"){
          RM5B2<-numeric()
          for(it in 1:(length(B2)-4)){
            RM5B2<-c(RM5B2,median(B2[it:(it+4)],na.rm=TRUE))
          }
          times5B2<-numeric()
          for(it in (length(A1)+length(B1)+length(A2)+1):(MT-4)){
            times5B2<-c(times5B2,median(x[it:(it+4)]))
          }
          for(it in 1:(length(B2)-4)){
            points(times5B2[it],RM5B2[it],pch=3)
            lines(c(times5B2[it],times5B2[it+1]),c(RM5B2[it],RM5B2[it+1]))
          }
        }
      }
      if(TREND=="RM42"){
        RM4A1<-numeric()
        for(it in 1:(length(A1)-3)){
          RM4A1<-c(RM4A1,median(A1[it:(it+3)],na.rm=TRUE))
        }
        RM42A1<-numeric()
        for(it in 1:(length(RM4A1)-1)){
          RM42A1<-c(RM42A1,mean(RM4A1[it:(it+1)],na.rm=TRUE))
        }
        times42A1<-numeric()
        for(it in 1:(length(A1)-4)){
          times42A1<-c(times42A1,median(x[it:(it+4)]))
        }
        for(it in 1:(length(A1)-4)){
          points(times42A1[it],RM42A1[it],pch=3)
          lines(c(times42A1[it],times42A1[it+1]),c(RM42A1[it],RM42A1[it+1]))
        }
        RM4B1<-numeric()
        for(it in 1:(length(B1)-3)){
          RM4B1<-c(RM4B1,median(B1[it:(it+3)],na.rm=TRUE))
        }
        RM42B1<-numeric()
        for(it in 1:(length(RM4B1)-1)){
          RM42B1<-c(RM42B1,mean(RM4B1[it:(it+1)],na.rm=TRUE))
        }
        times42B1<-numeric()
        for(it in (length(A1)+1):(length(A1)+length(B1)-4)){
          times42B1<-c(times42B1,median(x[it:(it+4)]))
        }
        for(it in 1:(length(B1)-4)){
          points(times42B1[it],RM42B1[it],pch=3)
          lines(c(times42B1[it],times42B1[it+1]),c(RM42B1[it],RM42B1[it+1]))
        }
        RM4A2<-numeric()
        for(it in 1:(length(A2)-3)){
          RM4A2<-c(RM4A2,median(A2[it:(it+3)],na.rm=TRUE))
        }
        RM42A2<-numeric()
        for(it in 1:(length(RM4A2)-1)){
          RM42A2<-c(RM42A2,mean(RM4A2[it:(it+1)],na.rm=TRUE))
        }
        times42A2<-numeric()
        for(it in (length(A1)+length(B1)+1):(length(A1)+length(B1)+length(A2)-4)){
          times42A2<-c(times42A2,median(x[it:(it+4)]))
        }
        for(it in 1:(length(A2)-4)){
          points(times42A2[it],RM42A2[it],pch=3)
          lines(c(times42A2[it],times42A2[it+1]),c(RM42A2[it],RM42A2[it+1]))
        }
        if(design=="ABAB"){
          RM4B2<-numeric()
          for(it in 1:(length(B2)-3)){
            RM4B2<-c(RM4B2,median(B2[it:(it+3)],na.rm=TRUE))
          }
          RM42B2<-numeric()
          for(it in 1:(length(RM4B2)-1)){
            RM42B2<-c(RM42B2,mean(RM4B2[it:(it+1)],na.rm=TRUE))
          }
          times42B2<-numeric()
          for(it in (length(A1)+length(B1)+length(A2)+1):(MT-4)){
            times42B2<-c(times42B2,median(x[it:(it+4)]))
          }
          for(it in 1:(length(B2)-4)){
            points(times42B2[it],RM42B2[it],pch=3)
            lines(c(times42B2[it],times42B2[it+1]),c(RM42B2[it],RM42B2[it+1]))
          }
        }        
      }
    }
  } 
  
  if(design=="MBD"){
    N<-ncol(data)/2
    par(mfrow=c(N,1))
    for(it in 1:N){
      A<-data[,it*2][data[,(it*2)-1]=="A"]
      B<-data[,it*2][data[,(it*2)-1]=="B"]
      if(TREND=="VLP"){
        if(CL=="mean"){
          CLA<-mean(A,na.rm=TRUE)
          CLB<-mean(B,na.rm=TRUE)
        }
        if(CL=="median"){
          CLA<-median(A,na.rm=TRUE)
          CLB<-median(B,na.rm=TRUE)
        }
        if(CL=="bmed"){
          aa<-sort(A)
          bb<-sort(B)
          if(length(aa)<5){
            CLA<-median(A,na.rm=TRUE)
          }
          if(length(aa)==5|length(aa)==7|length(aa)==9|length(aa)==11){
            CLA<-(aa[ceiling(length(aa)/2)-1]+aa[ceiling(length(aa)/2)]+aa[ceiling(length(aa)/2)+1])/3
          }
          if(length(aa)>=13&length(aa)%%2==1){
            CLA<-(aa[ceiling(length(aa)/2)-2]+aa[ceiling(length(aa)/2)-1]+aa[ceiling(length(aa)/2)]+aa[ceiling(length(aa)/2)+1] +aa[ceiling(length(aa)/2)+2])/5
          }
          if(length(aa)==6|length(aa)==8|length(aa)==10|length(aa)==12){
            CLA<-1/6*aa[length(aa)/2-1]+1/3*aa[length(aa)/2]+1/3*aa[length(aa)/2+1]+1/6*aa[length(aa)/2+2]
          }
          if(length(aa)>13&length(aa)%%2==0){
            CLA<-1/10*aa[length(aa)/2-2]+1/5*aa[length(aa)/2-1]+1/5*aa[length(aa)/2]+1/5*aa[length(aa)/2+1]+1/5*aa[length(aa)/2+2]+1/10*aa[length(aa)/2+3]
          }
          if(length(bb)<5){
            CLB<-median(B,na.rm=TRUE)
          }
          if(length(bb)==5|length(bb)==7|length(bb)==9|length(bb)==11){
            CLB<-(bb[ceiling(length(bb)/2)-1]+bb[ceiling(length(bb)/2)]+bb[ceiling(length(bb)/2)+1])/3
          }
          if(length(bb)>=13&length(bb)%%2==1){
            CLB<-(bb[ceiling(length(bb)/2)-2]+bb[ceiling(length(bb)/2)-1]+bb[ceiling(length(bb)/2)]+bb[ceiling(length(bb)/2)+1]+bb[ceiling(length(bb)/2)+2])/5
          }
          if(length(bb)==6|length(bb)==8|length(bb)==10|length(bb)==12){
            CLB<-1/6*bb[length(bb)/2-1]+1/3*bb[length(bb)/2]+1/3*bb[length(bb)/2+1]+1/6*bb[length(bb)/2+2]
          }
          if(length(bb)>13&length(bb)%%2==0){
            CLB<-1/10*bb[length(bb)/2-2]+1/5*bb[length(bb)/2-1]+1/5*bb[length(bb)/2]+1/5*bb[length(bb)/2+1]+1/5*bb[length(bb)/2+2]+1/10*bb[length(bb)/2+3]
          }
        }
        if(CL=="trimmean"){
          CLA<-mean(A,trim=tr,na.rm=TRUE)	
          CLB<-mean(B,trim=tr,na.rm=TRUE)
        }
        if(CL=="mest"){
          hpsi<-function(x,bend=1.28){
            hpsi<-ifelse(abs(x)<=bend,x,bend*sign(x))
            hpsi
          }
          mest<-function(x,bend=1.28,na.rm=TRUE){
            if(na.rm)x<-x[!is.na(x)]
            if(length(x)==0) return(NA)
            if(mad(x)==0)stop("MAD=0. The M-estimator cannot be computed.")
            y<-(x-median(x))/mad(x)
            A<-sum(hpsi(y,bend))
            B<-length(x[abs(y)<=bend])
            mest<-median(x)+mad(x)*A/B
            repeat{
              y<-(x-mest)/mad(x)
              A<-sum(hpsi(y,bend))
              B<-length(x[abs(y)<=bend])
              newmest<-mest+mad(x)*A/B
              if(abs(newmest-mest) <.0001)break
              mest<-newmest
            }
            mest
          }
          CLA<-mest(A,bend=tr)	
          CLB<-mest(B,bend=tr)
        }  
        residualsA<-A-CLA
        residualsB<-B-CLB
        residuals<-c(residualsA,residualsB)
        if (is.null(xlab)){
          xlab <- ""
        }
        if (is.null(ylab)){
          ylab <- "Residuals"
        }
        plot(x,residuals,xlab=xlab,ylab=ylab,ylim=ylim,type="n")
        lines(c(1,length(residualsA)),c(0,0))
        for(itr in 1:length(residualsA)){
          lines(c(itr,itr),c(residualsA[itr],0))
        }
        lines(c(length(residualsA)+1,length(residuals)),c(0,0))
        for(itr in 1:length(residualsB)){
          lines(c(length(residualsA)+itr,length(residualsA)+itr),c(residualsB[itr],0))
        }
        mtext(labels[1],side=3,at=(sum(data[,(it*2)-1]=="A")+1)/2)
        mtext(labels[2],side=3,at=(sum(data[,(it*2)-1]=="A")+(sum(data[,(it*2)-1]=="B")+1)/2))
      }
      if(TREND=="RTL"|TREND=="SM"|TREND=="LSR"|TREND=="RM3"|TREND=="RM5"|TREND=="RM42"){
        if (is.null(xlab)){
          xlab <- ""
        }
        if (is.null(ylab)){
          ylab <- "Scores"
        }
        plot(x,data[,it*2],xlab=xlab,ylab=ylab,ylim=ylim,pch=16)
        lines(c(sum(data[,(it*2)-1]=="A")+0.5,sum(data[,(it*2)-1]=="A")+0.5),
              c(min(data[,it*2],ylim[1],na.rm=TRUE)-5,max(data[,it*2],ylim[2],na.rm=TRUE)+5),lty=2)
        mtext(labels[1],side=3,at=(sum(data[,(it*2)-1]=="A")+1)/2)
        mtext(labels[2],side=3,at=(sum(data[,(it*2)-1]=="A")+(sum(data[,(it*2)-1]=="B")+1)/2))
        if(TREND=="LSR"){
          xa<-1:length(A)
          interceptLSRa<-coefficients(lm(A~xa,na.action=na.omit))[1]
          slopeLSRa<-coefficients(lm(A~xa,na.action=na.omit))[2]
          xb<-1:length(B)
          interceptLSRb<-coefficients(lm(B~xb,na.action=na.omit))[1]
          slopeLSRb<-coefficients(lm(B~xb,na.action=na.omit))[2]
          lines(c(1,length(A)),c(interceptLSRa+slopeLSRa,interceptLSRa+slopeLSRa*length(A)),lty=2)
          lines(c(length(A)+1,MT),c(interceptLSRb+slopeLSRb,interceptLSRb+slopeLSRb*length(B)),lty=2)
        } 
        if(TREND=="SM"){
          if(length(A)%%2==0){
            a1<-A[1:(length(A)/2)]
            a2<-A[(length(A)/2+1):length(A)]
            medTIMEa1<-median(x[1:length(a1)])
            medTIMEa2<-median(x[(length(a1)+1):length(A)])
            medVALUEa1<-median(a1,na.rm=TRUE)
            medVALUEa2<-median(a2,na.rm=TRUE)
            timeA<-c(medTIMEa1,medTIMEa2)
            valuesA<-c(medVALUEa1,medVALUEa2)
            interceptA<-coefficients(lm(valuesA~timeA))[1]
            slopeA<-coefficients(lm(valuesA~timeA))[2]
            lines(c(1,length(A)),c(interceptA+slopeA,interceptA+slopeA*length(A)),lty=3)
          }
          if(length(A)%%2==1){
            a11<-A[1:ceiling(length(A)/2)]
            a21<-A[(ceiling(length(A)/2)+1):length(A)]
            a12<-A[1:floor(length(A)/2)]
            a22<-A[(floor(length(A)/2)+1):length(A)]
            medTIMEa11<-median(x[1:length(a11)])
            medTIMEa21<-median(x[(length(a11)+1):length(A)])
            medTIMEa12<-median(x[1:length(a12)])
            medTIMEa22<-median(x[(length(a12)+1):length(A)])
            medVALUEa11<-median(a11,na.rm=TRUE)
            medVALUEa21<-median(a21,na.rm=TRUE)
            medVALUEa12<-median(a12,na.rm=TRUE)
            medVALUEa22<-median(a22,na.rm=TRUE)
            timeA1<-c(medTIMEa11,medTIMEa21)
            valuesA1<-c(medVALUEa11,medVALUEa21)
            interceptA1<-coefficients(lm(valuesA1~timeA1))[1]
            slopeA1<-coefficients(lm(valuesA1~timeA1))[2]
            lines(c(1,length(A)),c(interceptA1+slopeA1,interceptA1+slopeA1*length(A)),lty=3)
            timeA2<-c(medTIMEa12,medTIMEa22)
            valuesA2<-c(medVALUEa12,medVALUEa22)
            interceptA2<-coefficients(lm(valuesA2~timeA2))[1]
            slopeA2<-coefficients(lm(valuesA2~timeA2))[2]
            lines(c(1,length(A)),c(interceptA2+slopeA2,interceptA2+slopeA2*length(A)),lty=3)
          }
          if(length(B)%%2==0){
            b1<-B[1:(length(B)/2)]
            b2<-B[(length(B)/2+1):length(B)]
            medTIMEb1<-median(x[(length(A)+1):(length(A)+length(b1))])
            medTIMEb2<-median(x[(length(A)+length(b1)+1):(length(A)+length(B))])
            medVALUEb1<-median(b1,na.rm=TRUE)
            medVALUEb2<-median(b2,na.rm=TRUE)
            timeB<-c(medTIMEb1,medTIMEb2)
            valuesB<-c(medVALUEb1,medVALUEb2)
            interceptB<-coefficients(lm(valuesB~timeB))[1]
            slopeB<-coefficients(lm(valuesB~timeB))[2]
            lines(c(length(A)+1,MT),c(interceptB+slopeB*(length(A)+1),interceptB+slopeB*MT),lty=3)
          }
          if(length(B)%%2==1){
            b11<-B[1:ceiling(length(B)/2)]
            b21<-B[(ceiling(length(B)/2)+1):length(B)]
            b12<-B[1:floor(length(B)/2)]
            b22<-B[(floor(length(B)/2)+1):length(B)]
            medTIMEb11<-median(x[(length(A)+1):(length(A)+length(b11))])
            medTIMEb21<-median(x[(length(A)+length(b11)+1):(length(A)+length(B))])
            medTIMEb12<-median(x[(length(A)+1):(length(A)+length(b12))])
            medTIMEb22<-median(x[(length(A)+length(b12)+1):(length(A)+length(B))])
            medVALUEb11<-median(b11,na.rm=TRUE)
            medVALUEb21<-median(b21,na.rm=TRUE)
            medVALUEb12<-median(b12,na.rm=TRUE)
            medVALUEb22<-median(b22,na.rm=TRUE)
            timeB1<-c(medTIMEb11,medTIMEb21)
            valuesB1<-c(medVALUEb11,medVALUEb21)
            interceptB1<-coefficients(lm(valuesB1~timeB1))[1]
            slopeB1<-coefficients(lm(valuesB1~timeB1))[2]
            lines(c(length(A)+1,MT),c(interceptB1+slopeB1,interceptB1+slopeB1*MT),lty=3)
            timeB2<-c(medTIMEb12,medTIMEb22)
            valuesB2<-c(medVALUEb12,medVALUEb22)
            interceptB2<-coefficients(lm(valuesB2~timeB2))[1]
            slopeB2<-coefficients(lm(valuesB2~timeB2))[2]
            lines(c(length(A)+1,MT),c(interceptB2+slopeB2*(length(A)+1),interceptB2+slopeB2*MT),lty=3)
          }
        } 
        if(TREND=="RTL"){
          if(length(A)%%3==0){
            a1<-A[1:(length(A)/3)]
            a2<-A[(length(A)/3+1):(length(A)/3*2)]
            a3<-A[(length(A)/3*2+1):(length(A))]
          }
          if(length(A)%%3==1){
            a1<-A[1:floor(length(A)/3)]
            a2<-A[(floor(length(A)/3)+1):(floor(length(A)/3)*2+1)]
            a3<-A[(floor(length(A)/3)*2+2):length(A)]
          }
          if(length(A)%%3==2){
            a1<-A[1:(ceiling(length(A)/3))]
            a2<-A[(ceiling(length(A)/3)+1):(ceiling(length(A)/3)+floor(length(A)/3))]
            a3<-A[(ceiling(length(A)/3)+floor(length(A)/3)+1):length(A)]
          }
          if(length(B)%%3==0){
            b1<-B[1:(length(B)/3)]
            b2<-B[(length(B)/3+1):(length(B)/3*2)]
            b3<-B[(length(B)/3*2+1):(length(B))]
          }
          if(length(B)%%3==1){
            b1<-B[1:floor(length(B)/3)]
            b2<-B[(floor(length(B)/3)+1):(floor(length(B)/3)*2+1)]
            b3<-B[(floor(length(B)/3)*2+2):length(B)]
          }
          if(length(B)%%3==2){
            b1<-B[1:(ceiling(length(B)/3))]
            b2<-B[(ceiling(length(B)/3)+1):(ceiling(length(B)/3)+floor(length(B)/3))]
            b3<-B[(ceiling(length(B)/3)+floor(length(B)/3)+1):length(B)]
          }
          medTIMEa1<-median(x[1:length(a1)])
          medTIMEa2<-median(x[(length(a1)+1):(length(a1)+length(a2))])
          medTIMEa3<-median(x[(length(a1)+length(a2)+1):length(A)])
          medVALUEa1<-median(a1,na.rm=TRUE)
          medVALUEa2<-median(a2,na.rm=TRUE)
          medVALUEa3<-median(a3,na.rm=TRUE)
          medTIMEb1<-median(x[(length(A)+1):(length(A)+length(b1))])
          medTIMEb2<-median(x[(length(A)+length(b1)+1):(length(A)+length(b1)+length(b2))])
          medTIMEb3<-median(x[(length(A)+length(b1)+length(b2)+1):MT])
          medVALUEb1<-median(b1,na.rm=TRUE)
          medVALUEb2<-median(b2,na.rm=TRUE)
          medVALUEb3<-median(b3,na.rm=TRUE)
          slopeA<-(medVALUEa3-medVALUEa1)/(medTIMEa3-medTIMEa1)
          interceptA<-(1/3)*((medVALUEa1+medVALUEa2+medVALUEa3)-slopeA*(medTIMEa1+medTIMEa2+medTIMEa3))
          lines(c(1,length(A)),c(interceptA+slopeA,interceptA+slopeA*length(A)),lty=2)
          slopeB<-(medVALUEb3-medVALUEb1)/(medTIMEb3-medTIMEb1)
          interceptB<-(1/3)*((medVALUEb1+medVALUEb2+medVALUEb3)-slopeB*(medTIMEb1+medTIMEb2+medTIMEb3))
          lines(c(length(A)+1,MT),c(interceptB+slopeB*(length(A)+1),interceptB+slopeB*MT),lty=2)
          points(medTIMEa1,medVALUEa1,cex=1.2)
          points(medTIMEa2,medVALUEa2,cex=1.2)
          points(medTIMEa3,medVALUEa3,cex=1.2)
          lines(c(medTIMEa1,medTIMEa2),c(medVALUEa1,medVALUEa2),lty=3)
          lines(c(medTIMEa2,medTIMEa3),c(medVALUEa2,medVALUEa3),lty=3)
          points(medTIMEb1,medVALUEb1,cex=1.2)
          points(medTIMEb2,medVALUEb2,cex=1.2)
          points(medTIMEb3,medVALUEb3,cex=1.2)
          lines(c(medTIMEb1,medTIMEb2),c(medVALUEb1,medVALUEb2),lty=3)
          lines(c(medTIMEb2,medTIMEb3),c(medVALUEb2,medVALUEb3),lty=3)
        } 
        if(TREND=="RM3"){
          RM3a<-numeric()
          for(itr in 1:(length(A)-2)){
            RM3a<-c(RM3a,median(A[itr:(itr+2)],na.rm=TRUE))
          }
          times3A<-numeric()
          for(itr in 1:(length(A)-2)){
            times3A<-c(times3A,median(x[itr:(itr+2)]))
          }
          for(itr in 1:(length(A)-2)){
            points(times3A[itr],RM3a[itr],pch=3)
            lines(c(times3A[itr],times3A[itr+1]),c(RM3a[itr],RM3a[itr+1]))
          }
          RM3b<-numeric()
          for(itr in 1:(length(B)-2)){
            RM3b<-c(RM3b,median(B[itr:(itr+2)],na.rm=TRUE))
          }
          times3B<-numeric()
          for(itr in (length(A)+1):(MT-2)){
            times3B<-c(times3B,median(x[itr:(itr+2)]))
          }
          for(itr in 1:(length(B)-2)){
            points(times3B[itr],RM3b[itr],pch=3)
            lines(c(times3B[itr],times3B[itr+1]),c(RM3b[itr],RM3b[itr+1]))
          }
        }
        if(TREND=="RM5"){
          RM5a<-numeric()
          for(itr in 1:(length(A)-4)){
            RM5a<-c(RM5a,median(A[itr:(itr+4)],na.rm=TRUE))
          }
          times5A<-numeric()
          for(itr in 1:(length(A)-4)){
            times5A<-c(times5A,median(x[itr:(itr+4)]))
          }
          for(itr in 1:(length(A)-4)){
            points(times5A[itr],RM5a[itr],pch=3)
            lines(c(times5A[itr],times5A[itr+1]),c(RM5a[itr],RM5a[itr+1]))
          }
          RM5b<-numeric()
          for(itr in 1:(length(B)-4)){
            RM5b<-c(RM5b,median(B[itr:(itr+4)],na.rm=TRUE))
          }
          times5B<-numeric()
          for(itr in (length(A)+1):(MT-4)){
            times5B<-c(times5B,median(x[itr:(itr+4)]))
          }
          for(itr in 1:(length(B)-4)){
            points(times5B[itr],RM5b[itr],pch=3)
            lines(c(times5B[itr],times5B[itr+1]),c(RM5b[itr],RM5b[itr+1]))
          }
        }
        if(TREND=="RM42"){
          RM4a<-numeric()
          for(itr in 1:(length(A)-3)){
            RM4a<-c(RM4a,median(A[itr:(itr+3)],na.rm=TRUE))
          }
          RM42a<-numeric()
          for(itr in 1:(length(RM4a)-1)){
            RM42a<-c(RM42a,mean(RM4a[itr:(itr+1)],na.rm=TRUE))
          }
          times42A<-numeric()
          for(itr in 1:(length(A)-4)){
            times42A<-c(times42A,median(x[itr:(itr+4)]))
          }
          for(itr in 1:(length(A)-4)){
            points(times42A[itr],RM42a[itr],pch=3)
            lines(c(times42A[itr],times42A[itr+1]),c(RM42a[itr],RM42a[itr+1]))
          }
          RM4b<-numeric()
          for(itr in 1:(length(B)-3)){
            RM4b<-c(RM4b,median(B[itr:(itr+3)],na.rm=TRUE))
          }
          RM42b<-numeric()
          for(itr in 1:(length(RM4b)-1)){
            RM42b<-c(RM42b,mean(RM4b[itr:(itr+1)],na.rm=TRUE))
          }
          times42B<-numeric()
          for(itr in (length(A)+1):(MT-4)){
            times42B<-c(times42B,median(x[itr:(itr+4)]))
          }
          for(itr in 1:(length(B)-4)){
            points(times42B[itr],RM42b[itr],pch=3)
            lines(c(times42B[itr],times42B[itr+1]),c(RM42b[itr],RM42b[itr+1]))
          }
        }
      }       
    }
    if (is.null(xlab)){
      xlab <- "Measurement Times"
    }
    title(xlab=xlab,pch=16)

    par(mfrow=c(1,1))
  }
  
}

Try the SCVA package in your browser

Any scripts or data that you put into this service are public.

SCVA documentation built on Jan. 10, 2020, 1:06 a.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close