Lesson 3: R graphics

This lecture covers the following section in An Introduction to R:

Graphics input and output

One of the nice features of R is its graphics subroutines. Let us starts with the time series econ.dat

> X11()                            # To create a graphics window
> x<-scan(file="econ.dat",skip=2) 
> ts.plot(x)

econ.ps

To get hard copies, use (pp. 84-85)
> postscript("econ.ps", horizontal=F,height=6.width=4)
> ts.plot(x)
> dev.off()
The graph is in the file "econ.ps".

Histogram, CDF and density function

Histogram (related to density) and sample distribution function can be plotted by (pp. 41-42).

> par(mfrow=c(2,2))
> hist(x)                         # default output
> hist(x,nclass=20)               # with 20 interval
> hist(x,nclass=20,probability=T) # with density label in the y-axis
> library(stepfun)
> plot(ecdf(x))                   # CDF 

histogram.ps

Note that the par(mfrow=c(2,2)) instruction is the partition one page into a (2,2) compartments. It can be (m,n) for any m and n. Here are plots of some of the density functions (p.40, the preface d is for density)

postscript("density.ps")
par(mfrow=c(2,2))
x<-seq(-4,4,length=200)
y<-dnorm(x)
plot(x,y,type="l");title("Normal density")
x<-seq(0,10,length=200)
y<-dexp(x)
plot(x,y,type="l");title("Exponential density")
x<-seq(0,10,length=200)
y<-dgamma(x,2)
plot(x,y,type="l");title("Gamma(2) density")
x<-seq(0,10,length=200)
y<-dlnorm(x)
plot(x,y,type="l");title("Log-normal density")
dev.off()

density.ps

Scatter plot and curve fitting

The following data measures the realtion between temperqture the viscosity of two lubricants: Raji and Molt. It is saved as patel.dat.
Dr. Rajiv Patel's Data,  
 T  Raji Molt
40 0.840 1.150
38 0.930 1.150
36 1.050 1.210
34 1.020 1.130
32 1.050 1.180
30 1.150 1.290
28 1.140 1.350
26 1.250 1.430
24 1.290 1.530
22 1.410 1.660
20 1.530 1.720
18 1.590 1.770
16 1.680 1.930
14 1.820 2.070
12 1.890 2.180
10 2.020 2.240
 8 2.020 2.380
 6 2.220 2.510
 4 2.290 2.670
The following program is used to view the data.
 
postscript("patel.ps",height=5,width=6)
ex1<-scan("patel.dat",list(0,0,0),skip=2) 
TC<-ex1[[1]]
RAJ<-ex1[[2]]
MOLT<-ex1[[3]]

x<-c(2:42)    # producing theoretically fitted curves
cv1<- 3.0105*10**19/(273+x)**7.82529
cv2<- 3.46524*10**19/(273+x)**7.82529

x0<-5
y01<-0.8
y02<-0.5
x1<-c(8,14)
y1<-c(0.8,0.8)
y2<-c(0.5,0.5)

plot(TC,RAJ,xlim=range(0,50),ylim=range(0,3.5),
xlab="temperature (degree C)", 
ylab="apparent microviscosity",
type="p",pch=0)                     # pch =0  means empty square
                                    # Add more points (overlay)
points(TC,MOLT,type="p",pch=18)     # pch =18 means black diamond

title("Raji and Molt microviscosity plot")

lines(x,cv1)                    # Add a straight line
lines(x,cv2,lty=4)              # add another line

points(x0,y01,type="p",pch=0)   # add a point at (x0,y01)
points(x0,y02,type="p",pch=18)

lines(x1,y1)
lines(x1,y2,lty=4)

text(19,0.8,"Raji")            # text "Raji" at (19, 0.8)
text(19,0.5,"Molt")
dev.off()

patel.ps

If you wish to fit with a straight (linear regression) line. The computer will find the regression line for you.

postscript("linearfit.ps",height=4,width=5)
plot(MOLT,RAJ,xlim=range(0,3),ylim=range(0,3),
xlab="MOLT viscosity", ylab="RAJ viscosity",type="p",pch=3) # pch 3 = +

abline(lsfit(MOLT,RAJ),lty=2) # lsfit=least square fit
                              # lty (line type) = = => dashed line
dev.off()

linearfit.ps

Longitudinal data display (high-level overlay)

The following data were the weights over time (weeks) of rats under the treatemnt of thiourocil versus control. Both groups had 10 individuals. It is saved as longit.dat.
 
     1  2  3   4   5
CON 57 86 114 139 172
CON 60 93 123 146 177
CON 52 77 111 144 185
CON 49 67 100 129 164
 ..................
THI 56 78  95 103 108
THI 58 69  93 116 140
THI 46 61  78  90 107
THI 53 72  89 104 122
To disply them, we may use the following.
postscript("longit.ps",height=6,width=7.5)
x<-scan("longit.dat",list("",0,0,0,0,0),skip=1)
ym<-matrix(0,20,5)
for (i in c(1:5)) ym[,i]<-x[[i+1]]
ym  # To make sure that this ym recovers the original matrix
xaxis<-c(1:5)
plot(xaxis,ym[1,],xlim=range(0,6),ylim=range(50,180), xlab="time in weeks",
     ylab="Weight",type="both",lty=1,col="blue")

  # type="both" means both points "p" and line "l"

for (i in c(2:10)) points(xaxis,ym[i, ],type="both",lty=1,col="blue")
for (i in c(11:20)) points(xaxis,ym[i, ],type="both",lty=1,col="red")
title("Thyroxin (red) vs control (blue)")
dev.off()

longit.ps

Multivariate data display

If we want to know the time relation between the weights, we can use pairs() (p. 73). (Apparently, the relation losened when times were far away.)
postscript("pairs.ps")
pairs(ym)
dev.off()

pairs.ps

Three dimensional graphics

Here is a 3-D graph using persp() (p.74).
postscript("twoballs.ps",height=6,width=6)
x<-seq(0,1,length=100) 
y<-seq(0,1,length=100)
z<-matrix(0,100,100)
for(i in 1:100){for(j in 1:100){if ((x[i]-.5)^2+(y[j]-0.5)^2<0.09) z[i,j]<-(0.09-(x[i]-0.5)^2-(y[j]-0.5)^2)^0.5}}
for(i in 1:100){for(j in 1:100){if ((x[i]-.75)^2+(y[j]-0.25)^2<0.04) z[i,j]<-(0.04-(x[i]-0.75)^2-(y[j]-0.25)^2)^0.5}}
persp(x=seq(0,1,len=100),y=seq(0,1,len=100),z,xlim=range(0,1),
ylim=range(0,1),zlim=range(0,1.0))
dev.off()

twoballs.ps

Return to R home page