#SAS and R code  for Robust Multivariate Analysis homework problems;

#Copy and paste code one problem part at a time.


#The easiest way to get mpack.txt and mrobdata.txt is to copy and paste the following two commands into R.

source("http://parker.ad.siu.edu/Olive/mpack.txt")
source("http://parker.ad.siu.edu/Olive/mrobdata.txt")


#If mpack.txt and mrobdata.txt are on J drive, get rpack and robdata with
#the following commands:

source("J:/mpack.txt") 
source("J:/mrobdata.txt")

#####for Linux use left most port for flash drive

source("/media/190A-4001/mpack.txt") 
source("/media/190A-4001/mrobdata.txt")


##Problem 1.3

#a)

height <- rnorm(87, mean=1692, sd = 65)
height[61:65] <- 19.0 
#could use height <- buxy 

#b)

cci(height)

#c)

medci(height)


##Problem 2.10, HW1 D)

n3x <- matrix(rnorm(300),nrow=100,ncol=3)
ln3x <- exp(n3x)
library(MASS)

pairs(n3x)
pairs(ln3x)


##Problem 4.1, needs mpack

#b)

simx2 <- matrix(rnorm(200),nrow=100,ncol=2)
outx2 <- matrix(10 + rnorm(80),nrow=40,ncol=2)
outx2 <- rbind(outx2,simx2)
maha(outx2)


##Problem 4.2, needs mpack and get outx2 from Problem 4.1

library(MASS)
rmaha(outx2)

#Problem 4.3, needs mpack

rcovsim(100)

rcovsim(100)

rcovsim(100)


###R commands for 4.4, 4.5, and 4.6

source("J:/mpack.txt")  # assumes flash drive is on drive J
source("J:/mrobdata.txt") 

library(MASS)


##Problem 4.4, needs mpack and mrobdata

#a)

library(MASS)
ddcomp(buxx)  #left click on the outliers and right click Stop 4 times 

#b)

ddcomp(cbrainx)  #left click on the outliers and right click Stop 4 times 

#c)

ddcomp(museum[,-1])  #left click on the outliers and right click Stop 4 times 


##Problem 4.5 R commands, needs mpack

#repeat the following command several times

concmv() #right click Stop 5 times


##Problem 4.6, needs mpack

ddmv(p=2,gam=.4)     #right click Stop 5 times
ddmv(p=4,gam=.35)
ddmv(p=4,gam=.30)


##Problem 4.7, needs mpack, mrobdata

cmba2(buxx) #right click Stop on the plot 7 times

cmba2(cbrainx) #right click Stop 7 times

cmba2(museum[,-1]) # right click Stop 7 times


##Problem 4.8, HW4 F), needs mpack

#a)

covesim(n=100,p=4)

#b) 

covesim(n=100,p=4,outliers=1,pm=15)


##Problem 4.9, needs mrobdata

library(MASS)

z <- cbind(cbrainx,cbrainy)
z<-z[,-c(2,4,10)]
dim(z)
#Get the log determinant of the FMCD estimator
#applied to a permutation of the data. If the MCD
#estimator was computed, the log determinant would
#be a minimum and would not change under a permutation
#of the rows of the data matrix.
cov.mcd(z[sample(1:267),])$crit
cov.mcd(z[sample(1:267),])$crit
cov.mcd(z[sample(1:267),])$crit
cov.mcd(z[sample(1:267),])$crit
cov.mcd(z[sample(1:267),])$crit
cov.mcd(z[sample(1:267),])$crit
cov.mcd(z[sample(1:267),])$crit


###HW commands for 5.2-5.5

source("G:/mpack.txt")

library(MASS)


##Problem 5.2, needs mpack

#b)

library(MASS)
ddsim(n=20,p=2)
ddsim(n=30,p=2)


##Problem 5.3, needs mpack

library(MASS)
corrsim(n=20,p=2,nruns=10)
corrsim(n=30,p=2,nruns=10)


##Problem 5.4, needs mpack

#b)

