# Program to perform four analyses - 
#   1. Sqrt transform, all data
#   2. Log transform, all data
#   3. Sqrt transform, delete Palm Beach
#   4. Log transform, delete Palm Beach
#
# Results of count1-count8 at the end:
# > print(c(c1,d1,w1,a1))
# [1] 0.98715891 0.07622394 0.07648371 0.51404552
# > print(count1)
# [1] 0 0 0 0
# > print(count2)
# [1] 0
# > print(count3)
# [1]  0 58 22  8
# > print(count4)
# [1] 0
# > print(count5)
# [1] 882 708 812 783
# > print(count6)
# [1] 278
# > print(count7)
# [1] 156 433 246 203
# > print(count8)
# [1] 100
fla<-read.table('fldat1.txt',header=T,row.names=NULL)
attach(fla)
postscript(file='fl1.ps',horizo=F)
par(mfrow=c(3,4),bty='l',las=1,oma=c(0,0,0,0))
nsim<-1000
tvot<-bush+gore+brow+nade+harr+hage+mcre+phil+moor
pbush<-100*bush/tvot
pgore<-100*gore/tvot
pbrow<-100*brow/tvot
pnade<-100*nade/tvot
pbuch<-100*buch/tvot
inco<-inco/1000
u<-c(1:49,51:67)
plot(npop[u],pbuch[u],xlab='Total Population',ylab='Buchanan Vote',
xlim=c(min(npop),max(npop)),ylim=c(min(pbuch),max(pbuch)),
main='(a)',log='x')
points(npop[50],pbuch[50],pch='X')
plot(whit[u],pbuch[u],xlab='Percentage White',ylab='Buchanan Vote',
xlim=c(min(whit),max(whit)),ylim=c(min(pbuch),max(pbuch)),
main='(b)')
points(whit[50],pbuch[50],pch='X')
plot(blac[u],pbuch[u],xlab='Percentage Black',ylab='Buchanan Vote',
xlim=c(min(blac),max(blac)),ylim=c(min(pbuch),max(pbuch)),
main='(c)',log='x')
points(blac[50],pbuch[50],pch='X')
plot(hisp[u],pbuch[u],xlab='Percentage Hispanic',ylab='Buchanan Vote',
xlim=c(min(hisp),max(hisp)),ylim=c(min(pbuch),max(pbuch)),
main='(d)',log='x')
points(hisp[50],pbuch[50],pch='X')
plot(o65[u],pbuch[u],xlab='Percentage 65+',ylab='Buchanan Vote',
xlim=c(min(o65),max(o65)),ylim=c(min(pbuch),max(pbuch)),
main='(e)')
points(o65[50],pbuch[50],pch='X')
plot(hsed[u],pbuch[u],xlab='Percentage HS Grad',ylab='Buchanan Vote',
xlim=c(min(hsed),max(hsed)),ylim=c(min(pbuch),max(pbuch)),
main='(f)')
points(hsed[50],pbuch[50],pch='X')
plot(coll[u],pbuch[u],xlab='Percentage Coll Grad',ylab='Buchanan Vote',
xlim=c(min(coll),max(coll)),ylim=c(min(pbuch),max(pbuch)),
main='(g)')
points(coll[50],pbuch[50],pch='X')
plot(inco[u],pbuch[u],xlab='Mean Income',ylab='Buchanan Vote',
xlim=c(min(inco),max(inco)),ylim=c(min(pbuch),max(pbuch)),
main='(h)')
points(inco[50],pbuch[50],pch='X')
plot(pbush[u],pbuch[u],xlab='Bush Vote',ylab='Buchanan Vote',
xlim=c(min(pbush),max(pbush)),ylim=c(min(pbuch),max(pbuch)),
main='(i)')
points(pbush[50],pbuch[50],pch='X')
plot(pgore[u],pbuch[u],xlab='Gore Vote',ylab='Buchanan Vote',
xlim=c(min(pgore),max(pgore)),ylim=c(min(pbuch),max(pbuch)),
main='(j)')
points(pgore[50],pbuch[50],pch='X')
plot(pbrow[u],pbuch[u],xlab='Browne Vote',ylab='Buchanan Vote',
xlim=c(min(pbrow),max(pbrow)),ylim=c(min(pbuch),max(pbuch)),
main='(k)')
points(pbrow[50],pbuch[50],pch='X')
plot(pnade[u],pbuch[u],xlab='Nader Vote',ylab='Buchanan Vote',
xlim=c(min(pnade),max(pnade)),ylim=c(min(pbuch),max(pbuch)),
main='(l)')
points(pnade[50],pbuch[50],pch='X')
#mtext(' Gain Functions of Various Filters',
#line=1,side=3,cex=1.5,outer=T)
dev.off()
lpop<-log(npop)
lblac<-log(blac)
lhisp<-log(hisp)
lpop1<-tvot*lpop 
whit1<-tvot*whit 
lblac1<-tvot*lblac 
lhisp1<-tvot*lhisp 
o651<-tvot*o65 
hsed1<-tvot*hsed 
coll1<-tvot*coll 
inco1<-tvot*inco 
pbush1<-tvot*pbush 
pbrow1<-tvot*pbrow 
pnade1<-tvot*pnade 
lpop2<-sqrt(tvot)*lpop 
whit2<-sqrt(tvot)*whit 
lblac2<-sqrt(tvot)*lblac 
lhisp2<-sqrt(tvot)*lhisp 
o652<-sqrt(tvot)*o65 
hsed2<-sqrt(tvot)*hsed 
coll2<-sqrt(tvot)*coll 
inco2<-sqrt(tvot)*inco 
pbush2<-sqrt(tvot)*pbush 
pbrow2<-sqrt(tvot)*pbrow 
pnade2<-sqrt(tvot)*pnade 
v<-u
#v<-1:67
#v<-c(1:6,8:37,39:49,51:67)
# First do untransformed regression
y<-buch
nreg<-lm(y~lpop1+whit1+lblac1+lhisp1+o651+hsed1+inco1+pbrow1,subset=v)
rsd<-abs(studres(nreg))
rsd<-abs(nreg$resid)
rsd2<-rsd^2
prd<-predict(nreg)
# Residual plots
postscript(file='fl4a.ps',horizo=F)
par(mfrow=c(2,3),bty='l',las=1,oma=c(0,0,0,0))
plot(prd,rsd,xlab='Fitted Value',ylab='Absolute Residual',log='x',main='(a)')
sd1<-lowess(log(prd),rsd2,f=1.0)
lines(exp(sd1$x),sqrt(sd1$y))
plot(tvot[v],rsd,xlab='Total Vote',ylab='Absolute Residual',log='x',main='(b)')
sd1<-lowess(log(tvot[v]),rsd2,f=1.0)
lines(exp(sd1$x),sqrt(sd1$y))
plot(pbuch[v],rsd,xlab='Buchanan Vote',ylab='Absolute Residual',log='x',main='(c)')
sd1<-lowess(log(pbuch[v]),rsd2,f=1.0)
lines(exp(sd1$x),sqrt(sd1$y))
dev.off()
# Log-p transformation
y<-log(pbuch)
nreg<-lm(y~lpop+lhisp+hsed+inco+pbush,subset=v)
rsd<-abs(studres(nreg))
rsd<-abs(nreg$resid)
rsd2<-rsd*rsd
prd<-predict(nreg)
# Residual plots
postscript(file='fl4c.ps',horizo=F)
par(mfrow=c(2,3),bty='l',las=1,oma=c(0,0,0,0))
plot(prd,rsd,xlab='Fitted Value',ylab='Absolute Residual',main='(a)')
sd1<-lowess(prd,rsd2,f=1.0)
lines(sd1$x,sqrt(sd1$y))
plot(tvot[v],rsd,xlab='Total Vote',ylab='Absolute Residual',log='x',main='(b)')
sd1<-lowess(log(tvot[v]),rsd2,f=1.0)
lines(exp(sd1$x),sqrt(sd1$y))
plot(pbuch[v],rsd,xlab='Buchanan Vote',ylab='Absolute Residual',log='x',main='(c)')
sd1<-lowess(log(pbuch[v]),rsd2,f=1.0)
lines(exp(sd1$x),sqrt(sd1$y))
dev.off()
# Square root transformation
y<-sqrt(buch)
nreg<-lm(y~lpop2+whit2+lhisp2+hsed2+inco2+pbush2+pbrow2,subset=v)
rsd<-abs(nreg$resid)
rsd<-abs(studres(nreg))
prd<-predict(nreg)
# Residual plots
postscript(file='fl4b.ps',horizo=F)
par(mfrow=c(2,3),bty='l',las=1,oma=c(0,0,0,0))
plot(prd,rsd,xlab='Fitted Value',ylab='Absolute Residual',log='x',main='(a)')
sd1<-lowess(log(prd),rsd2,f=1.0)
lines(exp(sd1$x),sqrt(sd1$y))
plot(tvot[v],rsd,xlab='Total Vote',ylab='Absolute Residual',log='x',main='(b)')
sd1<-lowess(log(tvot[v]),rsd2,f=1.0)
lines(exp(sd1$x),sqrt(sd1$y))
plot(pbuch[v],rsd,xlab='Buchanan Vote',ylab='Absolute Residual',log='x',main='(c)')
sd1<-lowess(log(pbuch[v]),rsd2,f=1.0)
lines(exp(sd1$x),sqrt(sd1$y))
dev.off()
#
# deletion diagnostics for square root regression
#
count2<-0
#v<-u
v<-1:67
#v<-c(1:6,8:37,39:49,51:67)
y<-sqrt(buch)
nreg<-lm(y~lpop2+whit2+lhisp2+hsed2+inco2+pbush2+pbrow2,subset=v)
nn<-length(y[v])
nn2<-length(y)
df1<-sort(abs(dffits(nreg)))
rs1<-sort(abs(studres(nreg)))
xd1<-qnorm((((nn+2):(2*nn+1))-0.375)/(2*nn+1.25))
y2<-rnorm(nn2,mean=0,sd=1)
nreg1<-lm(y2~lpop2+whit2+lhisp2+hsed2+inco2+pbush2+pbrow2,subset=v)
dfsim<-sort(abs(dffits(nreg1)))
rssim<-sort(abs(studres(nreg1)))
if(max(abs(studres(nreg1)))>max(abs(studres(nreg))))count2<-count2+1
for(j in 2:nsim){
print(j)
y2<-rnorm(nn2,mean=0,sd=1)
nreg1<-lm(y2~lpop2+whit2+lhisp2+hsed2+inco2+pbush2+pbrow2,subset=v)
if(max(abs(studres(nreg1)))>max(abs(studres(nreg))))count2<-count2+1
dfsim<-c(dfsim,sort(abs(dffits(nreg1))))
rssim<-c(rssim,sort(abs(studres(nreg1))))
}
dfsim<-matrix(dfsim,ncol=nn,byrow=T)
rssim<-matrix(rssim,ncol=nn,byrow=T)
df2<-rep(0,nn)
df3<-rep(0,nn)
rs2<-rep(0,nn)
rs3<-rep(0,nn)
for(i in 1:nn){
t1<-sort(dfsim[,i])
df2[i]<-t1[0.05*nsim]
df3[i]<-t1[0.95*nsim]
t1<-sort(rssim[,i])
rs2[i]<-t1[0.05*nsim]
rs3[i]<-t1[0.95*nsim]
}
postscript(file='fl2.ps',horizo=F)
par(mfrow=c(2,2),bty='l',las=1,oma=c(0,0,0,0))
ymin<-min(c(rs1,rs2,rs3))
ymax<-max(c(rs1,rs2,rs3))
plot(xd1,rs1,xlab='Expected Value',ylab='Observed Value',
main='(a)',ylim=c(ymin,ymax))
points(xd1,rs2,pch='-')
points(xd1,rs3,pch='-')
ymin<-min(c(df1,df2,df3))
ymax<-max(c(df1,df2,df3))
plot(xd1,df1,xlab='Observed Value',ylab='Expected Value',
main='(b)',ylim=c(ymin,ymax))
points(xd1,df2,pch='-')
points(xd1,df3,pch='-')
dev.off()
# normality diagnostics
postscript(file='fl3.ps',horizo=F)
par(mfrow=c(2,1),bty='l',las=1,oma=c(0,0,0,0))
q1<-qqnorm(nreg$resid,plot=F)
plot(q1$x,q1$y/sqrt(var(q1$y)),
xlab='Expected value',ylab='Observed value')
lines(c(-1000,1000),c(-1000,1000))
dev.off()
c1<-cor(q1$x,q1$y)
u1<-pnorm(sort(nreg$resid),mean=0,sd=summary(nreg)$sigma)
d1<-max(c((1:nn)/nn-u1,u1-(0:(nn-1))/nn))
w1<-sum((u1-((1:nn)-0.5)/nn)^2)+1/(12*nn)
a1<--sum((2*(1:nn)-1)*log(u1)+(2*nn+1-2*(1:nn))*log(1-u1))/nn-nn
count1<-rep(0,4)
for(j in 1:nsim){
y2<-rnorm(nn2,mean=0,sd=1)
nreg1<-lm(y2~lpop2+whit2+lhisp2+hsed2+inco2+pbush2+pbrow2,subset=v)
q2<-qqnorm(nreg1$resid,plot=F)
c2<-cor(q2$x,q2$y)
if(c2<c1)count1[1]<-count1[1]+1
u2<-pnorm(sort(nreg1$resid),mean=0,sd=summary(nreg1)$sigma)
d2<-max(c((1:nn)/nn-u2,u2-(0:(nn-1))/nn))
w2<-sum((u2-((1:nn)-0.5)/nn)^2)+1/(12*nn)
a2<--sum((2*(1:nn)-1)*log(u2)+(2*nn+1-2*(1:nn))*log(1-u2))/nn-nn
if(d2>d1)count1[2]<-count1[2]+1
if(w2>w1)count1[3]<-count1[3]+1
if(a2>a1)count1[4]<-count1[4]+1
print(c(j,c2,d2,w2,a2))
}
print(c(c1,d1,w1,a1))
#
# deletion diagnostics for log regression
#
count4<-0
#v<-u
v<-1:67
#v<-c(1:6,8:37,39:49,51:67)
y<-log(pbuch)
nreg<-lm(y~lpop+lhisp+hsed+inco+pbush,subset=v)
nn<-length(y[v])
nn2<-length(y)
df1<-sort(abs(dffits(nreg)))
rs1<-sort(abs(studres(nreg)))
xd1<-qnorm((((nn+2):(2*nn+1))-0.375)/(2*nn+1.25))
y2<-rnorm(nn2,mean=0,sd=1)
nreg1<-lm(y2~lpop+lhisp+hsed+inco+pbush,subset=v)
dfsim<-sort(abs(dffits(nreg1)))
rssim<-sort(abs(studres(nreg1)))
if(max(abs(studres(nreg1)))>max(abs(studres(nreg))))count4<-count4+1
for(j in 2:nsim){
print(j)
y2<-rnorm(nn2,mean=0,sd=1)
nreg1<-lm(y2~lpop+lhisp+hsed+inco+pbush,subset=v)
if(max(abs(studres(nreg1)))>max(abs(studres(nreg))))count4<-count4+1
dfsim<-c(dfsim,sort(abs(dffits(nreg1))))
rssim<-c(rssim,sort(abs(studres(nreg1))))
}
dfsim<-matrix(dfsim,ncol=nn,byrow=T)
rssim<-matrix(rssim,ncol=nn,byrow=T)
df2<-rep(0,nn)
df3<-rep(0,nn)
rs2<-rep(0,nn)
rs3<-rep(0,nn)
for(i in 1:nn){
t1<-sort(dfsim[,i])
df2[i]<-t1[0.05*nsim]
df3[i]<-t1[0.95*nsim]
t1<-sort(rssim[,i])
rs2[i]<-t1[0.05*nsim]
rs3[i]<-t1[0.95*nsim]
}
postscript(file='fl2a.ps',horizo=F)
par(mfrow=c(2,2),bty='l',las=1,oma=c(0,0,0,0))
ymin<-min(c(rs1,rs2,rs3))
ymax<-max(c(rs1,rs2,rs3))
plot(xd1,rs1,xlab='Expected Value',ylab='Observed Value',
main='(a)',ylim=c(ymin,ymax))
points(xd1,rs2,pch='-')
points(xd1,rs3,pch='-')
ymin<-min(c(df1,df2,df3))
ymax<-max(c(df1,df2,df3))
plot(xd1,df1,xlab='Observed Value',ylab='Expected Value',
main='(b)',ylim=c(ymin,ymax))
points(xd1,df2,pch='-')
points(xd1,df3,pch='-')
dev.off()
# normality diagnostics
postscript(file='fl3a.ps',horizo=F)
par(mfrow=c(2,1),bty='l',las=1,oma=c(0,0,0,0))
q1<-qqnorm(nreg$resid,plot=F)
plot(q1$x,q1$y/sqrt(var(q1$y)),
xlab='Expected value',ylab='Observed value')
lines(c(-1000,1000),c(-1000,1000))
dev.off()
c1<-cor(q1$x,q1$y)
u1<-pnorm(sort(nreg$resid),mean=0,sd=summary(nreg)$sigma)
d1<-max(c((1:nn)/nn-u1,u1-(0:(nn-1))/nn))
w1<-sum((u1-((1:nn)-0.5)/nn)^2)+1/(12*nn)
a1<--sum((2*(1:nn)-1)*log(u1)+(2*nn+1-2*(1:nn))*log(1-u1))/nn-nn
count3<-rep(0,4)
for(j in 1:nsim){
y2<-rnorm(nn2,mean=0,sd=1)
nreg1<-lm(y2~lpop+lhisp+hsed+inco+pbush,subset=v)
q2<-qqnorm(nreg1$resid,plot=F)
c2<-cor(q2$x,q2$y)
if(c2<c1)count3[1]<-count3[1]+1
u2<-pnorm(sort(nreg1$resid),mean=0,sd=summary(nreg1)$sigma)
d2<-max(c((1:nn)/nn-u2,u2-(0:(nn-1))/nn))
w2<-sum((u2-((1:nn)-0.5)/nn)^2)+1/(12*nn)
a2<--sum((2*(1:nn)-1)*log(u2)+(2*nn+1-2*(1:nn))*log(1-u2))/nn-nn
if(d2>d1)count3[2]<-count3[2]+1
if(w2>w1)count3[3]<-count3[3]+1
if(a2>a1)count3[4]<-count3[4]+1
print(c(j,c2,d2,w2,a2))
}
print(c(c1,d1,w1,a1))
# 
# Repeat everything with Palm Beach deleted
#
#
# deletion diagnostics for square root regression
#
count6<-0
v<-u
#v<-1:67
#v<-c(1:6,8:37,39:49,51:67)
y<-sqrt(buch)
nreg<-lm(y~lpop2+whit2+lhisp2+hsed2+inco2+pbush2+pbrow2,subset=v)
df0<-dffits(nreg)
nn<-length(y[v])
nn2<-length(y)
df1<-sort(abs(dffits(nreg)))
rs1<-sort(abs(studres(nreg)))
xd1<-qnorm((((nn+2):(2*nn+1))-0.375)/(2*nn+1.25))
y2<-rnorm(nn2,mean=0,sd=1)
nreg1<-lm(y2~lpop2+whit2+lhisp2+hsed2+inco2+pbush2+pbrow2,subset=v)
dfsim<-sort(abs(dffits(nreg1)))
rssim<-sort(abs(studres(nreg1)))
if(max(abs(studres(nreg1)))>max(abs(studres(nreg))))count6<-count6+1
for(j in 2:nsim){
print(j)
y2<-rnorm(nn2,mean=0,sd=1)
nreg1<-lm(y2~lpop2+whit2+lhisp2+hsed2+inco2+pbush2+pbrow2,subset=v)
if(max(abs(studres(nreg1)))>max(abs(studres(nreg))))count6<-count6+1
dfsim<-c(dfsim,sort(abs(dffits(nreg1))))
rssim<-c(rssim,sort(abs(studres(nreg1))))
}
dfsim<-matrix(dfsim,ncol=nn,byrow=T)
rssim<-matrix(rssim,ncol=nn,byrow=T)
df2<-rep(0,nn)
df3<-rep(0,nn)
rs2<-rep(0,nn)
rs3<-rep(0,nn)
for(i in 1:nn){
t1<-sort(dfsim[,i])
df2[i]<-t1[0.05*nsim]
df3[i]<-t1[0.95*nsim]
t1<-sort(rssim[,i])
rs2[i]<-t1[0.05*nsim]
rs3[i]<-t1[0.95*nsim]
}
postscript(file='fl2b.ps',horizo=F)
par(mfrow=c(2,2),bty='l',las=1,oma=c(0,0,0,0))
ymin<-min(c(rs1,rs2,rs3))
ymax<-max(c(rs1,rs2,rs3))
plot(xd1,rs1,xlab='Expected Value',ylab='Observed Value',
main='(a)',ylim=c(ymin,ymax))
points(xd1,rs2,pch='-')
points(xd1,rs3,pch='-')
ymin<-min(c(df1,df2,df3))
ymax<-max(c(df1,df2,df3))
plot(xd1,df1,xlab='Observed Value',ylab='Expected Value',
main='(b)',ylim=c(ymin,ymax))
points(xd1,df2,pch='-')
points(xd1,df3,pch='-')
dev.off()
# normality diagnostics
postscript(file='fl3b.ps',horizo=F)
par(mfrow=c(2,1),bty='l',las=1,oma=c(0,0,0,0))
q1<-qqnorm(nreg$resid,plot=F)
plot(q1$x,q1$y/sqrt(var(q1$y)),
xlab='Expected value',ylab='Observed value')
lines(c(-1000,1000),c(-1000,1000))
dev.off()
c1<-cor(q1$x,q1$y)
u1<-pnorm(sort(nreg$resid),mean=0,sd=summary(nreg)$sigma)
d1<-max(c((1:nn)/nn-u1,u1-(0:(nn-1))/nn))
w1<-sum((u1-((1:nn)-0.5)/nn)^2)+1/(12*nn)
a1<--sum((2*(1:nn)-1)*log(u1)+(2*nn+1-2*(1:nn))*log(1-u1))/nn-nn
count5<-rep(0,4)
for(j in 1:nsim){
y2<-rnorm(nn2,mean=0,sd=1)
nreg1<-lm(y2~lpop2+whit2+lhisp2+hsed2+inco2+pbush2+pbrow2,subset=v)
q2<-qqnorm(nreg1$resid,plot=F)
c2<-cor(q2$x,q2$y)
if(c2<c1)count5[1]<-count5[1]+1
u2<-pnorm(sort(nreg1$resid),mean=0,sd=summary(nreg1)$sigma)
d2<-max(c((1:nn)/nn-u2,u2-(0:(nn-1))/nn))
w2<-sum((u2-((1:nn)-0.5)/nn)^2)+1/(12*nn)
a2<--sum((2*(1:nn)-1)*log(u2)+(2*nn+1-2*(1:nn))*log(1-u2))/nn-nn
if(d2>d1)count5[2]<-count5[2]+1
if(w2>w1)count5[3]<-count5[3]+1
if(a2>a1)count5[4]<-count5[4]+1
print(c(j,c2,d2,w2,a2))
}
print(c(c1,d1,w1,a1))
#
# deletion diagnostics for log regression
#
count8<-0
v<-u
#v<-1:67
#v<-c(1:6,8:37,39:49,51:67)
y<-log(pbuch)
nreg<-lm(y~lpop+lhisp+hsed+inco+pbush,subset=v)
#
# Plot dffits against population number
#
postscript(file='fl5.ps',horizo=F)
par(mfrow=c(2,2),bty='l',las=1,oma=c(0,0,0,0))
plot(tvot[v],df0,xlab='Total Votes',ylab='DFFITS',main='(a)')
plot(tvot[v],dffits(nreg),xlab='Total Votes',ylab='DFFITS',main='(b)')
dev.off()
nn<-length(y[v])
nn2<-length(y)
df1<-sort(abs(dffits(nreg)))
rs1<-sort(abs(studres(nreg)))
xd1<-qnorm((((nn+2):(2*nn+1))-0.375)/(2*nn+1.25))
y2<-rnorm(nn2,mean=0,sd=1)
nreg1<-lm(y2~lpop+lhisp+hsed+inco+pbush,subset=v)
dfsim<-sort(abs(dffits(nreg1)))
rssim<-sort(abs(studres(nreg1)))
if(max(abs(studres(nreg1)))>max(abs(studres(nreg))))count8<-count8+1
for(j in 2:nsim){
print(j)
y2<-rnorm(nn2,mean=0,sd=1)
nreg1<-lm(y2~lpop+lhisp+hsed+inco+pbush,subset=v)
if(max(abs(studres(nreg1)))>max(abs(studres(nreg))))count8<-count8+1
dfsim<-c(dfsim,sort(abs(dffits(nreg1))))
rssim<-c(rssim,sort(abs(studres(nreg1))))
}
dfsim<-matrix(dfsim,ncol=nn,byrow=T)
rssim<-matrix(rssim,ncol=nn,byrow=T)
df2<-rep(0,nn)
df3<-rep(0,nn)
rs2<-rep(0,nn)
rs3<-rep(0,nn)
for(i in 1:nn){
t1<-sort(dfsim[,i])
df2[i]<-t1[0.05*nsim]
df3[i]<-t1[0.95*nsim]
t1<-sort(rssim[,i])
rs2[i]<-t1[0.05*nsim]
rs3[i]<-t1[0.95*nsim]
}
postscript(file='fl2c.ps',horizo=F)
par(mfrow=c(2,2),bty='l',las=1,oma=c(0,0,0,0))
ymin<-min(c(rs1,rs2,rs3))
ymax<-max(c(rs1,rs2,rs3))
plot(xd1,rs1,xlab='Expected Value',ylab='Observed Value',
main='(a)',ylim=c(ymin,ymax))
points(xd1,rs2,pch='-')
points(xd1,rs3,pch='-')
ymin<-min(c(df1,df2,df3))
ymax<-max(c(df1,df2,df3))
plot(xd1,df1,xlab='Observed Value',ylab='Expected Value',
main='(b)',ylim=c(ymin,ymax))
points(xd1,df2,pch='-')
points(xd1,df3,pch='-')
dev.off()
# normality diagnostics
postscript(file='fl3c.ps',horizo=F)
par(mfrow=c(2,1),bty='l',las=1,oma=c(0,0,0,0))
q1<-qqnorm(nreg$resid,plot=F)
plot(q1$x,q1$y/sqrt(var(q1$y)),
xlab='Expected value',ylab='Observed value')
lines(c(-1000,1000),c(-1000,1000))
dev.off()
c1<-cor(q1$x,q1$y)
u1<-pnorm(sort(nreg$resid),mean=0,sd=summary(nreg)$sigma)
d1<-max(c((1:nn)/nn-u1,u1-(0:(nn-1))/nn))
w1<-sum((u1-((1:nn)-0.5)/nn)^2)+1/(12*nn)
a1<--sum((2*(1:nn)-1)*log(u1)+(2*nn+1-2*(1:nn))*log(1-u1))/nn-nn
count7<-rep(0,4)
for(j in 1:nsim){
y2<-rnorm(nn2,mean=0,sd=1)
nreg1<-lm(y2~lpop+lhisp+hsed+inco+pbush,subset=v)
q2<-qqnorm(nreg1$resid,plot=F)
c2<-cor(q2$x,q2$y)
if(c2<c1)count7[1]<-count7[1]+1
u2<-pnorm(sort(nreg1$resid),mean=0,sd=summary(nreg1)$sigma)
d2<-max(c((1:nn)/nn-u2,u2-(0:(nn-1))/nn))
w2<-sum((u2-((1:nn)-0.5)/nn)^2)+1/(12*nn)
a2<--sum((2*(1:nn)-1)*log(u2)+(2*nn+1-2*(1:nn))*log(1-u2))/nn-nn
if(d2>d1)count7[2]<-count7[2]+1
if(w2>w1)count7[3]<-count7[3]+1
if(a2>a1)count7[4]<-count7[4]+1
print(c(j,c2,d2,w2,a2))
}
print(c(c1,d1,w1,a1))
print(count1)
print(count2)
print(count3)
print(count4)
print(count5)
print(count6)
print(count7)
print(count8)