Demonstrations for teaching and learning

Demonstrations for teaching and learning [R logo]

Help pages for package `TeachingDemos' version 1.4

bct | Box-Cox Transforms
char2seed | Convert a character string into a random seed
chisq.detail | Print details of a chi-squared test
ci.examp | Plot examples of Confidence Intervals
clipplot | Clip plotting to a rectangular region
clt.examp | Plot Examples of the Central Limit Theorem
cnvrt.coords | Convert between the 5 different coordinate sytems on a graphical device
coin.faces | Designs for coin faces for use with plot.rgl.coin
dice | Simulate rolling dice
dots | Create a quick dotchart (histogram)
dots2 | Create a quick dotchart (histogram)
emp.hpd | Compute Highest Posterior Density Intervals
face2.plot | Chernoff Faces
faces | Chernoff Faces
faces2 | Chernoff Faces
fagan.plot | Create a Fagan plot to demonstrate Bayes Theorem and screening tests
flip.rgl.coin | Animated die roll or coin flip
gp.close | Alpha version functions to send plotting commands to GnuPlot
gp.open | Alpha version functions to send plotting commands to GnuPlot
gp.plot | Alpha version functions to send plotting commands to GnuPlot
gp.send | Alpha version functions to send plotting commands to GnuPlot
hpd | Compute Highest Posterior Density Intervals
identify.Map | Identify polygons within a plotted map from the maptools packages.
identify.polylist | Identify polygons within a plotted map from the maptools packages.
lattice.demo | Interactively explore the conditioned panels in lattice plots.
limits.g | Functions to plot Shewhart g chart
loess.demo | Demonstrate the internals of loess curve fits
mle.demo | Demonstrate the basic concept of Maximum Likelihood Estimation
panel.dice | Simulate rolling dice
plot.dice | Simulate rolling dice
plot.rgl.coin | Animated die roll or coin flip
plot.rgl.die | Animated die roll or coin flip
plotFagan | Create a Fagan plot to demonstrate Bayes Theorem and screening tests
plotFagan2 | Create a Fagan plot to demonstrate Bayes Theorem and screening tests
power.examp | Graphically illustrate the concept of power.
power.refresh | Graphically illustrate the concept of power.
prepanel.dice | Simulate rolling dice
project.Map | Apply projections to Map objects (from maptools package)
put.points.demo | Demonstrate Correlation and Regression by placing and moving data points
recenter.Map | Recenter a Map object from the maptools package.
rgl.Map | Plot a map in an rgl window
roc.demo | Demonstrate ROC curves by interactively building one
roll.rgl.die | Animated die roll or coin flip
rotate.cloud | Interactively rotate 3D plots
rotate.persp | Interactively rotate 3D plots
rotate.wireframe | Interactively rotate 3D plots
run.ci.examp | Plot examples of Confidence Intervals
run.cor.examp | Interactively demonstrate correlations
run.cor2.examp | Interactively demonstrate correlations
run.hist.demo | Create a histogram and interactively change the number of bars.
run.power.examp | Graphically illustrate the concept of power.
sd.g | Functions to plot Shewhart g chart
slider | slider / button control widgets
sliderv | Create a Tk slider window
squishplot | Squish the plotting area to a specified aspect ratio
state.vbm | Map object for the US State Visibility Based Map
stats.g | Functions to plot Shewhart g chart
strip.shingle | strip function for trellis plots
subplot | Embed a new plot within an existing plot
tkBrush | Change the Color and Styles of points interactively
tree.demo | Interactively demonstrate regression trees
triplot | Create or add to a Trilinear Plot
vis.binom | Plot various distributions then interactivly adjust the parameters.
vis.boxcox | Interactively visualize Box-Cox transformations
vis.boxcoxu | Interactively visualize Box-Cox transformations
vis.gamma | Plot various distributions then interactivly adjust the parameters.
vis.normal | Plot various distributions then interactivly adjust the parameters.
vis.t | Plot various distributions then interactivly adjust the parameters.
z.test | Z test for known population standard deviation


Box-Cox Transforms

Description

Computes the Box-Cox transform of the data for a given value of lambda. Includes the scaling factor.

Usage

bct(y, lambda)

Arguments

y | Vector of data to be transformed.
lambda | Scalar exponent for transform (1 is linear, 0 is log).

Details

bct computes the Box-Cox family of transforms: y = (y^lambda - 1)/(lambda*gm^(lambda-1)), where gm is the geometric mean of the y's. returns log(y)*gm when lambda equals 0.

Value

A vector of the same length as y with the corresponding transformed values.

Author(s)

Greg Snow greg.snow@intermountainmail.org

See Also

vis.boxcox, vis.boxcoxu, boxcox in package MASS, other implementations in various packages

Examples

y <- rlnorm(500, 3, 2)
par(mfrow=c(2,2))
qqnorm(y)
qqnorm(bct(y,1/2))
qqnorm(bct(y,0))
hist(bct(y,0))

## The function is currently defined as
function(y,lambda){

  gm <- exp( mean( log(y) ) )

  if(lambda==0) return( log(y)*gm )

  yt <- (y^lambda - 1)/( lambda * gm^(lambda-1) )
  return(yt)
  }
[Package TeachingDemos version 1.4 Index]


Convert a character string into a random seed

Description

This function creates a seed for the random number generator from a character string. Character strings can be based on student names so that every student has a different random sample, but the teacher can generate the same datasets.

Usage

char2seed(x, set = TRUE, ...)

Arguments

x | A character string
set | Logical, should the seed be set or just returned
... | Additional parameters passed on to set.seed

Details

Simulations or other situations call for the need to have repeatable random numbers, it is easier to remember a word or string than a number, so this function converts words or character strings to an integer and optionally sets the seed based on this.

Teachers can assign students to generate a random dataset using their name to seed the rng, this way each student will have a different dataset, but the teacher can generate the same set of data to check values.

Any characters other than letters (a-zA-Z) or digits (0-9) will be silently removed. This function is not case sensitive, so "ABC" and "abc" will generate the same seed.

This is a many to one function, so it is possible to find different words that generate the same seed, but this is unlikely by chance alone.

Value

This returns an integer (but mode numeric) to use as a seed for the RNG. If set is true then it is returned invisibly.

Author(s)

Greg Snow greg.snow@intermountainmail.org

See Also

set.seed

Examples

char2seed('Snow')
x <- rnorm(100)
rnorm(10)
tmp <- char2seed('Snow',set=FALSE)
set.seed(tmp)
y <- rnorm(100)

all.equal(x,y) # should be true


Print details of a chi-squared test

Description

Prints out the details of the computations involved in a chi-squared test on a table. Includes the expected values and the chi-squared contribution of each cell.

Usage

chisq.detail(tab)

Arguments

tab Matrix or table to be analyzed

Details

This function prints out the input table along with the expected value for each cell under the null hypothesis. It also prints out the chi-squared contribution of each cell in the same pattern as the table. This shows the computations involved and one rule of thumb is to look for these values that are greater than 4 as a post-hoc analysis.

Value

This function is used primarily for its side effect of printing the results, but does return invisibly a list with the following components:

obs | A matrix of the observed values, same as tab.
expected | A matrix of the expected values under the null hypothesis.
chi.table | A matrix of the chi-squared contributions of each cell.
chi2 | The chi-squared test statistic.

Author(s)

Greg Snow, greg.snow@intermountainmail.org

References

~put references to the literature/web site here ~ Moore, bps

See Also

chisq.test,loglin, xtabs, table, prop.table, CrossTable from the gmodels package.

Examples

chisq.detail(HairEyeColor[,,1])
chisq.detail(HairEyeColor[,,2])


Plot examples of Confidence Intervals

Description

Generate reps samples from a normal distribution then compute and plot confidence intervals for each sample along with information about the population to demonstrate confidence intervals. Optionally change the confidence level using a Tk slider.

Usage

ci.examp(mean.sim = 100, sd = 10, n = 25, reps = 50, conf.level = 0.95, method = "z", lower.conf = (1 - conf.level)/2, upper.conf = 1 - (1 - conf.level)/2)
run.ci.examp(reps = 100, seed, method="z", n=25)

Arguments

mean.sim | The mean of the population.
sd | The standard deviation of the population.
n | The sample size for each sample.
reps | The number of samples/intervals to create.
conf.level | The confidence level of the intervals.
method | 'z', 't', or 'both', should the intervals be based on the normal, the t, or both distributions.
lower.conf | Quantile for lower confidence bound.
upper.conf | Quantile for upper confidence bound.
seed | The seed to use for the random number generation.

Details

These functions demonstrate the concept of confidence intervals by taking multiple samples from a known normal distribution and calculating a confidence interval for each sample and plotting the interval relative to the true mean. Intervals that contain the true mean will be plotted in black and those that do not include the true mean will be plotted in different colors.

The method argument determines the type of interval: 'z' will use the normal distribution and the known population standard deviation, 't' will use the t distribution and the sample standard deviations, 'both' will compute both for each sample for easy comparison (it is best to reduce reps to about 25 when using 'both').

The optional arguments lower.conf and upper.conf can be used to plot non-symmetric or 1 sided confidence intervals.

The function run.ci.examp also creates a Tk slider that will allow you to interactively change the confidence level and replot the intervals to show how the interval widths change with the confidence level.

Value

These functions are run solely for the side effect of plotting the intervals, there is no meaningfull return value.

Author(s)

Greg Snow greg.snow@intermountainmail.org

See Also

z.test, t.test

Examples

ci.examp()
## Not run: run.ci.examp()

# 1 sided confidence intervals
ci.examp(lower.conf=0, upper.conf=0.95)

# non-symmetric intervals
ci.examp(lower.conf=0.02, upper.conf=0.97)


Clip plotting to a rectangular region

Description

Clip plotting to a rectangular region that is a subset of the plotting area

Usage

clipplot(fun, xlim = par("usr")[1:2], ylim = par("usr")[3:4])

Arguments

fun | The function or expression to do the plotting.
xlim | A vector of length 2 representing the x-limits to clip plotting to, defaults to the entire width of the plotting region.
ylim | A vector of length 2 representing the y-limits to clip the plot to, defaults to the entire height of the plotting region.

Details

This function resets the active region for plotting to a rectangle within the plotting area and turns on clipping so that any points, lines, etc. that are outside the rectange are not plotted.

A side effect of this function is a call to the box() command, most of the time this will have no visible effect unless you suppressed the drawing of the box initially.

Value

Nothing meaningful is returned

Author(s)

Greg Snow greg.snow@intermountainmail.org

See Also

par, lines

Examples

x <- seq(1,100)
y <- rnorm(100)
plot(x,y, type='b', col='blue')
clipplot( lines(x,y, type='b', col='red'), ylim=c(par('usr')[3],0))

attach(iris)

tmp <- c('red','green','blue')
names(tmp) <- levels(Species)
plot(Petal.Width,Petal.Length, col=tmp[Species])
for(s in levels(Species)){
  clipplot( abline(
    lm(Petal.Length~Petal.Width, data=iris, subset=Species==s),
    col=tmp[s]),
    xlim=range(Petal.Width[Species==s]))
}

detach(iris)


Plot Examples of the Central Limit Theorem

Description

Takes samples of size n from 4 different distributions and plots histograms of the means along with a normal curve with matching mean and standard deviation. Creating the plots for different values of n demonstrates the Central Limit Theorem.

Usage

clt.examp(n = 1, reps = 10000, nclass = 16)

Arguments

n | size of the individual samples
reps | number of samples to take from each distribution
nclass | number of bars in the histograms

Details

The 4 distributions sampled from are a Normal with mean 0 and standard deviation 1, an exponential with lambda 1/3 (mean = 3), a uniform distribution from 0 to 1, and a beta distribution with alpha 0.35 and beta 0.25 (U shaped left skewed).

Running the function with n=1 will show the populations. Run the function again with n at higher values to show that the sampling distribution of the uniform quickly becomes normal and the exponential and beta distributions eventually become normal (but much slower than the uniform).

Value

This function is run for its side effect of creating plots. It returns NULL invisibly.

Author(s)

Greg Snow greg.snow@intermountainmail.org

See Also

rnorm, rexp, runif, rbeta

Examples

clt.examp()
clt.examp(5)
clt.examp(30)
clt.examp(50)


Convert between the 5 different coordinate sytems on a graphical device

Description

Takes a set of coordinates in any of the 5 coordinate systems (usr, plt, fig, dev, or tdev) and returns the same points in all 5 coordinate systems.

Usage

cnvrt.coords(x, y = NULL, input = c("usr", "plt", "fig", "dev","tdev"))

Arguments

x | Vector, Matrix, or list of x coordinates (or x and y coordinates), NA's allowed.
y | y coordinates (if x is a vector), NA's allowed.
input | Character scalar indicating the coordinate system of the input points.

Details

Every plot has 5 coordinate systems:

usr (User): the coordinate system of the data, this is shown by the tick marks and axis labels.

plt (Plot): Plot area, coordinates range from 0 to 1 with 0 corresponding to the x and y axes and 1 corresponding to the top and right of the plot area. Margins of the plot correspond to plot coordinates less than 0 or greater than 1.

fig (Figure): Figure area, coordinates range from 0 to 1 with 0 corresponding to the bottom and left edges of the figure (including margins, label areas) and 1 corresponds to the top and right edges. fig and dev coordinates will be identical if there is only 1 figure area on the device (layout, mfrow, or mfcol has not been used).

dev (Device): Device area, coordinates range from 0 to 1 with 0 corresponding to the bottom and left of the device region within the outer margins and 1 is the top and right of the region withing the outer margins. If the outer margins are all set to 0 then tdev and dev should be identical.

tdev (Total Device): Total Device area, coordinates range from 0 to 1 with 0 corresponding to the bottom and left edges of the device (piece of paper, window on screen) and 1 corresponds to the top and right edges.

Value

A list with 5 components, each component is a list with vectors named x and y. The 5 sublists are:

usr | The coordinates of the input points in usr (User) coordinates.
plt | The coordinates of the input points in plt (Plot) coordinates.
fig | The coordinates of the input points in fig (Figure) coordinates.
dev | The coordinates of the input points in dev (Device) coordinates.
tdev | The coordinates of the input points in tdev (Total Device) coordinates.

Note

You must provide both x and y, but one of them may be NA.

Author(s)

Greg Snow greg.snow@intermountainmail.org

See Also

par specifically 'usr','plt', and 'fig'. Also 'xpd' for plotting outside of the plotting region and 'mfrow' and 'mfcol' for multi figure plotting. subplot

Examples

old.par <- par(no.readonly=TRUE)

par(mfrow=c(2,2),xpd=NA)

# generate some sample data
tmp.x <- rnorm(25, 10, 2)
tmp.y <- rnorm(25, 50, 10)
tmp.z <- rnorm(25, 0, 1)

plot( tmp.x, tmp.y)

# draw a diagonal line across the plot area
tmp1 <- cnvrt.coords( c(0,1), c(0,1), input='plt' )
lines(tmp1$usr, col='blue')

# draw a diagonal line accross figure region
tmp2 <- cnvrt.coords( c(0,1), c(1,0), input='fig')
lines(tmp2$usr, col='red')

# save coordinate of point 1 and y value near top of plot for future plots
tmp.point1 <- cnvrt.coords(tmp.x[1], tmp.y[1])
tmp.range1 <- cnvrt.coords(NA, 0.98, input='plt')

# make a second plot and draw a line linking point 1 in each plot
plot(tmp.y, tmp.z)

tmp.point2 <- cnvrt.coords( tmp.point1$dev, input='dev' )
arrows( tmp.y[1], tmp.z[1], tmp.point2$usr$x, tmp.point2$usr$y,
 col='green')

# draw another plot and add rectangle showing same range in 2 plots

plot(tmp.x, tmp.z)
tmp.range2 <- cnvrt.coords(NA, 0.02, input='plt')
tmp.range3 <- cnvrt.coords(NA, tmp.range1$dev$y, input='dev')
rect( 9, tmp.range2$usr$y, 11, tmp.range3$usr$y, border='yellow')

# put a label just to the right of the plot and
#  near the top of the figure region.
text( cnvrt.coords(1.05, NA, input='plt')$usr$x,
        cnvrt.coords(NA, 0.75, input='fig')$usr$y,
        "Label", adj=0)

par(mfrow=c(1,1))

## create a subplot within another plot (see also subplot)

plot(1:10, 1:10)

tmp <- cnvrt.coords( c( 1, 4, 6, 9), c(6, 9, 1, 4) )

par(plt = c(tmp$dev$x[1:2], tmp$dev$y[1:2]), new=TRUE)
hist(rnorm(100))

par(fig = c(tmp$dev$x[3:4], tmp$dev$y[3:4]), new=TRUE)
hist(rnorm(100))

par(old.par)


Designs for coin faces for use with plot.rgl.coin

Description

This is a list of matricies where each matrix represents a design for drawing lines on the face of a coin.

Usage

data(coin.faces)

Format

The format is: List of 4 $ qh: num [1:57, 1:2] 0.387 0.443 0.515 0.606 0.666 ... $ qt: num [1:62, 1:2] 0.862 0.873 0.875 0.857 0.797 ... $ H : num [1:28, 1:2] 0.503 0.506 0.548 0.548 0.500 ... $ T : num [1:18, 1:2] 0.506 0.520 0.569 0.626 0.626 ...

Details

The current options are a capitol "H", a capitol "T", a design representing George Washingtons head traced from the heads of a US quarter, and a design representing an eagle traced from the tails of a US quarter.

The tracings here have pretty much exhausted my artistic ability, if you can do better, please do, I will be happy to include it in future versions. It would also be nice to include some designs representing faces of non-US coins, please submit your contributions (the design should fit within a circle inscribed within the unit square).

Examples

## Not run:
plot.rgl.coin(heads=coin.faces$H, tails=coin.faces$T)
## End(Not run)


Simulate rolling dice

Description

Simulate and optionally plot rolls of dice.

Usage

dice(rolls = 1, ndice = 2, sides = 6, plot.it = FALSE, load = rep(1, sides))
plot.dice(x, ...)

Arguments

rolls | Scalar, the number of times to roll the dice.
ndice | Scalar, the number of dice to roll each time.
sides | Scalar, the number of sides per die.
plot.it | Logical, Should the results be plotted.
load | Vector of length sides, how the dice should be loaded.
x | Data frame, return value from dice.
... | Additional arguments passed to lattice plotting function.

Details

Simulates the rolling of dice. By default it will roll 2 dice 1 time and the dice will be fair. Internally the sample function is used and the load option is passed to sample. load is not required to sum to 1, but the elements will be divided by the sum of all the values.

Value

A data frame with rolls rows and ndice columns representing the results from rolling the dice.
If only 1 die is rolled, then the return value will be a vector.
If plot.it is TRUE, then the return value will be invisible.

Note

If the plot function is used or if plot.it is TRUE, then a plot will be created on the current graphics device.

Author(s)

Greg Snow greg.snow@intermountainmail.org

See Also

sample

Examples

# 10 rolls of 4 fair dice
dice(10,4, plot.it=TRUE)

# or

plot(dice(10,4))

# or

tmp <- dice(10,4)
plot(tmp)

# a loaded die
table(tmp <- dice(100,1,plot.it=TRUE, load=6:1 ) )
mean(tmp)

## redo De Mere's question

demere1 <- dice(10000,4)
demere2 <- dice(10000,24,sides=36)

mean(apply( demere1, 1, function(x) 6 %in% x ))

mean(apply( demere2, 1, function(x) 36 %in% x))

plot(demere1[1:10,])

## plot all possible combinations of 2 dice

plot.dice( expand.grid(1:6,1:6), layout=c(6,6) )


Create a quick dotchart (histogram)

Description

Create a quick dotchart of 1 or 2 datasets. These dotcharts are a poor man's histogram, not the trellis dotplot.

Usage

dots(x,...)
dots2(x, y, colx = "green", coly = "blue", lab1 =
deparse(substitute(x)), lab2 = deparse(substitute(y)),...)

Arguments

x | Vector, data to be plotted (should be rounded).
y | Vector, second dataset to be plotted.
colx | Color of points for x.
coly | Color of points for y.
lab1 | Label for x.
lab2 | Label for y.
... | Additional arguments passed to plotting functions.

Details

These functions create basic dotcharts that are quick "back of the envelope" approximations to histograms. Mainly intended for demonstration.

Value

No meaninful value. These functions are run for the side effect of creating a plot.

Author(s)

Greg Snow greg.snow@intermountainmail.org

See Also

dotplot in the lattice package, hist

Examples

dots( round( rnorm(50, 10,3) ) )
dots2( round( rnorm(20, 10,3) ), round(rnorm(20,12,2)) )


Chernoff Faces

Description

faces represent the rows of a data matrix by faces

Usage

faces(xy, which.row, fill = FALSE, nrow, ncol, scale = TRUE, byrow = FALSE, main, labels)

Arguments

xy | xy data matrix, rows represent individuals and columns attributes
which.row | defines a permutation of the rows of the input matrix
fill | if(fill==TRUE), only the first nc attributes of the faces are transformed, nc is the number of columns of xy
nrow | number of columns of faces on graphics device
ncol | number of rows of faces
scale | if(scale==TRUE), attributes will be normalized
byrow | if(byrow==TRUE), xy will be transposed
main | title
labels | character strings to use as names for the faces

Details

The features paramters of this implementation are: 1-height of face, 2-width of face, 3-shape of face, 4-height of mouth, 5-width of mouth, 6-curve of smile, 7-height of eyes, 8-width of eyes, 9-height of hair, 10-width of hair, 11-styling of hair, 12-height of nose, 13-width of nose, 14-width of ears, 15-height of ears. For details look at the literate program of faces

Value

a plot of faces is created on the graphics device, no numerical results

Note

version 12/2003

Author(s)

H. P. Wolf

References

Chernoff, H. (1973): The use of faces to represent statistiscal assoziation, JASA, 68, pp 361–368. The smooth curves are computed by an algorithm found in Ralston, A. and Rabinowitz, P. (1985): A first course in numerical analysis, McGraw-Hill, pp 76ff. http://www.wiwi.uni-bielefeld.de/~wolf/ : S/R - functions : faces

See Also

—

Examples

##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--    or do  help(data=index)  for the standard data sets.

faces(rbind(1:3,5:3,3:5,5:7))

data(longley)
faces(longley[1:9,])

set.seed(17)
faces(matrix(sample(1:1000,128,),16,8),main="random faces")

## The function is currently defined as
function(xy=rbind(1:3,5:3,3:5,5:7),which.row,fill=FALSE,nrow,ncol,
                scale=TRUE,byrow=FALSE,main,labels){

spline<-function(a,y,m=200,plot=FALSE){
    n<-length(a)
  h<-diff(a)
  dy<-diff(y)
  sigma<-dy/h
  lambda<-h[-1]/(hh<-h[-1]+h[-length(h)])
  mu<-1-lambda
  d<-6*diff(sigma)/hh
  tri.mat<-2*diag(n-2)
  tri.mat[2+  (0:(n-4))*(n-1)] <-mu[-1]
  tri.mat[    (1:(n-3))*(n-1)] <-lambda[-(n-2)]
  M<-c(0,solve(tri.mat)%*%d,0)
  x<-seq(from=a[1],to=a[n],length=m)
  anz.kl <- hist(x,breaks=a,plot=FALSE)$counts
  adj<-function(i) i-1
  i<-rep(1:(n-1),anz.kl)+1
  S.x<-  M[i-1]*(a[i]-x          )^3 / (6*h[adj(i)])  +
         M[i]  *(x        -a[i-1])^3 / (6*h[adj(i)])  +
         (y[i-1] - M[i-1]*h[adj(i)]^2 /6) * (a[i]-x)/ h[adj(i)] +
        (y[i]   - M[i]  *h[adj(i)]^2 /6) * (x-a[i-1]) / h[adj(i)]
  if(plot){ plot(x,S.x,type="l"); points(a,y)    }
  return(cbind(x,S.x))
  }

n.char<-15
xy<-rbind(xy)
if(byrow) xy<-t(xy)
if(!missing(which.row)&& all(  !is.na(match(which.row,1:dim(xy)[2]))  ))
       xy<-xy[,which.row,drop=FALSE]
mm<-dim(xy)[2];  n<-dim(xy)[1]
xnames<-dimnames(xy)[[1]]
if(is.null(xnames)) xnames<-as.character(1:n)
if(!missing(labels)) xnames<-labels
if(scale){
   xy<-apply(xy,2,function(x){
           x<-x-min(x); x<-if(max(x)>0) 2*x/max(x)-1 else x })
} else xy[]<-pmin(pmax(-1,xy),1)
xy<-rbind(xy);n.c<-dim(xy)[2]
xy<-xy[,(h<-rep(1:mm,ceiling(n.char/mm))),drop=FALSE]
if(fill) xy[,-(1:n.c)]<-0

face.orig<-list(
      eye  =rbind(c(12,0),c(19,8),c(30,8),c(37,0),c(30,-8),c(19,-8),c(12,0))
     ,iris =rbind(c(20,0),c(24,4),c(29,0),c(24,-5),c(20,0))
     ,lipso=rbind(c(0,-47),c( 7,-49),lipsiend=c( 16,-53),c( 7,-60),c(0,-62))
     ,lipsi=rbind(c(7,-54),c(0,-54))                  # add lipsiend
     ,nose =rbind(c(0,-6),c(3,-16),c(6,-30),c(0,-31))
     ,shape =rbind(c(0,44),c(29,40),c(51,22),hairend=c(54,11),earsta=c(52,-4),
                  earend=c(46,-36),c(38,-61),c(25,-83),c(0,-89))
     ,ear  =rbind(c(60,-11),c(57,-30))                # add earsta,earend
     ,hair =rbind(hair1=c(72,12),hair2=c(64,50),c(36,74),c(0,79)) # add hairend
)
lipso.refl.ind<-4:1
lipsi.refl.ind<-1
nose.refl.ind<-3:1
hair.refl.ind<-3:1
shape.refl.ind<-8:1
shape.xnotnull<-2:8
nose.xnotnull<-2:3

nr<-n^0.5; nc<-n^0.5
if(!missing(nrow)) nr<-nrow
if(!missing(ncol)) nc<-ncol
opar<-par(mfrow=c(ceiling(c(nr,nc))),oma=rep(6,4), mar=rep(.7,4))
on.exit(par(opar))

for(ind in 1:n){

factors<-xy[ind,]
face <- face.orig

m<-mean(face$lipso[,2])
face$lipso[,2]<-m+(face$lipso[,2]-m)*(1+0.7*factors[4])
face$lipsi[,2]<-m+(face$lipsi[,2]-m)*(1+0.7*factors[4])
face$lipso[,1]<-face$lipso[,1]*(1+0.7*factors[5])
face$lipsi[,1]<-face$lipsi[,1]*(1+0.7*factors[5])
face$lipso["lipsiend",2]<-face$lipso["lipsiend",2]+20*factors[6]

m<-mean(face$eye[,2])
face$eye[,2] <-m+(face$eye[,2] -m)*(1+0.7*factors[7])
face$iris[,2]<-m+(face$iris[,2]-m)*(1+0.7*factors[7])
m<-mean(face$eye[,1])
face$eye[,1] <-m+(face$eye[,1] -m)*(1+0.7*factors[8])
face$iris[,1]<-m+(face$iris[,1]-m)*(1+0.7*factors[8])

m<-min(face$hair[,2])
face$hair[,2]<-m+(face$hair[,2]-m)*(1+0.2*factors[9])
m<-0
face$hair[,1]<-m+(face$hair[,1]-m)*(1+0.2*factors[10])
m<-0
face$hair[c("hair1","hair2"),2]<-face$hair[c("hair1","hair2"),2]+50*factors[11]

m<-mean(face$nose[,2])
face$nose[,2]<-m+(face$nose[,2]-m)*(1+0.7*factors[12])
face$nose[nose.xnotnull,1]<-face$nose[nose.xnotnull,1]*(1+factors[13])

m<-mean(face$shape[c("earsta","earend"),1])
face$ear[,1]<-m+(face$ear[,1]-m)* (1+0.7*factors[14])
m<-min(face$ear[,2])
face$ear[,2]<-m+(face$ear[,2]-m)* (1+0.7*factors[15])

face<-lapply(face,function(x){ x[,2]<-x[,2]*(1+0.2*factors[1]);x})
face<-lapply(face,function(x){ x[,1]<-x[,1]*(1+0.2*factors[2]);x})
face<-lapply(face,function(x){ x[,1]<-ifelse(x[,1]>0,
                                        ifelse(x[,2] > -30, x[,1],
                  pmax(0,x[,1]+(x[,2]+50)*0.2*sin(1.5*(-factors[3])))),0);x})

invert<-function(x) cbind(-x[,1],x[,2])
face.obj<-list(
     eyer=face$eye
    ,eyel=invert(face$eye)
    ,irisr=face$iris
    ,irisl=invert(face$iris)
    ,lipso=rbind(face$lipso,invert(face$lipso[lipso.refl.ind,]))
    ,lipsi=rbind(face$lipso["lipsiend",],face$lipsi,
                 invert(face$lipsi[lipsi.refl.ind,,drop=FALSE]),
                 invert(face$lipso["lipsiend",,drop=FALSE]))
    ,earr=rbind(face$shape["earsta",],face$ear,face$shape["earend",])
    ,earl=invert(rbind(face$shape["earsta",],face$ear,face$shape["earend",]))
    ,nose=rbind(face$nose,invert(face$nose[nose.refl.ind,]))
    ,hair=rbind(face$shape["hairend",],face$hair,invert(face$hair[hair.refl.ind,]),
                invert(face$shape["hairend",,drop=FALSE]))
    ,shape=rbind(face$shape,invert(face$shape[shape.refl.ind,]))
)

plot(1,type="n",xlim=c(-105,105)*1.1, axes=FALSE,
     ylab="",ylim=c(-105,105)*1.3)
title(xnames[ind])
for(ind in seq(face.obj)) {
       x <-face.obj[[ind]][,1]; y<-face.obj[[ind]][,2]
       xx<-spline(1:length(x),x,40,FALSE)[,2]
       yy<-spline(1:length(y),y,40,FALSE)[,2]
       lines(xx,yy)
  }
  }

if(!missing(main)){
  par(opar);par(mfrow=c(1,1))
  mtext(main, 3, 3, TRUE, 0.5)
  title(main)
  }
  }


Chernoff Faces

Description

Plot Chernoff Faces of the dataset, rows represent subjects/observations, columns represent variables.

Usage

faces2(mat, which = 1:ncol(mat), labels = rownames(mat),
  nrows = ceiling(nrow(mat)/ncols), ncols = ceiling(sqrt(nrow(mat))),
  byrow = TRUE, scale = c("columns", "all", "center", "none"),
  fill = c(0.5, 0.5, 1, 0.5, 0.5, 0.3, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
           0.5, 0.5, 0.5, 0.5, 1, 0.5), ...)

Arguments

mat | Matrix containing the data to plot.
which | Which columns correspond to which features (see details).
labels | Labels for the individual faces
nrows | Number of rows in the graphical layout
ncols | Number of columns in the graphical layout
byrow | Logical, should the faces be drawn rowwise or columnwise.
scale | Character, how should the data be scaled.
fill | What value to use for features not assocaiated with a column of data.
... | Additional arguments passed on to plotting functions.

Details

The features are: 1 Width of center 2 Top vs. Bottom width (height of split) 3 Height of Face 4 Width of top half of face 5 Width of bottom half of face 6 Length of Nose 7 Height of Mouth 8 Curvature of Mouth (abs < 9) 9 Width of Mouth 10 Height of Eyes 11 Distance between Eyes (.5-.9) 12 Angle of Eyes/Eyebrows 13 Circle/Ellipse of Eyes 14 Size of Eyes 15 Position Left/Right of Eyeballs/Eyebrows 16 Height of Eyebrows 17 Angle of Eyebrows 18 Width of Eyebrows

The face plotting routine needs the data values to be between 0 and 1 (inclusive). The scale option controls how scaling will be done on mat: "columns" scales each column to range from 0 to 1, "all" scales the entire dataset to vary from 0 to 1, "center" scales each column so that the mean of the column becomes 0.5 and all other values are between 0 and 1, and "none" does no scaling assuming that the data has already been scaled.

Value

This function is run for its side effect of plotting and does not return anything.

Note

If you choose to not scale the data and any data values are outside of the 0 to 1 range, then strange things may happen.

This function is based on code available at , the good things come from there, any problems are likely due to my (Greg's) tweaking.

Author(s)

Original code by ; current implementation by Greg Snow greg.snow@intermountainmail.org

References

Chernoff, H. (1973): The use of faces to represent statistiscal assoziation, JASA, 68, pp 361–368.

See Also

faces

Examples

  faces2(matrix( runif(18*10), nrow=10), main='Random Faces')


Create a Fagan plot to demonstrate Bayes Theorem and screening tests

Description

These functions create a plot showing the relationship between the prior probability, the LR (combination of sensitivity and specificity), and the posterior probability.

Usage

fagan.plot(probs.pre.test, LR, test.result="+")
plotFagan()
plotFagan2()

Arguments

probs.pre.test | The prior probability
LR | the likelihood ratio (sensitivity/(1-specificity))
test.result | either '+' or '-' indicating whether you want the probability of the event or of not seeing the event

Details

When Bayes theorem is expressed in terms of log-odds it turns out that the posterior log-odds are a linear function of the prior log-odds and the log likelihood ratio. These functions plot an axis on the left with the prior log-odds, an axis in the middle representing the log likelihood ratio and an axis on the right representing the posterior log-odds. A line is then drawn from the prior probability on the left through the LR in the center and extended to the posterior probability on the right. The fagan.plot creates the plot based on input to the function. The plotFagan and plotFagan2 functions set up Tk windows with sliders representing the possible inputs and show how the plot and the posterior probability changes when you adjust the inputs. The plotFagan function creates sliders for the prior probability and the LR, while the plotFagan2 function replaces the LR slider with 2 sliders for the sensitivity and specificity.

More detail on the plots and the math behind them can be found at the websites below.

Value

These functions are run for their side effects and do not return a meaningful value

Author(s)

Guazzetti Stefano and Greg Snow greg.snow@intermountainmail.org

References

Fagan TJ. Nomogram for Bayes theorem. N Engl J Med 1975;293(5):257-61. http://www.cmh.edu/stats/definitions/fagan.htm http://ebm.bmjjournals.com/cgi/content/full/6/6/164

See Also

slider

Examples

## Not run:
 fagan.plot(0.8, 2)

 fagan.plot(0.8, 0.95/(1-0.90) )

 plotFagan()

 plotFagan2()
## End(Not run)


Alpha version functions to send plotting commands to GnuPlot

Description

These functions allow you to open a connection to a gnuplot process, send data and possibly other information to gnuplot for it to plot, then close gnuplot and clean up temporary files and variables. These functions are alpha level at best, use at your own risk.

Usage

gp.open(where='c:/progra~1/GnuPlot/bin/pgnuplot.exe')
gp.close(pipe=.gp)
gp.send(cmd='replot',pipe=.gp)
gp.plot(x,y,type='p',add=FALSE, title=deparse(substitute(y)),pipe=.gp)

Arguments

where | Path to GnuPlot Executable
pipe | The pipe object connected to GnuPlot (returned from gp.open), warning: changing this from the default will probably break things
cmd | Text string, the command to be sent verbatim to the GnuPlot process
x | The x coordinates to plot
y | the y coordinates to plot
type | Either 'p' or 'l' for plotting points or lines
add | Logical, should the data be added to the existing plot or start a new plot
title | The title or legend entry

Details

These functions provide a basic interface to the GnuPlot program (you must have GnuPlot installed (separate install)), gp.open runs GnuPlot and establishes a pipe connection, gp.close sends a quite command to gnuplot and cleans up temporary variables and files, gp.send sends a command to the GnuPlot process verbatim, and gp.plot sends data and commands to the process to create a standard scatterplot or line plot.

Value

gp.open returns and invisible copy of the pipe connection object (to pass to other functions, but don't do this because it doesn't work right yet).
The other 3 functions don't return anything meaningful. All functions are run for their side effects.

Note

These functions create some temporary files and 2 temporary global variables (.gp and .gp.tempfiles), running gp.close will clean these up (so use it).

These functions are still alpha level.

Author(s)

Greg Snow greg.snow@intermountainmail.org

References

http://www.gnuplot.info/

See Also

plot

Examples

 ## Not run:

x <- 1:10
y <- 3-2*x+x*x+rnorm(10)

gp.open()
gp.plot(x,y)
gp.send('replot 3-2*x+x**2')
gp.close()
 ## End(Not run)


Compute Highest Posterior Density Intervals

Description

Compute the Highest Posterior Density Interval (HPD) from an inverse density function (hpd) or a vector of realizations of the distribution (emp.hpd).

Usage

hpd(posterior.icdf, conf=0.95, tol=0.00000001,...)

emp.hpd(x, conf=0.95)

Arguments

posterior.icdf | Function, the inverse cdf of the posterior distribution (usually a function whose name starts with 'q').
x | A vector of realizations from the posterior distribution.
conf | Scalar, the confidence level desired.
tol | Scalar, the tolerance for optimize.
... | Additional arguments to posterior.icdf.

Details

These functions compute the highest posterior density intervals (sometimes called minimum length confidence intervals) for a Bayesian posterior distribution. The hpd function is used when you have a function representing the inverse cdf (the common case with conjugate families). The emp.hpd function is used when you have realizations of the posterior (when you have results from an MCMC run).

Value

A vector of length 2 with the lower and upper limits of the interval.

Note

These functions assume that the posterior distribution is unimodal, they compute only 1 interval, not the set of intervals that are appropriate for multimodal distributions.

Author(s)

Greg Snow greg.snow@intermountainmail.org

See Also

hdr in the hdrcde package.

Examples

hpd(qbeta, shape1=50, shape2=250)

tmp <- rbeta(10000, 50, 250)
emp.hpd(tmp)


Identify polygons within a plotted map from the maptools packages.

Description

This function allows you to click in a polygon in a map plotted using the maptools package and labels the polygon.

Usage

identify.Map(x, labels = as.character(Map$att.data$NAME), n = Inf, ...)

identify.polylist(...)

Arguments

x | The map object plotted using plot.Map.
labels | Vector of labels used to label the polygons.
n | Maximum number of polygons to label.
... | Additional arguments passed to identify.Map or text.

Details

The function identify.polylist is just a wrapper function for identify.Map which handles both types of data.

The point you click on is checked against every polygon in the map, so this function may be slow on very complicated maps.

Right click on the map to stop before n.

Value

A vector of integers corresponding to the polygons clicked on.

Note

These functions are beta level at best and will probably soon be seperceeded by functions in maptools, sp, or spproj. The plot.Map function has been depricated so use plot.polylist instead (convert Map objects to polylist objects before plotting), identify.Map works on polylist objects as well.

Author(s)

Greg Snow greg.snow@intermountainmail.org

See Also

plot.Map in the maptools package, identify

Examples

## Not run:
library(maptools)
plot(state.vbm)
identify(state.vbm) # now click on the map a few times.

plot(state.vbm)
identify(state.vbm, state.abb, n=5)
## End(Not run)


Interactively explore the conditioned panels in lattice plots.

Description

Plot 1 panel from an xyplot, and optionally a 3d graph highligting the shown points, then allow you to interactively set the conditioning set of data to see the effects and help you better understand how xyplot works.

Usage

lattice.demo(x, y, z, show3d = TRUE)

Arguments

x | The x variable to plot (numeric).
y | The y variable to plot (numeric).
z | The variable to condition on (numeric).
show3d | Logical, should a 3D cloud be shown as well.

Details

This function is intended to for demonstration purposes to help understand what is happening in an xyplot (lattice). When you run the demo it will create a single panel from a conditioned xyplot and optionally a 3D cloud with the points included in the panel highlighted. The function then opens a tcl/tk dialog box that allows you to choose which points are included in the panel (based on the conditioning variable). You can choose the center and width of the shingle displayed and the graph will update to show the new selection.

The intent for this function is for a teacher to show a class how lattice graphics take slices of a 3d plot and show each slice seperately. Students could then work through some examples on their own to better understand what functions like xyplot are doing automatically.

Value

No meaningful return value, this function is run for the side effects.

Author(s)

Greg Snow greg.snow@intermountainmail.org

See Also

xyplot

Examples

## Not run:
require(stats)
lattice.demo(quakes$long, quakes$lat, quakes$depth)
## End(Not run)


Functions to plot Shewhart g chart

Description

These functions are used to compute statistics required by the g chart (geometric distribution) for use with the qcc package.

Usage

stats.g(data, sizes)
sd.g(data, sizes)
limits.g(center, std.dev, sizes, conf)

Arguments

data | the observed data values
center | sample center statistic
sizes | sample sizes (not used)
std.dev | standard deviation of geometric distribution
conf | a numeric value used to compute control limits, specifying the number of standard deviations (if 'conf' > 1) or the confidence level (if 0 < 'conf' < 1).

Details

These functions work with the qcc package and may someday be moved into the qcc package.

The g chart plots the number of non-events between events. np charts do not work well when the probability of an event is rare (see example below). Instead of plotting the number of events, the g chart plots the number of non-events between events.

Value

The function 'stats.g' returns a list with components 'statistics' and 'center'.
The function 'sd.g' returns 'std.dev' the standard deviation (sqrt(1-p)/p).
The function 'limits.g' returns a matrix with lower and upper control limits.

Note

The geometric distribution is quite skewed so it is best to set conf at the required confidence interval (0 < conf < 1) rather than as a multiplier of sigma.

Author(s)

Greg Snow greg.snow@intermountainmail.org

References

Kaminsky, FC et. al. (1992): Statistical Control Charts Based on a Geometric Distribution, Journal of Quality Technology, 24, pp 63–69.

Yang, Z et. al. (2002): On the Performance of Geometric Charts with Estimated Control Limits, Journal of Quality Technology, 34, pp 448–458.

See Also

qcc

Examples

if(require(qcc)){
  success <- rbinom(1000, 1, 0.01)
  num.noevent <- diff(which(c(1,success)==1))-1
  qcc(success, type='np', sizes=1)
  qcc(num.noevent, type="g")
}


Demonstrate the internals of loess curve fits

Description

Creates a scatterplot with a loess fit, then interactively shows the window and case weights used to create the curve at the selected value of x.

Usage

loess.demo(x, y, span = 2/3, degree = 1, nearest = F, xlim = numeric(0), ylim = numeric(0), verbose = F)

Arguments

x | The x coordinates to be plotted.
y | The y coordinates to be plotted.
span | The relative width of the window, passed on to loess.
degree | Degree of polynomial to use (0, 1, or 2), passed on to loess.
nearest | Logical, should predictions be made at the point where you clicked (FALSE), or at the nearest x value of the data to where you clicked (TRUE).
xlim | Limits of the Horizonal axis.
ylim | Limits of the Vertical axis.
verbose | If true then print the x coordinate being predicted.

Details

This function demonstrates the underlying calculations of loess curves.

Given x and y vectors it will create a scatterplot and add 2 loess fit lines (one using straight loess smooth with linear interpolation and one that does a spline interpolation of the loess fit).

The function then waits for the user to click on the plot. The function then shows the window of points (centered at the x value clicked on) used in the weighting for predicting that point and shows a circle around each point in the window where the area of the circle is proportional to the weight of that point in the linear fit. The function also shows the linear (or quadratic) fit used to predict at the selected point.

The basic steps of the loess algorithm (as demonstrated by the function) is that to predict the y-value for a given x-value the computer:

1. Find all the points within a window around the x-value (the width of the window is based on the parameter span). 2. Weight the points in the window with points nearest the x-value having the highest weight. 3. Fit a weighted linear (quadratic) line to the points in the window. 4. Use the y-value of the fitted line (curve) at the x-value to give loess prediction at that x-value.

Clicking on another point in the graph will replot with the new situation.

Right click and select 'stop' to end the demonstration.

Value

This function does not return anything, it is run purely for its side effects.

Author(s)

Greg Snow greg.snow@intermountainmail.org

See Also

loess, locator

Examples

## Not run:
data(ethanol, package='lattice')
attach(ethanol)
loess.demo(E, NOx)
# now click a few places, right click to end
loess.demo(E, NOx, span=1.5)
loess.demo(E, NOx, span=0.25)
loess.demo(E, NOx, degree=0)
loess.demo(E, NOx, degree=2)
detach()
## End(Not run)


Demonstrate the basic concept of Maximum Likelihood Estimation

Description

This function graphically shows log likelihoods for a set of data and the normal distribution and allows you to interactively change the parameter estimates to see the effect on the log likelihood.

Usage

mle.demo(x = rnorm(10, 10, 2), start.mean = mean(x) - start.sd, start.sd = 1.2 * sqrt(var(x)))

Arguments

x | A vector of data (presumably from a normal distribution).
start.mean | The initial value for estimating the mean.
start.sd | The initial value for estimating the standard deviation.

Details

The function creates a plot with 3 panels: the top panel shows a normal curve based on the current values of the mean and standard deviation along with a vertical line for each point in x (the product of the heights of these lines is the likelihood, the sum of the logs of their heights is the log likelihood).

The lower 2 plots show the profiles of the mean and standard deviation. The y-axis is the likelihoods of the parameters tried so far, and the x-axes are the mean and standard deviation tried. The point corresponding to the current parameter estimates will be solid red.

A Tk slider box is also created that allows you to change the current estimates of the mean and standard deviation to show the effect on the log likelihood and find the maximum likelihood estimate.

Value

This function is run for its side effects and returns NULL.

Author(s)

Greg Snow greg.snow@intermountainmail.org

See Also

fitdistr in package MASS, mle in package stats4, slider

Examples

## Not run:
mle.demo()

m <- runif(1, 50,100)
s <- runif(1, 1, 10)
x <- rnorm(15, m, s)

mm <- mean(x)
ss <- sqrt(var(x))
ss2 <- sqrt(var(x)*11/12)
mle.demo(x)
# now find the mle from the graph and compare it to mm, ss, ss2, m, and s
## End(Not run)


Animated die roll or coin flip

Description

Open an rgl window, plot either a representation of a coin or a die then animate the flipping/rolling.

Usage

plot.rgl.coin(x, col = "black", heads = x[[1]], tails = x[[2]], ...)

plot.rgl.die(x=1:6, col.cube = "white", col.pip = "black", sides = x, ...)

flip.rgl.coin(side = sample(2, 1), steps = 150)

roll.rgl.die(side = sample(6, 1), steps = 250)

Arguments

x | for plot.rgl.coin a list with information for drawing the faces of the coin, defaults to coin.faces. For plot.rgl.die a vector with the number of pips to put on the sides of the die (alternative way of specifying sides).
col | Color of lines on the coin faces.
heads | Design to use as "heads" side of coin.
tails | Design to use as "tails" side of coin.
col.cube | Color of the cube for the die.
col.pip | Color of the pips (spots) on the die
sides | Vector of length 6 indicating which numbers to show on the die.
side | Which side of the coin (1 or 2) or die (1 through 6) should end up face up.
steps | The number of steps in each part of the animation, higher values will be smoother and slower, lower values will be faster but more jumpy.
... | Currently any additional options are silently ignored.

Details

You must use the plot function first to create the coin or die, then use the flip or roll function to start the animation. You can animate multiple times for a single use of the plotting function.

You can manually rotate the image as well, see the rgl package for details.

The defaults plot a regular coin and die, but arguments are available to create special casses (2 headed coin, die with 2 6's and no 1, ...).

The data list coin.faces contains information on designs for the faces of the coins in case you want to choose a different design.

The default rolling and flipping options ranomly choose which side will be face up following a uniform distribution. You can specify the side yourself, or use the sample function to do a biased random flip/roll.

Value

Which side ended up face up (1 or 2 for coin, 1 through 6 for die). This is the internal numbering and does not match a change in the sides argument.

Note

The current algorithm for animating the die roll shows all the sides, but I am not satisfied with it. Please suggest improvements.

Author(s)

Greg Snow greg.snow@intermountainmail.org

See Also

dice, plot.dice, coin.faces, sample

Examples

## Not run:
plot.rgl.coin()
flip.rgl.coin()
flip.rgl.coin(1)
flip.rgl.coin(2)

rgl.clear()

# two-headed coin
plot.rgl.coin(tails=coin.faces$qh)

rgl.clear()

# letters instead of pictures
plot.rgl.coin(heads=coin.faces$H, tails=coin.faces$T)

# biased flip
flip.rgl.coin( sample(2,1, prob=c(0.65, 0.35) ) )

rgl.clear()

plot.rgl.die()
roll.rgl.die()
roll.rgl.die(6)

# biased roll
roll.rgl.die( sample(6,1, prob=c(1,2,3,3,2,1) ) )
## End(Not run)


Graphically illustrate the concept of power.

Description

Create graphs of a normal test statistic under the null and alternative hypotheses to graphically show the idea of power.

Usage

power.examp(n = 1, stdev = 1, diff = 1, alpha = 0.05, xmin = -2, xmax = 4)
run.power.examp()

Arguments

n | The sample size for the test statistic.
stdev | The standard deviation of the population.
diff | The true difference in means (alternate hypothesis).
alpha | The type I error rate to use for the test.
xmin | The minimum x value to show on the graph.
xmax | The maximum x value to show on the graph.

Details

This function will draw 2 graphs representing an upper-tailed test of hypothesis.

The upper panel represents the test statistic under the null hypothesis that the true mean (or mean difference) is 0. It then also shows the upper tail area equal to alpha and the rejection region for the test statistic.

The lower panel shows the normal distribution for the test statistic under the alternative hypothesis where the true mean (or mean difference) is diff. Using the rejection region from the upper panel it shades the upper tail area that corresponds to the power of the test.

Both curves are affected by the specified stdev and sample size n.

The function run.power.examp will in addition create a Tk slider box that will allow you to interactively change the values of stdev, diff, alpha, and n to dynamically see the effects of the change on the graphs and on the power of the test.

This can be used to demonstrate the concept of power, show the effect of sample size on power, show the inverse relationship between the type I and type II error rates, and show how power is dependent on the true mean (or difference) and the population standard deviation.

Value

power.examp invisibly returns the power computed.
run.power.examp does not return anything meaningful.

Author(s)

Greg Snow greg.snow@intermountainmail.org

See Also

power.t.test

Examples

power.examp()
power.examp(n=25)
power.examp(alpha=0.1)


Apply projections to Map objects (from maptools package)

Description

Apply the mapproject function from the mapproj library to Map objects from the maptools library.

Usage

project.Map(Map, projection = "", parameters = NULL, orientation = c(90, 0, 0))

Arguments

Map | A Map object.
projection | See mapproject
parameters | See mapproject
orientation | See mapproject

Details

This function takes a Map object and applies mapproject from the mapproject library to it.

Value

An object of class Map with the coordinates transformed.

Note

This is beta level software and will someday be superceded by the spproj package.

Author(s)

Greg Snow greg.snow@intermountainmail.org

See Also

mapproject from the mapproj library.

Examples

## Not run:

library(mapproj)
# assumes that the time zone shape files have been downloaded
# from: http://openmap.bbn.com/data/shape/timezone/

tz <- read.shape('WRLDTZA')
plot(Map2poly(project.Map(tz,'bonne',param=45)))

plot(Map2poly(project.Map(tz,'bonne',param=45)),
  col=rainbow(48)[tz$att.data$OFFSET])

plot(Map2poly(project.Map(tz,'albers',param=c(30,40))))
## End(Not run)


Demonstrate Correlation and Regression by placing and moving data points

Description

Place data points on a graph to demonstrate concepts related to correlation and regression.

Usage

put.points.demo(x = NULL, y = NULL, lsline = TRUE)

Arguments

x | x-coordinates for initial points.
y | y-coordinates for initial points.
lsline | Logical, should the ls regresion line be included.

Details

The plot area is divided into 2 sections, the left section shows a scatterplot of your points, the right panel controls what happens when you click in the left panel.

The top of the right panel has an "end" button that you click on to end the demonstration.

The middle right panel toggles the least squares line and information.

The bottom right panel has radio buttons that determine what clicking in the left panel will do, the options are to add a point, delete a point, or move a point.

To move a point click on the point you want to move, it will become solid, then click in the place you want it to move to.

When deleting or moving points, the closest point to where you click will be deleted or moved, even if you click in an empty area.

Whenever you add, delete, or move a point the correlation, r^2, and regression line will be updated. You can start with a set of points then demonstrate what happens to the correlation and regression line when outliers are added or important points are moved or deleted.

Value

This function does not return anything.

Author(s)

Greg Snow greg.snow@intermountainmail.org

See Also

plot, cor

Examples

## Not run:
put.points.demo()

x <- rnorm(25, 5, 1)
y <- x + rnorm(25)
put.points.demo(x,y)
## End(Not run)


Recenter a Map object from the maptools package.

Description

Most maps are centered at a longitude of 0 with points of east longitude plotted on the right side of the map and points with west longitude plotted to the left. This function transforms the points to have a new center moving polygons from the right to left or left to right.

Usage

recenter.Map(Map, center = 0, range = c(-180, 180))

Arguments

Map | A Map object from the maptools package.
center | The new longitude to use as the center.
range | The range of possible longitudes on the original scale.

Details

Plotting a map of the United States can place some of the Alaskan islands to the far right of the map (they have east longitude). This function shifts the longitudes to bring the seperate pieces together.

Polygons will only be shifted if all points can move together, it will not split a polygon between sides of the map.

Value

An object of class Map (from the maptools package) that is a copy of the input Map but with the longitudes shifted.

Note

This is beta level software and will probably be superceded by the spproj or sp packages. The function plot.Map has been depricated, use plot.polylist instead (use Map2poly to convert).

Author(s)

Greg Snow greg.snow@intermountainmail.org

See Also

plot.polylist

Examples

## Not run:
## this assumes that the time zone shape files have been downloaded
## from: http://openmap.bbn.com/data/shape/timezone/
## and the US states map (st99_d00) files have been downloaded from
## http://www.census.gov/geo/www/maps/

us <- read.shape('st99_d00')
tz <- read.shape('WRLDTZA')
plot(Map2poly(us))
plot(Map2poly(recenter.Map(us,-100)))
plot(Map2poly(recenter.Map(tz,-165)))
## End(Not run)


Plot a map in an rgl window

Description

Plots a map (from a Map object from package maptools) on a unit sphere in an rgl window that can then be interactively rotated.

Usage

rgl.Map(Map, which, ...)

Arguments

Map | A Map object
which | Vector indicating the subset of polygons to plot.
... | Additional arguments passed on to rgl.lines.

Details

This assumes that the map is cordinates in degrees and plots the map on a unit sphere in an rgl window making a globe. You can then rotate the globe by clicking and dragging in the window.

Value

There is no return value, this function is run for its side effect.

Note

This function is still beta level software (some extra lines show up).

Author(s)

Greg Snow greg.snow@intermountainmail.org

See Also

rgl in package rgl, plot.Map in package maptools

Examples

## Not run:
# assumes that the time zone shape files have been downloaded
# from: http://openmap.bbn.com/data/shape/timezone/

tz <- read.shape('WRLDTZA')
rgl.Map(tz)
rgl.spheres(0,0,0,.999, col='darkblue')
## End(Not run)


Demonstrate ROC curves by interactively building one

Description

This demonstration allows you to interactively build a Receiver Operator Curve to better understand what goes into creating them.

Usage

roc.demo(x = rnorm(25, 10, 1), y = rnorm(25, 11, 1.5))

Arguments

x | Data values for group 1 (controls).
y | Data values for group 2 (cases).

Details

Density plots for the 2 groups will be created in the lower panel of the plot (colored red (group 1) and blue (group 2)) along with rug plots of the actual datapoints. There is also a green vertical line that represents a decision rule cutoff, any points higher than the cutoff are predicted to be in group 2 and points less than the cuttoff are predicted to be in group 1. The sensitivity and specificity for the current cuttoff value are printed below the plot.

A Tk slider box is also created that allows you to move the cuttoff value and update the plots. As the cutoff value changes, the different combinations of sensitivity and specificity are added to the ROC curve in the top panel (the point corresponding to the current cuttoff value is highlighted in red). A line is also drawn from the point representing sensitivity and specificity both equal to 1 to the point closest to it.

Value

No meaninful value is returned, this function is run solely for the side effects.

Author(s)

Greg Snow greg.snow@intermountainmail.org

See Also

slider, ROC function in package Epi, auROC in package limma, package ROC

Examples

## Not run:
roc.demo()
with(CO2,
  roc.demo(uptake[Type=='Mississippi'],
           uptake[Type=='Quebec']       )
)
## End(Not run)


Interactively rotate 3D plots

Description

Interactively rotate common 3d plots: cloud, persp, and wireframe.

Usage

rotate.cloud(formula, ...)
rotate.persp(x, y, z)
rotate.wireframe(formula, ...)

Arguments

formula | formula, see cloud and wireframe
x | x, see persp
y | y, see persp
z | z, see persp
... | additional arguments passed on to cloud or persp

Details

Use these functions just like cloud, persp, and wireframe. In addition to the default plot a Tk slider window will be created that will allow you to rotate the plot.

The rotations parameters are passed the screen argument of cloud and wireframe and the theta, phi, r, d, ltheta, lphi, and shade arguments of persp.

For cloud and wireframe plots the order of the x, y, and z argumets can be rearanged, just type the appropriate letters in the boxes on the left, then press the "refresh" button (changing the order changes the plot for these 2 plots).

Value

These functions are run for the side effects of the plots and Tk windows, nothing meaninful is returned.

Author(s)

Greg Snow greg.snow@intermountainmail.org

See Also

cloud it the lattice package, persp, wireframe in the lattice package

Examples

## Not run:
rotate.cloud(Sepal.Length ~ Petal.Length*Petal.Width, data=iris)

rotate.wireframe(volcano)

z <- 2 * volcano        # Exaggerate the relief
x <- 10 * (1:nrow(z))   # 10 meter spacing (S to N)
y <- 10 * (1:ncol(z))   # 10 meter spacing (E to W)
rotate.persp(x,y,z)
## End(Not run)


Interactively demonstrate correlations

Description

Make a scatterplot and a Tk slider window that allows you to interactively set the correlation and/or R^2.

Usage

run.cor.examp(n = 100, seed)
run.cor2.examp(n = 100, seed)

Arguments

n | Number of points to plot.
seed | What seed to ues.

Details

The function run.cor.examp draws a scatterplot and allows you to set the correlation using a Tk slider window.

The function run.cor2.examp does the same, but has a slider for R^2 as well as the correlation, when either slider is moved the other one will update to match.

Value

These functions are run for their side effect, nothing meaningful is returned.

Author(s)

Greg Snow greg.snow@intermountainmail.org

See Also

cor, slider

Examples

## Not run:
run.cor.examp()
run.cor2.examp()
## End(Not run)


Create a histogram and interactively change the number of bars.

Description

Create a histogram then use a Tk slider window to change the number of bars, the minimum, and the maximum.

Usage

run.hist.demo(x)

Arguments

x Data to plot.

Details

Draws a histogram and creates a Tk slider window that allows you to explore how changing the parameters affects the appearance of the plot.

Value

No meaninful value is returned.

Author(s)

Greg Snow greg.snow@intermountainmail.org

See Also

hist, slider

Examples

## Not run:
run.hist.demo( rnorm(250, 100, 5) )
## End(Not run)


slider / button control widgets

Description

slider constructs a Tcl/Tk-widget with sliders and buttons automated calculation and plotting. For example slider allows complete all axes rotation of objects in a plot.

Usage

slider(sl.functions, sl.names, sl.mins, sl.maxs, sl.deltas, sl.defaults, but.functions, but.names, no, set.no.value, obj.name, obj.value, reset.function, title)

Arguments

sl.functions | set of functions or function connected to the slider(s)
sl.names | labels of the sliders
sl.mins | minimum values of the sliders' ranges
sl.maxs | maximum values of the sliders' ranges
sl.deltas | change of step per click
sl.defaults | default values for the sliders
but.functions | function or list of functions that are assigned to the button(s)
but.names | labels of the buttons
no | slider(no=i) requests slider i
set.no.value | slider(set.no.value=c(i,val)) sets slider i to value val
obj.name | slider(obj.name=name) requests the value of variable name from environment slider.env
obj.value | slider(obj.name=name,obj.value=value) assigns value to variable name in environment slider.env
reset.function | function that comprises the commands of the reset.button
title | title of the control window

Details

With slider you can: a) define (multiple) sliders and buttons, b) request or set slider values, and c) request or set variables in the environment slider.env. Slider function management takes place in the environment slider.env. If slider.env is not found it is generated.

Definition
... of sliders: First of all you have to define sliders, buttons and the attributes of them. Sliders are established by six arguments: sl.functions, sl.names, sl.minima, sl.maxima,sl.deltas, and sl.defaults. The first argument, sl.functions, is either a list of functions or a single function that entails the commands for the sliders. If there are three sliders and slider 2 is moved with the mouse the function stored in sl.functions[[2]] (or in case of one function for all sliders the function sl.functions) is called.
Definition
... of buttons: Buttons are defined by a vector of labels but.names and a list of functions: but.functions. If button i is pressed the function stored in but.functions[[i]] is called.
Requesting
... a slider: slider(no=1) returns the actual value of slider 1, slider(no=2) returns the value of slider 2, etc. You are allowed to include expressions of the type slider(no=i) in functions describing the effect of sliders or buttons.
Setting
... a slider: slider(set.no.value=c(2,333)) sets slider 2 to value 333. slider(set.no.value=c(i,value)) can be included in the functions defining the effects of moving sliders or pushing buttons.
Variables
... of the environment slider.env: Sometimes information has to be trransferred back and forth between functions defining the effects of sliders and buttons. Imagine for example two sliders: one to control p and another one to control q, but they should satisfy: p+q=1. Consequently, you have to correct the value of the first slider after the second one was moved. To prevent the creation of global variables store them in the environment slider.env. Use slider(obj.name="p.save",obj.value=1-slider(no=2)) to assign value 1-slider(no=2) to the variable p.save . slider(obj.name=p.save) returns the value of variable p.save.

Value

Using slider in definition mode slider returns the value of new created the top level widget. slider(no=i) returns the actual value of slider i. slider(obj.name=name) returns the value of variable name in environment slider.env.

Note

You can move the slider in 3 different ways: You can left click and drag the slider itself, you can left click in the trough to either side of the slider and the slider will move 1 unit in the direction you clicked, or you can right click in the trough and the slider will jump to the location you clicked at.

Author(s)

Hans Peter Wolf

See Also

sliderv

Examples

# example 1, sliders only
## Not run:
## This example cannot be run by examples() but should work in an interactive R session
plot.sample.norm<-function(){
 refresh.code<-function(...){
   mu<-slider(no=1); sd<-slider(no=1); n<-slider(no=3)
   x<-rnorm(n,mu,sd)
   plot(x)
 }
 slider(refresh.code,sl.names=c("value of mu","value of sd","n number of observations"),
       sl.mins=c(-10,.01,5),sl.maxs=c(+10,50,100),sl.deltas=c(.01,.01,1),sl.defaults=c(0,1,20))
}
plot.sample.norm()
## End(Not run)

# example 2, sliders and buttons
## Not run:
## This example cannot be run by examples() but should work in an interactive R session
plot.sample.norm.2<-function(){
 refresh.code<-function(...){
   mu<-slider(no=1); sd<-slider(no=2); n<-slider(no=3)
   type=  slider(obj.name="type")
   x<-rnorm(n,mu,sd)
   plot(seq(x),x,ylim=c(-20,20),type=type)
 }
 slider(refresh.code,sl.names=c("value of mu","value of sd","n number of observations"),
       sl.mins=c(-10,.01,5),sl.maxs=c(10,10,100),sl.deltas=c(.01,.01,1),sl.defaults=c(0,1,20),
       but.functions=list(
              function(...){slider(obj.name="type",obj.value="l");refresh.code()},
              function(...){slider(obj.name="type",obj.value="p");refresh.code()},
              function(...){slider(obj.name="type",obj.value="b");refresh.code()}
       ),
       but.names=c("lines","points","both"))
  slider(obj.name="type",obj.value="l")
}
plot.sample.norm.2()
## End(Not run)

# example 3, dependent sliders
## Not run:
## This example cannot be run by examples() but should work in an interactive R session
print.of.p.and.q<-function(){
 refresh.code<-function(...){
   p.old<-slider(obj.name="p.old")
   p<-slider(no=1); if(abs(p-p.old)>0.001) {slider(set.no.value=c(2,1-p))}
   q<-slider(no=2); if(abs(q-(1-p))>0.001) {slider(set.no.value=c(1,1-q))}
   slider(obj.name="p.old",obj.value=p)
   cat("p=",p,"q=",1-p,"\n")
 }
 slider(refresh.code,sl.names=c("value of p","value of q"),
       sl.mins=c(0,0),sl.maxs=c(1,1),sl.deltas=c(.01,.01),sl.defaults=c(.2,.8))
 slider(obj.name="p.old",obj.value=slider(no=1))
}
print.of.p.and.q()
## End(Not run)

# example 4, rotating a surface
## Not run:
## This example cannot be run by examples() but should work in an interactive R session
R.veil.in.the.wind<-function(){
  # Mark Hempelmann / Peter Wolf
  par(bg="blue4", col="white", col.main="white",
      col.sub="white", font.sub=2, fg="white") # set colors and fonts
  samp  <- function(N,D) N*(1/4+D)/(1/4+D*N)
  z<-outer(seq(1, 800, by=10), seq(.0025, 0.2, .0025)^2/1.96^2, samp) # create 3d matrix
  h<-100
  z[10:70,20:25]<-z[10:70,20:25]+h; z[65:70,26:45]<-z[65:70,26:45]+h
  z[64:45,43:48]<-z[64:45,43:48]+h; z[44:39,26:45]<-z[44:39,26:45]+h
  x<-26:59; y<-11:38; zz<-outer(x,y,"+"); zz<-zz*(65<zz)*(zz<73)
  cz<-10+col(zz)[zz>0];rz<-25+row(zz)[zz>0]; z[cbind(cz,rz)]<-z[cbind(cz,rz)]+h
  refresh.code<-function(...){
    theta<-slider(no=1); phi<-slider(no=2)
    persp(x=seq(1,800,by=10),y=seq(.0025,0.2,.0025),z=z,theta=theta,phi=phi,
          scale=T, shade=.9, box=F, ltheta = 45,
          lphi = 45, col="aquamarine", border="NA",ticktype="detailed")
  }
  slider(refresh.code, c("theta", "phi"), c(0, 0),c(360, 360),c(.2, .2),c(85, 270)  )
}
R.veil.in.the.wind()
## End(Not run)

## The function is currently defined as
function(sl.functions,sl.names,sl.mins,sl.maxs,sl.deltas,sl.defaults,
                  but.functions,but.names,
                  no,set.no.value,obj.name,obj.value,
                  reset.function,title){
  # slider, version2, pw 040107
  if(!missing(no)) return(as.numeric(tclvalue(get(paste("slider",no,sep=""),env=slider.env))))
  if(!missing(set.no.value)){ try(eval(parse(text=paste("tclvalue(slider",set.no.value[1],")<-",
                            set.no.value[2],sep="")),env=slider.env)); return(set.no.value[2]) }
  if(!exists("slider.env")) slider.env<<-new.env()
  if(!missing(obj.name)){
    if(!missing(obj.value)) assign(obj.name,obj.value,env=slider.env) else
      obj.value<-get(obj.name,env=slider.env)
    return(obj.value)
  }
  if(missing(title)) title<-"slider control widget"
  require(tcltk); nt<-tktoplevel(); tkwm.title(nt,title); tkwm.geometry(nt,"+0+0")
  if(missing(sl.names)) sl.names<-NULL
  if(missing(sl.functions)) sl.functions<-function(...){}
  for(i in seq(sl.names)){
    eval(parse(text=paste("assign('slider",i,"',tclVar(sl.defaults[i]),env=slider.env)",sep="")))
    tkpack(fr<-tkframe(nt));  lab<-tklabel(fr, text=sl.names[i], width="25")
    sc<-tkscale(fr,from=sl.mins[i],to=sl.maxs[i],showvalue=T,resolution=sl.deltas[i],orient="horiz")
    tkpack(lab,sc,side="right"); assign("sc",sc,env=slider.env)
    eval(parse(text=paste("tkconfigure(sc,variable=slider",i,")",sep="")),env=slider.env)
    sl.fun<-if(length(sl.functions)>1) sl.functions[[i]] else sl.functions
    if(!is.function(sl.fun)) sl.fun<-eval(parse(text=paste("function(...){",sl.fun,"}")))
    tkconfigure(sc,command=sl.fun)
  }
  assign("slider.values.old",sl.defaults,env=slider.env)
  tkpack(f.but<-tkframe(nt),fill="x")
  tkpack(tkbutton(f.but, text="Exit", command=function()tkdestroy(nt)),side="right")
  if(missing(reset.function)) reset.function<-function(...) print("relax")
  if(!is.function(reset.function))
    reset.function<-eval(parse(text=paste("function(...){",reset.function,"}")))
  tkpack(tkbutton(f.but, text="Reset", command=function(){
         for(i in seq(sl.names))
            eval(parse(text=paste("tclvalue(slider",i,")<-",sl.defaults[i],sep="")),env=slider.env)
         reset.function()  }  ),side="right")
  if(missing(but.names)) but.names<-NULL
  for(i in seq(but.names)){
    but.fun<-if(length(but.functions)>1) but.functions[[i]] else but.functions
    if(!is.function(but.fun))but.fun<-
       eval(parse(text=paste("function(...){",but.fun,"}")))
    tkpack(tkbutton(f.but, text=but.names[i], command=but.fun),side="left")
    cat("button",i,"eingerichtet")
  }
  invisible(nt)
}


Create a Tk slider window

Description

Create a Tk slider window with the sliders positioned vertically instead of horizontally.

Usage

sliderv(refresh.code, names, minima, maxima, resolutions, starts, title = "control", no = 0, set.no.value = 0)

Arguments

refresh.code | Function to be called when sliders are moved.
names | Labels for the sliders.
minima | Vector of minimum values for the sliders.
maxima | Vector of maximum values for the sliders.
resolutions | Vector of resolutions for the sliders.
starts | Vector of starting values for the sliders.
title | Title to put at the top of the Tk box.
no | The number of the slider whose value you want.
set.no.value | Vector of length 2 with the number of slider to set and the new value.

Details

This is a variation on the slider function with vertical sliders arranged in a row rather than horizontal sliders arranged in a column.

This is based on an early version of slider and therefore does not have as many bells and whistles (but sometimes fits the screen better).

Value

Returns the value of a given slider when used as: slider(no=i).

Note

You can move the slider in 3 different ways: You can left click and drag the slider itself, you can left click in the trough to either side of the slider and the slider will move 1 unit in the direction you clicked, or you can right click in the trough and the slider will jump to the location you clicked at.

Author(s)

Greg Snow greg.snow@intermountainmail.org

See Also

slider

Examples

## Not run:
face.refresh <- function(...){
        vals <- sapply(1:15, function(x) slider(no=x))
        faces( rbind(0, vals, 1), scale=F)
}

sliderv( face.refresh, as.character(1:15), rep(0,15), rep(1,15),
        rep(0.05, 15), rep(0.5,15), title='Face Demo')
## End(Not run)


Squish the plotting area to a specified aspect ratio

Description

Adjusts the plotting area to a specific aspect ratio. This is different from using the asp argument in that it puts the extra space in the margins rather than inside the plotting region.

Usage

squishplot(xlim, ylim, asp = 1)

Arguments

xlim | The x limits of the plot, or the entire x vector.
ylim | The y limits of the plot, or the entire y vector.
asp | The y/x aspect ratio.

Details

This function sets the plot area of the current graph device so that the following plot command will plot with the specified aspect ratio.

This is different from using the asp argument to plot.default in where the created white space goes (see the example). Using plot.default will place the whitespace within the plotting region and can result in the axes and annotations being quite far from the actual data. This command sets up the plotting region so that the extra whitespace is in the margin areas and moves the axes and annotations close to the data.

Any other desired parameter settings or resizing of the graphics device should be set before calling squishplot, especially settings dealing with multiple figures or margin areas.

After plotting, the parameters need to be reset or later plots may come out wrong.

Value

Invisible list containing the 'plt' values from par that were in place before the call to squishplot that can be used to reset the graphical parameters after plotting is finished.

Note

Remember to set other graphical parameters, then call squishplot, then call the plotting function(s), then reset the parameters.

Author(s)

Greg Snow greg.snow@intermountainmail.org

See Also

plot.default, plot.window, par

Examples

x <- rnorm(25, 10, 2 )
y <- 5 + 1.5*x + rnorm(25,0,2)

par(mfrow=c(1,3))
plot(x,y)

op <- squishplot(x,y,1)
plot(x,y)
par(op)

plot(x,y, asp=1)


Map object for the US State Visibility Based Map

Description

A Map object (for use with the maptools package) to plot a Visibility Based Map.

Usage

data(state.vbm)

Details

A polylist object (for use with the maptools package) to plot a map of the US states where the sizes of the states have been adjusted to be more equal.

This map can be useful for plotting state data using colors patterns without the larger states dominating and the smallest states being lost.

The original map is copyrighted by Mark Monmonier. Official publications based on this map should acknowledge him. Comercial publications of maps based on this probably need permission from him to use.

Source

The data was converted from the maps library for S-PLUS. S-PLUS uses the map with permission from the author. This version of the data has not received permission from the author (no attempt made, not that it was refused), most of my uses I feel fall under fair use and do not violate copyright, but you will need to decide for yourself and your applications.

References

http://www.markmonmonier.com/index.htm, http://www.math.yorku.ca/SCS/Gallery/bright-ideas.html

Examples

library(maptools)
data('state.vbm',package='TeachingDemos')
plot(state.vbm,axes=FALSE)

tmp <- state.x77[,'HS Grad']
tmp2 <- cut(tmp, seq(min(tmp),max(tmp), length.out=11),
  include.lowest=TRUE)
plot(state.vbm,col=cm.colors(10)[tmp2],axes=FALSE)


strip function for trellis plots

Description

The strip.shingle function is a modification of strip.default that shows the range of each shingle in the panel strips as well as the name of the shingle variable.

This function is now depricated, you can do the same thing with strip.custom now, see the example below.

Usage

strip.shingle(which.given, which.panel, var.name, factor.levels, shingle.intervals, strip.names = c(FALSE, TRUE), style = 1, bg = trellis.par.get("strip.background")$col[which.given], fg = trellis.par.get("strip.shingle")$col[which.given], par.strip.text = trellis.par.get("add.text"),...)

Arguments

which.given | See strip.default.
which.panel | See strip.default.
var.name | See strip.default.
factor.levels | See strip.default.
shingle.intervals | See strip.default.
strip.names | See strip.default.
style | See strip.default.
bg | See strip.default.
fg | See strip.default.
par.strip.text | See strip.default.
... | Arguments passed through

Details

Does the same as strip.default with the exception that by default when a shingle is used as the conditioning variable, the strip at the top of the panel will inclued the name of the conditioning variable and the range of the shingle for that plot.

This function is now depricated, see the example below for the way to do the same thing using core lattice code.

Value

This function is only usefull for its side effects.

Author(s)

Greg Snow greg.snow@intermountainmail.org

See Also

strip.default, strip.custom

Examples

library(lattice)
Depth <- equal.count(quakes$depth, number=8, overlap=.1)

## Not run:
# this is the old way, don't do this anymore use the example below
xyplot(lat ~ long | Depth, data = quakes, strip=strip.shingle)
## End(Not run)

## Not run:
xyplot(lat ~ long | Depth, data=quakes,
  strip=strip.custom(strip.names=TRUE, strip.levels=TRUE),
  par.strip.text=list(cex=0.75))
## End(Not run)


Embed a new plot within an existing plot

Description

Subplot will embed a new plot within an existing plot at the coordinates specified (in user units of the existing plot).

Usage

subplot(fun, x, y, size=c(1,1), vadj=0.5, hadj=0.5, pars=NULL)

Arguments

fun | an expression defining the new plot to be embedded.
x | x-coordinate(s) of the new plot (in user coordinates of the existing plot).
y | y-coordinate(s) of the new plot, x and y can be specified in any of the ways understood by xy.coords.
size | The size of the embedded plot in inches if x and y have length 1.
vadj | vertical adjustment of the plot when y is a scalar, the default is to center vertically, 0 means place the bottom of the plot at y, 1 places the top of the plot at y.
hadj | horizontal adjustment of the plot when x is a scalar, the default is to center horizontally, 0 means place the left edge of the plot at x, and 1 means place the right edge of the plot at x.
pars | a list of parameters to be passed to par before running fun.

Details

The coordinates x and y can be scalars or vectors of length 2. If vectors of length 2 then they determine the opposite corners of the rectangle for the embedded plot (and the parameters size, vadj, and hadj are all ignored.

If x and y are given as scalars then the plot position relative to the point and the size of the plot will be determined by the arguments size, vadj, and hadj. The default is to center a 1 inch by 1 inch plot at x,y. Setting vadj and hadj to (0,0) will position the lower left corner of the plot at (x,y).

The rectangle defined by x, y, size, vadj, and hadj will be used as the plotting area of the new plot. Any tick marks, axis labels, main and sub titles will be outside of this rectangle.

Any graphical parameter settings that you would like to be in place before fun is evaluated can be specified in the pars argument (warning: specifying layout parameters here (plt, mfrow, etc.) may cause unexpected results).

After the function completes the graphical parameters will have been reset to what they were before calling the function (so you can continue to augment the original plot).

Value

An invisible list with the graphical parameters that were in effect when the subplot was created. Passing this list to par will enable you to augment the embedded plot.

Author(s)

Greg Snow greg.snow@intermountainmail.org

See Also

cnvrt.coords, par, symbols

Examples

# make an original plot
plot( 11:20, sample(51:60) )

# add some histograms

subplot( hist(rnorm(100)), 15, 55)
subplot( hist(runif(100),main='',xlab='',ylab=''), 11, 51, hadj=0, vadj=0)
subplot( hist(rexp(100, 1/3)), 20, 60, hadj=1, vadj=1, size=c(0.5,2) )
subplot( hist(rt(100,3)), c(12,16), c(57,59), pars=list(lwd=3,ask=FALSE) )

# augment a map
if( require(TeachingDemos) && require(maptools) ){
        plot(state.vbm,fg=NULL)
        tmp <- t(sapply(state.vbm, function(x) attr(x,"centroid")))
        for( i in 1:50 ){
                tmp2 <- as.matrix(USArrests[i,c(1,4)])
                tmp3 <- max(USArrests[,c(1,4)])
                subplot( barplot(tmp2, ylim=c(0,tmp3),names=c('',''),yaxt='n'),
                        x=tmp[i,1], y=tmp[i,2], size=c(.1,.1))
        }
}

tmp <- rnorm(25)
qqnorm(tmp)
qqline(tmp)
tmp2 <- subplot( hist(tmp,xlab='',ylab='',main=''),
                cnvrt.coords(0.1,0.9,'plt')$usr, vadj=1, hadj=0 )
abline(v=0, col='red') # wrong way to add a reference line to histogram

# right way to add a reference line to histogram
op <- par(no.readonly=TRUE)
par(tmp2)
abline(v=0, col='green')
par(op)

# scatter-plot using images
if(require(rimage)){
        data(logo)

        x <- runif(10)
        y <- runif(10)

        plot(x,y,type='n')
        for(i in 1:10){
                subplot(plot(logo),x[i],y[i],size=c(.3,.3))
        }
}


Change the Color and Styles of points interactively

Description

Creates a Tk window with a scatterplot matrix, then allows you to "brush" the points to change their color and/or style.

Usage

tkBrush(mat,hscale=1.75,vscale=1.75,wait=TRUE,...)

Arguments

mat | A matrix of the data to plot, columns are variables, rows are observations, same as pairs
hscale | Passed to tkrplot
vscale | Passed to tkrplot
wait | Should the function wait for you to finish, see below
... | Additional arguments passed to the panel functions

Details

This function creates a Tk window with a pairs plot of mat, then allows you to interactively move a rectangle (the brush) over the points to change their color and plotting character.

The arrow keys can be used to change the size and shape of the brush. The left arrow makes the rectangle wider, the right makes it narrower. The up arrow key makes it taller, the right makes it shorter.

When the mouse button is not pressed the points inside the brush will change while in the brush, but return to their previous state when the brush moves off them. If the mouse button is pressed then the points inside the brush will be changed and the change will remain until a different set of conditions is brushed on them.

The style of the brushed points is determined by the values of the 2 entry boxes on the right side of the plot. You can specify the plotting character in the pch box, this can be anything that you would regularly pass to the pch argument of points, e.g. an integer or single character. You can specify the color of the brushed points using the color entry box, specify the name of any color recognized by R (see colors), if this box does not contain a legal color name then black will be used.

If wait is FALSE then the Tk window will exist independently of R and you can continue to do other things in the R window, in this case the function returns NULL. If wait is TRUE then R waits for you to close the Tk window (using the quit button) then returns a list with the colors and plotting characters resulting from your brushing, this information can be used to recreate the plot using pairs on a new graphics device (for printing or saving).

Value

Either NULL (if Wait=FALSE) or a list with components col and pch corresponding to the state of the points.

Author(s)

Greg Snow greg.snow@intermountainmail.org

See Also

pairs,colors,points, the iplots package

Examples

## Not run:

# Iris dataset

out1 <- tkBrush(iris)

#  Now brush the points

pairs(iris, col=out1$col, pch=out1$pch)

# or

colhist <- function(x,...){
    tmp <- hist(x,plot=F)
    br <- tmp$breaks
    w <- as.numeric(cut(x,br,include.lowest=TRUE))
    sy <- unlist(lapply(tmp$counts,function(x)seq(length=x)))
    my <- max(sy)
    sy <- sy/my
    my <- 1/my
    sy <- sy[order(order(x))]
    tmp.usr <- par('usr'); on.exit(par(usr=tmp.usr))
    par(usr=c(tmp.usr[1:2],0,1.5))
    rect(br[w], sy-my, br[w+1], sy,
       col=out1$col, # note out1$col is hardcoded here.
       border=NA)
    rect(br[-length(br)], 0, br[-1], tmp$counts*my)
}
pairs(iris, col=out1$col, pch=out1$pch, diag.panel=colhist)

# some spheres

s1 <- matrix(nrow=0,ncol=3)

while( nrow(s1) < 1000 ){
        tmp <- rnorm(3)
        if( sum(tmp^2) <= 1 ){
                s1 <- rbind(s1,tmp)
        }
}

s2 <- matrix(rnorm(3000), ncol=3)
s2 <- s2/apply(s2,1,function(x) sqrt(sum(x^2)))

tkBrush(s1, wait=FALSE)
tkBrush(s2, wait=FALSE)

# now paint values where var 2 is close to 0 in both plots
# and compare the var 1 and var 3 relationship

## End(Not run)


Interactively demonstrate regression trees

Description

Interactively recursively partition a dataset to demonstrate regression trees.

Usage

tree.demo(x, y)

Arguments

x | The predictor variable.
y | The response variable.

Details

This function first creates a scatterplot of x and y and shows the residual sum of squares from fitting a horizontal line to the y-values.

Clicking anywhere on the graph will show an updated graph where the data is partitioned into 2 groups based on the x-value where you clicked with a horizontal line fit to each group (including showing the updated residual sum of squares). Clicking again will move the partitioning value based on the new click.

When you have found a good partitioning (reduces the RSS), right click and choose 'stop' and that partition will become fixed. Now you can click to do a second set of partions (breaking the data into 3 groups).

To finish the demo, right click and choose 'stop', then right click again and choose 'stop' again.

Value

A vector with the x-values of the cut points that you selected (sorted).

Author(s)

Greg Snow greg.snow@intermountainmail.org

See Also

The rpart and tree packages

Examples

## Not run:
data('ethanol', package='lattice')
print(with(ethanol, tree.demo(E,NOx)))
## End(Not run)


Create or add to a Trilinear Plot

Description

Create (or add to) a trilinear plot of 3 proportions that sum to 1.

Usage

triplot(x, y = NULL, z = NULL, labels = dimnames(x)[[2]], txt =
dimnames(x)[[1]], legend = NULL, legend.split = NULL, inner = TRUE, inner.col = c("lightblue", "pink"), inner.lty = c(2, 3), add = FALSE, main = "", ...)

Arguments

x | Vector or matrix of up to 3 columns.
y | Vector (if x is a vector).
z | Vector (if x is a vector).
labels | Labels for the 3 components (printed at corners).
txt | Vector of text strings to be plotted instead of points.
legend | Labels for the data points
legend.split | What proportion of the labels will go on the left.
inner | Logical, should the inner reference lines be plotted.
inner.col | Colors for the 2 inner triangles.
inner.lty | Line types for the 2 inner triangles.
add | Add points to existing plot (TRUE), or create a new plot (FALSE).
main | Main title for the plot.
... | Additional arguments passed on to points or text.

Details

Trilinear plots are useful for visualizing membership in 3 groups by plotting sets of 3 proportions that sum to 1 within each set.

The data can be passed to the function as a matrix with either 2 or 3 columns, or as seperate vectors to x, y, and optionaly z. If 2 columns are passed in, then they must be between 0 and 1 and the 3rd column will be created by subtracting both from 1. If 3 columns of data are given to the function then each will be divided by the sum of the 3 columns (they don't need to sum to 1 before being passed in).

Value

An invisible matrix with 2 columns and the same number of rows as x corresponding to the points plotted (after transforming to 2 dimensions).
The return matrix can be passed to identify for labeling of individual points.
Using type='n' and add=FALSE will return the transformed points without doing any plotting.

Author(s)

Greg Snow greg.snow@intermountainmail.org

References

Allen, Terry. Using and Interpreting the Trilinear Plot. Chance. 15 (Summer 2002).

See Also

triangle.plot in package ade4, ternaryplot in package vcd, tri in package cwhtool, soil.texture and triax.plot in package plotrix.

Examples

triplot(USArrests[c(1,4,2)])
tmp <- triplot(USArrests[c(1,4,2)],txt=NULL)
## Not run:  identify(tmp, lab=rownames(USArrests) )

tmp <- rbind( HairEyeColor[,,'Male'], HairEyeColor[,,'Female'])
tmp[,3] <- tmp[,3] + tmp[,4]
tmp <- tmp[,1:3]
triplot(tmp, legend=rep(c('Male','Femal'),each=4),
 col=rep(c('black','brown','red','yellow'),2))


Plot various distributions then interactivly adjust the parameters.

Description

Plot a curve of a distribution, then using a Tk slider window adjust the parameters and see how the distribution changes. Optionally also plots reference distributions.

Usage

vis.binom()
vis.gamma()
vis.normal()
vis.t()

Details

These functions plot a distribution, then create a Tk slider box that allows you to adjust the parameters of the distribution to see how the curve changes.

Check boxes are available in some cases to also show reference distributions (normal and poisson for the binomial, exponential and chi-squared for gamma, and normal for t). The exponential and chi-squared distributions are those with the same mean as the plotted gamma.

If you change the plotting ranges then you need to click on the 'refresh' button to update the plot.

Value

These functions are run for their side effects and do not return anything meaningful.

Author(s)

Greg Snow greg.snow@intermountainmail.org

See Also

dnorm, dgamma, etc.

Examples

## Not run:
vis.binom()
vis.normal()
vis.gamma()
vis.t()
## End(Not run)


Interactively visualize Box-Cox transformations

Description

Explore the Box-Cox family of distributions by plotting data transformed and untransformed and interactively choose values for lambda.

Usage

vis.boxcox(lambda = sample(c(-1, -0.5, 0, 1/3, 1/2, 1, 2), 1))

vis.boxcoxu(lambda = sample(c(-1, -0.5, 0, 1/3, 1/2, 1, 2), 1))

Arguments

lambda The true value of lambda to use.

Details

These functions will generate a sample of data and plot the untrasformed data (left panels) and the transformed data (right panels). Initially the value of lambda is 1 and the 2 sets of plots will be identical.

You then adjust the transformation parameter lambda to see how the right panels change.

The function vis.boxcox shows the effect of transforming the y-variable in a simple linear regression.

The function vis.boxcoxu shows a single variable compared to the normal distribution.

Value

These functions are run for the side effects. There is no meaningful return value.

Author(s)

Greg Snow greg.snow@intermountainmail.org

References

GEP Box; DR Cox. An Analysis of Transformations. Journal of the Royal Statitical Society. Series B, Vol. 26, No. 2 (1964) 211-252

See Also

bct, boxcox in package MASS

Examples

## Not run:
vis.boxcoxu()
vis.boxcox()
## End(Not run)


Z test for known population standard deviation

Description

Compute the test of hypothesis and compute confidence interval on the mean of a population when the standard deviation of the population is known.

Usage

z.test(x, mu = 0, stdev, alternative = c("two.sided", "less", "greater"), sd = stdev, conf.level = 0.95, ...)

Arguments

x | Vector of data values.
mu | Hypothesized mean of the population.
stdev | Known standard deviation of the population.
alternative | Direction of the alternative hypothesis.
sd | Alternative to stdev
conf.level | Confidence level for the interval computation.
... | Additional arguments are silently ignored.

Details

Most introductory statistical texts introduce inference by using the Z test and Z based confidence intervals based on knowing the population standard deviation. Most statistical packages do not include functions to do Z tests since the T test is usually more appropriate for real world situations. This function is meant to be used during that short period of learning when the student is learning about inference using Z procedures, but has not learned the T based procedures yet. Once the student has learned about the T distribution the t.test function should be used instead of this one (but the syntax is very similar, so this function should be an appropriate introductory step to learning t.test).

Value

An object of class htest containing the results

Note

This function should be used for learning only, real data should generally use t.test.

Author(s)

Greg Snow greg.snow@intermountainmail.org

See Also

t.test, print.htest

Examples

x <- rnorm(25, 100, 5)
z.test(x, 99, 5)