library(MASS)
n <- 400
p <- 3
eps <- 0.4
x <- matrix(rnorm(n * p), ncol = p, nrow = n)
zu <- runif(n)
x[zu < eps,] <- x[zu < eps,]*5

#c)

ddplot(x) #right click Stop once on the plot


##Problem 5.5, needs mpack

#b)

#plot for classical covering ellipsoid
simx2 <- matrix(rnorm(200),nrow=100,ncol=2)
outx2 <- matrix(10 + rnorm(80),nrow=40,ncol=2)
outx2 <- rbind(outx2,simx2)
ellipse(outx2)

#plot for RMVN covering ellipsoid
zout <- covrmvn(outx2)
ellipse(outx2,center=zout$center,cov=zout$cov)


##Problem 5.6, needs mpack 

#c)
# get x using Problem 5.4


mplot(x)  #right click Stop  once


##Problem 5.7, needs mpack

#c)

# get x using Problem 5.4

wddplot(x)  


##Problem 5.8, HW4 G) needs mpack

source("G:/mrobdata.txt")
ddplot4(buxx,alpha=0.2) #right click Stop once


##Problem 5.9, needs mpack

predsim()


##Problem 5.10, needs mpack


#a) diagnostic for sphericity

x <- matrix(rnorm(1000), nrow=100, ncol=10)
x <- 10*x
MDSQ <- mahalanobis(x,mean(x),var(x))
DSQ <- mahalanobis(x,mean(x),diag(10))
plot(MDSQ,DSQ)
#abline(0,100)


##Problem 6.9, needs mpack

pcasim()


##Problem 6.10, needs mpack

#a) 

library(MASS)
help(cpus)
z <- cpus[,-1]
z <- z[,-c(3,4)]
z <- z + 1
zz<-z
pairs(zz)

#b)

zz <- log(z)
pairs(zz)

ddplot(zz) #right click Stop on the plot once

#c)

#uses the log transformation on the 1st 4 variables
#currently variables 5 and 6 also use the log transformation
zz[,-c(5,6)] <- log(z[,-c(5,6)])

#possible transformations for the 5th and 6th variables are below
zz[,6] <- z[,6]^{-1/3}
pairs(zz)

ddplot(zz) #right click Stop once

zz[,6] <- 1/sqrt(z[,6])
pairs(zz)
ddplot(zz) #right click Stop once

zz[,6] <- 1/z[,6]
pairs(zz)

ddplot(zz) #right click Stop once

zz[,5] <- log(z[,5])
zz[,6] <- log(z[,6])
pairs(zz)

ddplot(zz) #right click Stop once

zz[,5] <- z[,5]^{-1/3}
zz[,6] <- z[,6]^{-1/3}
pairs(zz)

ddplot(zz) #right click Stop once

zz[,5] <- 1/sqrt(z[,5])
zz[,6] <- 1/sqrt(z[,6])
pairs(zz)

ddplot(zz) #right click Stop once

zz[,5] <- 1/z[,5]
zz[,6] <- 1/z[,6]
pairs(zz)

ddplot(zz) #right click Stop once

#d) 

pairs(zz)

ddplot(zz) #right click Stop once

#e)

prcomp(zz, scale = T)

#f) 

rprcomp(zz)


##Problem 6.11

#a) 

out<-prcomp(USArrests)
summary(out)

out<-rprcomp(USArrests,corr=F)
summary(out$out)

#b)

out<-prcomp(USArrests,scale=T)
summary(out)

out<-rprcomp(USArrests)
summary(out$out)


##Problem 6.12, needs mpack and mrobdata
##scree plots and biplots 
#robust biplot only shows data from U

#a)

out <- prcomp(buxx,corr=T)
plot(out)

biplot(out)

#b)

out <- rprcomp(buxx)
plot(out$out)

biplot(out$out)

#Problem 6.13

#a)

out <- prcomp(bodfat,corr=T)
plot(out)

biplot(out)

#b)

rout <- rprcomp(bodfat)
plot(rout$out)

biplot(rout$out)


