
## qqplot for testing the distributions

## for linear model
qq.linear<-function(T, N, r){

T=4
N=300
r=0.1

prob<-seq(.01,.99,.01)

plotname<-paste('QQ-Linear-', T, '-N-', N, '-R-', r, '.pdf', sep='')

pdf(file=plotname, width=5, height=5)

par(mar=c(5, 4, 4, 5) + 0.1)

## complete data
file.comp<-paste('C:\\zzy\\LongPower\\Results\\qqplot\\linear-MCAR','T-',T,'-N',N,'-R',0,'.txt',sep='')
y.comp<-scan(file.comp)
plot(qchisq(prob, 1), quantile(y.comp, prob=prob), ylim=c(0,7), xlim=c(0,7),type='l', lwd=2, ylab='Sample quantiles', xlab='Theoretical quantiles')
#abline(0,1)

## MCAR data
file.mcar<-paste('C:\\zzy\\LongPower\\Results\\qqplot\\linear-MCAR','T-',T,'-N',N,'-R',r,'.txt',sep='')
y.mcar<-scan(file.mcar)
lines(qchisq(prob, 1), quantile(y.mcar, prob=prob), lty=2, lwd=2)

## MAR data
file.mar<-paste('C:\\zzy\\LongPower\\Results\\qqplot\\linear-MAR','T-',T,'-N',N,'-R',r,'.txt',sep='')
y.mar<-scan(file.mar)
lines(qchisq(prob, 1), quantile(y.mar, prob=prob), lty=3, lwd=2)

par(new=T, mar=c(5, 4, 4, 5) + 0.1)
plot(1,1,type='n',ylim=c(0,55), xlim=c(0,7), axes=F, ylab='', xlab='')

## MNAR data
file.mnar<-paste('C:\\zzy\\LongPower\\Results\\qqplot\\linear-MNAR','T-',T,'-N',N,'-R',r,'.txt',sep='')
y.mnar<-scan(file.mnar)
quantile(y.mnar, .99)
lines(qchisq(prob, 1), quantile(y.mnar, prob=prob), lty=4, lwd=2, ylim=c(0,160), xlim=c(0,7))

axis(4)
mtext('MNAR sample quantitles', 4, line=3)

legend('bottomright', c('Complete','MCAR', 'MAR', 'MNAR'), lty=c(1,2,3,4), lwd=c(3,3,3,3))

dev.off()

}

N.v<-c(100,200,300,400,500)
T.v<-3:5
r.v<-c(.1,.3)

 for (T in T.v){
  for (N in N.v){
    for (r in r.v){
         qq.linear(T, N, r)
       }
    }
 }


## Q-Q plot for quadratic model

qq.quad<-function(T, N, r){

T=4
N=300
r=0.3

prob<-seq(.01,.99,.01)

plotname<-paste('QQ-Quad-', T, '-N-', N, '-R-', r, '.pdf', sep='')

pdf(file=plotname, width=5, height=5)

par(mar=c(5, 4, 4, 5) + 0.1)


## complete data
file.comp<-paste('C:\\zzy\\LongPower\\Results\\qqplot\\quad-MCAR','T-',T,'-N',N,'-R',0,'.txt',sep='')
y.comp<-scan(file.comp)
qqplot(qchisq(prob, 1), quantile(y.comp, prob=prob), ylim=c(0,7), xlim=c(0,7),type='l', lwd=2, ylab='Sample quantiles', xlab='Theoretical quantiles')
#abline(0,1)

## MCAR data
file.mcar<-paste('C:\\zzy\\LongPower\\Results\\qqplot\\quad-MCAR','T-',T,'-N',N,'-R',r,'.txt',sep='')
y.mcar<-scan(file.mcar)
lines(qchisq(prob, 1), quantile(y.mcar, prob=prob), lty=2, lwd=2)

## MAR data
file.mar<-paste('C:\\zzy\\LongPower\\Results\\qqplot\\quad-MAR','T-',T,'-N',N,'-R',r,'.txt',sep='')
y.mar<-scan(file.mar)
lines(qchisq(prob, 1), quantile(y.mar, prob=prob), lty=3, lwd=2)

## MNAR data
file.mnar<-paste('C:\\zzy\\LongPower\\Results\\qqplot\\quad-MNAR','T-',T,'-N',N,'-R',r,'.txt',sep='')
y.mnar<-scan(file.mnar)
lines(qchisq(prob, 1), quantile(y.mnar, prob=prob), lty=4, lwd=2)

legend('bottomright', c('Complete','MCAR', 'MAR', 'MNAR'), lty=c(1,2,3,4), lwd=c(3,3,3,3))

dev.off()

}

N.v<-c(100,200,300,400,500)
T.v<-4:5
r.v<-c(.1,.3)

 for (T in T.v){
  for (N in N.v){
    for (r in r.v){
         qq.quad(T, N, r)
       }
    }
 }