##Problem 7.6, needs rcancor

m <- 3
q <- 5
ni <- 500

w <- matrix(rnorm(ni*m), ni, m)
y <- matrix(rnorm(ni*q), ni, q)
cancor(w,y)$cor[1]
rcancor(w,y)$out$cor[1]

ni <- 1000 
w <- matrix(rnorm(ni*m), ni, m)
y <- matrix(rnorm(ni*q), ni, q)
cancor(w,y)$cor[1]
rcancor(w,y)$out$cor[1]

ni <- 1500 
w <- matrix(rnorm(ni*m), ni, m)
y <- matrix(rnorm(ni*q), ni, q)
cancor(w,y)$cor[1]
rcancor(w,y)$out$cor[1]

ni <- 2000 
w <- matrix(rnorm(ni*m), ni, m)
y <- matrix(rnorm(ni*q), ni, q)
cancor(w,y)$cor[1]
rcancor(w,y)$out$cor[1]

ni <- 2500 
w <- matrix(rnorm(ni*m), ni, m)
y <- matrix(rnorm(ni*q), ni, q)
cancor(w,y)$cor[1]
rcancor(w,y)$out$cor[1]

ni <- 3000 
w <- matrix(rnorm(ni*m), ni, m)
y <- matrix(rnorm(ni*q), ni, q)
cancor(w,y)$cor[1]
rcancor(w,y)$out$cor[1]

ni <- 3500 
w <- matrix(rnorm(ni*m), ni, m)
y <- matrix(rnorm(ni*q), ni, q)
cancor(w,y)$cor[1]
rcancor(w,y)$out$cor[1]

ni <- 4000 
w <- matrix(rnorm(ni*m), ni, m)
y <- matrix(rnorm(ni*q), ni, q)
cancor(w,y)$cor[1]
rcancor(w,y)$out$cor[1]

ni <- 4500 
w <- matrix(rnorm(ni*m), ni, m)
y <- matrix(rnorm(ni*q), ni, q)
cancor(w,y)$cor[1]
rcancor(w,y)$out$cor[1]

ni <- 5000 
w <- matrix(rnorm(ni*m), ni, m)
y <- matrix(rnorm(ni*q), ni, q)
cancor(w,y)$cor[1]
rcancor(w,y)$out$cor[1]


##problem 8.6, needs mrobdata

#a)

group <- pottery[pottery[,1]!=5,1]
group <- (as.integer(group!=1)) + 1
x <- pottery[pottery[,1]!=5,-1]

#b)

z <- x[,1:9]
library(MASS)
out <- qda(x=z,group)
1-mean(predict(out,z)$class==group)

#c)

#delete the first variable
out <- qda(x=z[,-c(1)],group)
1-mean(predict(out,z[,-c(1)])$class==group)

#delete the first two variables
out <- qda(x=z[,-c(1,2)],group)
1-mean(predict(out,z[,-c(1,2)])$class==group)

#deleting the second variable caused a classification error
#so delete the first and third variable
out <- qda(x=z[,-c(1,3)],group)
1-mean(predict(out,z[,-c(1,3)])$class==group)

#deleting the third variable caused a classification error
#so delete the first and fourth variable
out <- qda(x=z[,-c(1,4)],group)
1-mean(predict(out,z[,-c(1,4)])$class==group)

#do not need the 4th variable so continue deleting variables
out <- qda(x=z[,-c(1,4,5)],group)
1-mean(predict(out,z[,-c(1,4,5)])$class==group)


##Problem 8.7
#Get x as in Problem 8.6 a).

out<-lda(x[,c(6,11,14,18)],group)
1-mean(predict(out,x[,c(6,11,14,18)])$class==group)


##Problem 8.8, needs mpack

#a)

n <- 100
p <- 2
x <- matrix(rnorm(n*p),nrow=n,ncol=p)
group <- 1 + 0*1:n
covv <- diag(p)
mns<- apply(x, 2, mean)
md2 <- mahalanobis(x, center = mns, covv)
group[md2>qchisq(0.5,p)] <- 2

#b) 

out1 <- ddiscr(x,w=x,group,xwflag=T)
out1$err
out1$toterr

#c)

out2<-ddiscr(x=out1$mdx[,1],w=out1$mdw[,1],group,xwflag=T)
out2$err
out2$toterr

#d)

out3 <- ddiscr2(x,w=x,group,xwflag=T)
out3$err
out3$toterr

#e)

out4<-ddiscr2(x=out1$mdx[,1],w=out1$mdw[,1],group,xwflag=T)
out4$err
out4$toterr

#f)

library(MASS)

out5  <- lda(x,group)
1-mean(predict(out5,x)$class==group)

#g)

out6  <- qda(x,group)
1-mean(predict(out6,x)$class==group)


##Problem 9.4, needs mpack

n<-4000; p<-30
zout <- rhotsim(n=4000,p=30)
SRHOT <- zout$rhot/(1.04 + 0.12/p + (40+p)/n)
HOT <- zout$hot
plot(SRHOT,HOT)
abline(0,1)


##Problem 9.5, needs mpack


source("G:/mpack.txt") 
source("G:/mrobdata.txt")


#a)

rhotsim(n=150,p=10,runs=5000,delta=0)

#b) 

rhotsim(n=150,p=10,runs=5000,delta=0.2)


##Problem 10.3, needs mpack, mrobdata


#a)

out<- manova1w(y=turtle[,1:3],p=2,group=turtle[,4]) #right click Stop 6 times

#b)

ddplot4(out$res) #right click Stop on the plot once

##Problem 10.4

tear <- c(6.5, 6.2, 5.8, 6.5, 6.5, 6.9, 7.2, 6.9, 6.1, 
6.3, 6.7, 6.6, 7.2, 7.1, 6.8, 7.1, 7.0, 7.2, 7.5, 7.6)
gloss <- c(9.5, 9.9, 9.6, 9.6, 9.2, 9.1, 10.0, 9.9, 
9.5, 9.4, 9.1, 9.3, 8.3, 8.4, 8.5, 9.2, 8.8, 9.7, 
10.1, 9.2)
opacity <- c(4.4, 6.4, 3.0, 4.1, 0.8, 5.7, 2.0, 3.9, 
1.9, 5.7, 2.8, 4.1, 3.8, 1.6, 3.4, 8.4, 5.2, 6.9, 
2.7, 1.9)

Y <- cbind(tear, gloss, opacity)

rate <- factor(gl(2,10), labels=c("Low", "High"))

additive <- factor(gl(2, 5, length=20), labels=c("Low", "High"))

#one way MANOVA model

fit <- manova(Y ~ rate )

summary.aov(fit)           #univariate ANOVA tables

summary(fit)               #MANOVA table with Pillai

summary(fit, test="Wilks") #Wilks' lambda

summary(fit,test = "Hotelling-Lawley")

summary(fit,test = "Roy")

#for one way MANOVA with df = 1 (p = 2 groups), 
#the 4 tests are the same = 
#two sample Hotelling's T^2 test
grp <- as.integer(rate)
out<-manova1w(y=Y,p=2,group=grp)  #right click Stop 6 times
summary(out$out,test="Hotelling-Lawley")

#for two way MANOVA, the 4 tests seem to be the same
#for df = 1 or p = 2 groups
#two way MANOVA model
fit <- manova(Y ~ rate * additive)

summary.aov(fit)   # univariate two way ANOVA tables

summary(fit, test="Wilks") # Wilks' lambda

summary(fit)            # Pillai's test: the default

                                      
#delete the interaction to get the additive model
fit <- manova(Y ~ rate + additive)

summary.aov(fit)   # univariate two way ANOVA tables

summary(fit, test="Wilks") #MANOVA withWilks' lambda

summary(fit)               #Pillai
          


##problem 11.4, needs mpack, mrobdata

#a) 

library(MASS)

z <- cbind(buxx,buxy)
covhat <- var(z)
out1 <- factanal(factors = 2, covmat=covhat,rotation="promax")
out1

#b) 

rcovhat <- covrmvn(z)$cov
out2 <- factanal(factors = 2, covmat=rcovhat,rotation="promax")
out2


##Problem 12.8, needs mpack, mrobdata

#a)

y <- log(mussels)[,4:5]
x <- mussels[,1:3]
x[,2] <- log(x[,2])
z<-cbind(x,y)
out <- mltreg(x,y,indices=c(3,4)) #right click Stop 4 times

#b)

out

#c)

ddplot4(out$res)  #right click Stop once


##Problem 12.9, needs mpack and mrobdata

#a) 

x <- fitness[,1:3]
y <- fitness[,4:6]
out <- mltreg(x,y,indices=c(2,4)) #right click Stop on the plot 6 times


#b)

ddplot4(out$res) #right click Stop once


##Problem 12.10

#a)

mregsim(n=20,m=2,p=4,nruns=5000,etype=1)


#b)

mregsim(n=20,m=2,p=4,nruns=5000,etype=1,mnull=F)


##Problem 12.11

mpredsim(n=100,m=2,p=4,nruns=5000,mnull=F,alpha=0.1)


##Ch 12 problem

y <- log(mussels)[,4:5]
x <- mussels[,1:3]
x[,2] <- log(x[,2])
z<-cbind(x,y)
pairs(z, labels=c("L","log(W)","H","log(S)","log(M)"))

ddplot4(z)

out <- mltreg(x,y)

#for Y2, cases 8, 25 and 48 are not fit well
ddplot4(out$res)


out <- rmltreg(x,y)

ddplot4(out$res)


##ch 12 problem
##outliers

x <- fitness[,1:3]
y <- fitness[,4:6]

z<-cbind(x,y)
pairs(z, labels=c("Wt","Waist","Pulse","ChUps","Situps","Jumps"))

ddplot4(z)

out <- mltreg(x,y)

ddplot4(out$res)

out <- rmltreg(x,y)

ddplot4(out$res)


##Problem 13.1

x <- cbind(buxx,buxy)
out<-kmeans(x,2,nstart=25)
plot(x, col = out$cluster)
points(out$centers, col = 1:2, pch = 8, cex=2)

##Problem 13.2

x <- cbind(buxx,buxy)
out <- hclust(dist(x),"complete") 
#complete is the default
plot(out)
#plot(out,hang=-1)


##Problem 13.3

x1 <- c(-0.6,  0.1,-1.5,-1.4, 1.1,-0.9, 1.4, 0.6,0  )
x2 <- c(-1  ,-0.75,-0.4,-1.6,-0.3,-1.2,   0,-0.2,0.7)
x <- cbind(x1,x2)
##out<-hclust(x) #errors
out <- hclust(dist(x))
plot(out)
plot(x[,1],x[,2])
library(cluster)
out<-agnes(x)
plot(out) #right click twice


### Problem 13.4

#a)

mbb1415 <- read.csv(file.choose(), header=T)
#select mbb1415 from F drive, or where you stored the file

#b) 

out <- hclust(dist(mbb1415[,c(7,13)]))
plot(out) #uses index as label
plot(out, labels=as.character(mbb1415[,1]))
#uses position as label

#c) 

quantile(mbb1415[,3])
mcut <- mbb1415[mbb1415[,3]>182,]
out <- hclust(dist(mcut[,c(7,13)]))
plot(out, labels=as.character(mcut[,1]))


### Problem 13.5

#a)

wbb1415 <- read.csv(file.choose(), header=T)
#select wbb1415 from F drive, or where you stored the file

#b) 

out <- hclust(dist(wbb1415[,c(7,13)]))
plot(out) #uses index as label
plot(out, labels=as.character(wbb1415[,1]))
#uses position as label

#c) 

quantile(wbb1415[,3])
wcut <- wbb1415[wbb1415[,3]>182,]
out <- hclust(dist(wcut[,c(7,13)]))
plot(out, labels=as.character(wcut[,1]))


##Problem 14.4

zgam <- c(.01,.05,.1,.15,.2,.25,.3,.35,.4,.45,.5)
pifclean(3000,zgam)


##Problem 14.5

zh <- c(10,20,30,40,50,60,70,80,90,100)
for(i in 1:10) gamper(zh[i])


#Problem 14.6 R commands, needs mpack

#b)

nltv(0.5)
nltv(0.75) 
nltv(0.9)
nltv(0.9999)

#Problem 14.7 R commands, needs mpack

#b)

deltv(0.5)
deltv(0.75) 
deltv(0.9)
deltv(0.9999)

#Problem 14.8 R commands, needs mpack

#b)

cltv(0.5)
cltv(0.75) 
cltv(0.9)
cltv(0.9999)

#Problem 14.9 R commands, needs mpack and mrobdata

#b)

mbamv(belx,bely) #right click Stop seven times

#c)

#for better looking plots than the default

cmar <- par("mar")
par(mfrow = c(2, 1))
par(mar=c(4.0,4.0,2.0,0.5))
#main command is below

mbamv2(buxx,buxy) #right click Stop fourteen times


#restore default

par(mfrow=c(1,1))
par(mar=cmar)

#Problem 14.10 R commands, needs mpack and mrobdata

#b)

mlrplot2(belx,bely) #right click Stop four times

#c)

mlrplot2(cbrainx,cbrainy) #right click Stop four times

#d)

mlrplot2(museum[,3:11],museum[,2]) #right click Stop four times

#e)

mlrplot2(buxx,buxy) #right click Stop four times

#f)

#use the following command several times
mlrplot2(hx,hy) #right click Stop four times

library(MASS) 
ffplot2(hx,hy)


##Problem 14.11, needs mpack

n3x <- matrix(rnorm(300),nrow=100,ncol=3)
ln3x <- exp(n3x)
library(MASS)

#a)

pairs(n3x)
pairs(ln3x)

#b)

ncy <- (n3x%*%1:3)^3 + 0.1*rnorm(100)
plot(n3x%*%(1:3),ncy)

#c)

trviews(n3x, ncy) #right click Stop 10 times

#d)

#change betahat coefficients below

plot(n3x%*%c(22.469,61.242,75.285),n3x%*%1:3)

#e)

lncy <- (ln3x%*%1:3)^3 + 0.1*rnorm(100)
trviews(ln3x,lncy) #right click Stop 10 times

#f)

#change betahat coefficients below
plot(ln3x%*%c(94.848,216.719,328.444),ln3x%*%1:3)


##Problem 14.12, needs mpack

#a)

library(MASS)
nx <- matrix(rnorm(300),nrow=100,ncol=3)
lnx <- exp(nx)
SP <- lnx%*%1:3
lnsincy <- sin(SP)/SP + 0.01*rnorm(100)

#b)

trviews(lnx,lnsincy)     #right click Stop 10 times
essp(lnx,lnsincy,M=40)   #right click Stop

#c)

ctrviews(lnx,lnsincy)    #right click Stop 10 times

#d)

lmsviews(lnx,lnsincy)     #right click Stop 10 times


##Problem 14.13, needs mpack

#a)

nx <- matrix(rnorm(300),nrow=100,ncol=3)
SP <- nx%*%1:3
ncuby <- SP^3 + rnorm(100)
nexpy <- exp(SP) + rnorm(100)
nlinsy <- SP + 4*sin(SP) + 0.1*rnorm(100)
nsincy <- sin(SP)/SP + 0.01*rnorm(100)
nsiny <- sin(SP) + 0.1*rnorm(100)
nsqrty <- sqrt(abs(SP)) + 0.1*rnorm(100)
nsqy <- SP^2 + rnorm(100)

#b)

plot(SP,ncuby)
plot(-SP,ncuby)

#c)

library(MASS)
trviews(nx,ncuby)     #right click Stop 10 times
essp(nx,ncuby, M=40)  #right click Stop

#d)

tem <- c(12.60514, 25.06613, 37.25504)
ESP <- nx%*%tem
plot(ESP,SP)

#e)

lnx <- exp(nx)
SP <- lnx%*%1:3
lncuby <- (SP/3)^3 + rnorm(100)
lnlinsy <- SP + 10*sin(SP) + 0.1*rnorm(100)
lnsincy <- sin(SP)/SP + 0.01*rnorm(100)
lnsiny <- sin(SP/3) + 0.1*rnorm(100)
ESP <- lnx%*%tem


##Problem 14.15, needs mpack and mrobdata

library(MASS)
tvreg(buxx,buxy,ii=1)     #right click Stop 10 times
tvreg2(buxx,buxy, M = 0)  #right click Stop


##Problem 14.16
#a)

cmar <- par("mar")
par(mfrow = c(2, 1))
par(mar=c(4.0,4.0,2.0,0.5))
#main command is below

MLRplot(buxx,buxy)   #right click Stop twice

#restore default

par(mfrow=c(1,1))
par(mar=cmar)

#c)

library(MASS)
rrplot2(buxx,buxy)


#d)

ffplot2(buxx,buxy)



###############SAS problem below


****HW5 problem F) SAS commands ;
*SAS Institute (1985, p. 626-630);
options ls = 70;
data Temperature;
      title 'Mean Temperature in January and July for Selected Cities ';
      input City $1-15 January July;
      datalines;
Mobile          51.2 81.6
Phoenix         51.2 91.2
Little Rock     39.5 81.4
Sacramento      45.1 75.2
Denver          29.9 73.0
Hartford        24.8 72.7
Wilmington      32.0 75.8
Washington DC   35.6 78.7
Jacksonville    54.6 81.0
Miami           67.2 82.3
Atlanta         42.4 78.0
Boise           29.0 74.5
Chicago         22.9 71.9
Peoria          23.8 75.1
Indianapolis    27.9 75.0
Des Moines      19.4 75.1
Wichita         31.3 80.7
Louisville      33.3 76.9
New Orleans     52.9 81.9
Portland, ME    21.5 68.0
Baltimore       33.4 76.6
Boston          29.2 73.3
Detroit         25.5 73.3
Sault Ste Marie 14.2 63.8
Duluth           8.5 65.6
Minneapolis     12.2 71.9
Jackson         47.1 81.7
Kansas City     27.8 78.8
St Louis        31.3 78.6
Great Falls     20.5 69.3
Omaha           22.6 77.2
Reno            31.9 69.3
Concord         20.6 69.7
Atlantic City   32.7 75.1
Albuquerque     35.2 78.7
Albany          21.5 72.0
Buffalo         23.7 70.1
New York        32.2 76.6
Charlotte       42.1 78.5
Raleigh         40.5 77.5
Bismarck         8.2 70.8
Cincinnati      31.1 75.6
Cleveland       26.9 71.4
Columbus        28.4 73.6
Oklahoma City   36.8 81.5
Portland, OR    38.1 67.1
Philadelphia    32.3 76.8
Pittsburgh      28.1 71.9
Providence      28.4 72.1
Columbia        45.4 81.2
Sioux Falls     14.2 73.3
Memphis         40.5 79.6
Nashville       38.3 79.6
Dallas          44.8 84.8
El Paso         43.6 82.3
Houston         52.1 83.3
Salt Lake City  28.0 76.7
Burlington      16.8 69.8
Norfolk         40.5 78.3
Richmond        37.5 77.9
Spokane         25.4 69.7
Charleston, WV  34.5 75.0
Milwaukee       19.4 69.9
Cheyenne        26.6 69.1
   ;

 title 'Mean Temperature in January and July for Selected Cities';
   title2 'Plot of Raw Data';
   proc plot;
      plot July*January=City/vpos=36;
   run;

 
   title 'Mean Temperature in January and July for Selected Cities';
*leave cov out to use correlation matrix;
   proc princomp data=Temperature cov out=prin;
      var July January;
   run;
   proc plot;
     plot prin2*prin1/vpos=26;
   run;