Package 'TeachingDemos'

Description

This package provides various demonstrations that can be used in classes or by individuals to better learn statistical concepts and usage of R. Various utility functions are also included

Details

Demonstration functions in this package include:

Utility functions include:

Author(s)

Greg Snow [email protected]

See Also

The tkrplot package

Examples

ci.examp()

clt.examp()
clt.examp(5)

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

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

plot(1:10, 1:10)
my.symbols( 1:10, 1:10, ms.polygram, n=1:10, inches=0.3 )

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))

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

Description

Comparison operators that can be chained together into something like 0 %<% x %<% 1 instead of 0 < x && x < 1.

Usage

x %<% y
x %<=% y

Arguments

Details

These functions/operators allow chained inequalities. To specify that you want the values between two values (say 0 and 1) you can use 0 %<% x %<% 1 rather than 0 < x \&\& x < 1 .

Value

A logical vector is returned that can be used for subsetting like <, but the original values are included as attributes to be used in additional comparisons.

Note

This operator is not fully associative and has different precedence than < and <=, so be careful with parentheses. See the examples.

Author(s)

Greg Snow, [email protected]

Examples

x <- -3:3

 -2 %<% x %<% 2
c( -2 %<% x %<% 2 )
x[ -2 %<% x %<% 2 ]
x[ -2 %<=% x %<=% 2 ]


x <- rnorm(100)
y <- rnorm(100)

x[ -1 %<% x %<% 1 ]
range( x[ -1 %<% x %<% 1 ] )


cbind(x,y)[ -1 %<% x %<% y %<% 1, ]
cbind(x,y)[ (-1 %<% x) %<% (y %<% 1), ]
cbind(x,y)[ ((-1 %<% x) %<% y) %<% 1, ]
cbind(x,y)[ -1 %<% (x %<% (y %<% 1)), ]
cbind(x,y)[ -1 %<% (x %<% y) %<% 1, ] # oops


3 
3

Description

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

Usage

bct(y, lambda)

Arguments

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 [email protected]

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))

Description

Plot a calendar of the specified year or month. Monthly calendars can have additional information (text/plots) added to the individual cells.

Usage

cal(month, year)

Arguments

Details

This function plots on the current (or default) graphics device a yearly or monthly calendar. It tries to guess what you want, if both year and month are ommitted then it will plot the current month. If month is an integer greater than 12 and no year is specified then that value will be used as the year for a yearly calendar. The month can be either an integer from 1 to 12 or a character string that will be matched against month.name using pmatch.

Each day of the monthly calendar is a plotting frame that can be added to using stardard low level functions, the coordinates of the plotting region (the entire box) are from 0 to 1 in both dimensions. The updateusr function can be used to change the coordinates. The return from the function (when creating a monthly calendar) can be used to select the day.

Value

Nothing is returned when a whole year calendar is created. When the month calendar is created a function is returned invisibly that if passed an integer corresponding to a day of the month will set the graphics parameters so the corresponding day in the calendar becomes the current plotting figure. See the examples below.

Author(s)

Greg Snow, [email protected]

See Also

Sys.time, as.POSIXlt, par, updateusr

Examples

cal(2011)

cal('May')

setday <- cal(11, 2011)

setday(3)
text(0.5,0.5, 'Some\nCentered\nText')

setday(8)
text(1,1,'Top Right',adj=c(1,1))

setday(18)
text(0,0,'Bottom Left', adj=c(0,0) )

setday(21)
tmp.x <- runif(25)
tmp.y <- rnorm(25, tmp.x, .1)
mrgn.x <- 0.04*diff(range(tmp.x))
mrgn.y <- 0.04*diff(range(tmp.y))
updateusr( 0:1, 0:1, range(tmp.x)+c(-1,1)*mrgn.x, range(tmp.y)+c(-1,1)*mrgn.y)
points(tmp.x, tmp.y)

setday(30)
tmp <- hist(rnorm(100), plot=FALSE)
updateusr( 0:1, 0:1, range(tmp$breaks), range(tmp$counts*1.1,0) )
lines(tmp)

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

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 [email protected]

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

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

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:

Author(s)

Greg Snow, [email protected]

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])

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

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 [email protected]

See Also

z.test, t.test

Examples

ci.examp()

if(interactive()) {
  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)

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

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, it is called with a fully transparent color so if your graphics device honors transparency then you will probably see no effect.

Value

Nothing meaningful is returned

Note

This function abuses some of the intent of what par(plt=...) is supposed to mean. In R2.7.0 and beyond there is a new funcntion clip with the intended purpose of doing this in a more proper manner (however as of my last test it is not working perfectly either, so clipplot will remain undepricated for now).

It uses some hacks to make sure that the clipping region is set, but it does this by plotting some tranparent boxes, therefore you should not use this on devices where tranparency is not supported (or you may see extra boxes).

Author(s)

Greg Snow [email protected]

See Also

par, lines, clip in R2.7.0 and later

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)

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, norm.param=list(mean=0,sd=1),
          gamma.param=list(shape=1, rate=1/3), unif.param=list(min=0,max=1),
          beta.param=list(shape1=0.35, shape2=0.25))

Arguments

Details

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

The norm.param, gamma.param, unif.param, and beta.param arguments can be used to change the parameters of the generating distributions.

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 [email protected]

See Also

rnorm, rexp, runif, rbeta

Examples

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

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

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:

Note

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

This function is now depricated with the new functions grconvertX and grconvertY in R version 2.7.0 and beyond. These new functions use the correct coordinate system names and have more coordinate systems available, you should start using them instead.

Author(s)

Greg Snow [email protected]

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, grconvertX and grconvertY in R2.7.0 and later

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)

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)

Description

Convert colors to grey/grayscale so that you can see how your plot will look after photocopying or printing to a non-color printer.

Usage

col2grey(cols)
col2gray(cols)

Arguments

Details

converts colors to greyscale using the formula grey=0.3*red + 0.59*green + 0.11*blue. This allows you to see how your color plot will approximately look when printed on a non-color printer or photocopied.

Value

A vector of colors (greys) corresponding to the input colors.

Author(s)

Greg Snow [email protected]

See Also

grey, col2rgb, dichromat package

Examples

par(mfcol=c(2,2))
tmp <- 1:3
names(tmp) <- c('red','green','blue')

barplot( tmp, col=c('red','green','blue') )
barplot( tmp, col=col2gray( c('red','green','blue') ) )

barplot( tmp, col=c('red','#008100','#3636ff') )
barplot( tmp, col=col2grey( c('red','#008100','#3636ff') ) )

Description

This function creates a scatterplot of the data, then adds colored rectangles between the points and the mean of x and y to represent the idea of the correlation coefficient.

Usage

cor.rect.plot(x, y, corr = TRUE, xlab = deparse(substitute(x)),
  ylab = deparse(substitute(y)), col = c("#ff000055", "#0000ff55"),
  ...)

Arguments

Details

This will create a scatterplot of the data, draw refrence lines at the mean of x and the mean of y, then draw rectangles from the mean point to the data points. The right and top axes will show the centered (and possibly scaled if corr=TRUE) values.

The idea is that the correlation/covariance is based on summing the area of the "positive" rectangles and subtracting the sum of the areas of the "negative" rectangles (then dividing by n-1). If the positive and negative areas are about the same then the correlation/covariance is near 0, if there is more area in the positive rectangles then the correlation/covariance will be positive.

Value

This function returns an invisible NULL, it is run for its side effects.

Author(s)

Greg Snow, [email protected]

See Also

cor

Examples

## low correlation
x <- rnorm(25)
y <- rnorm(25)
cor(x,y)
cor.rect.plot(x,y)

## Positive correlation
x <- rnorm(25)
y <- x + rnorm(25,3, .5)
cor(x,y)
cor.rect.plot(x,y)

## negative correlation
x <- rnorm(25)
y <- rnorm(25,10,1.5) - x
cor(x,y)
cor.rect.plot(x,y)

## zero correlation but a definite relationship
x <- -5:5
y <- x^2
cor(x,y)
cor.rect.plot(x,y)

Description

Simulate and optionally plot rolls of dice.

Usage

dice(rolls = 1, ndice = 2, sides = 6, plot.it = FALSE, load = rep(1, sides))
## S3 method for class 'dice'
plot(x, ...)

Arguments

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 [email protected]

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 ) )
colMeans(tmp)

# Efron's dice
ed <- list( rep( c(4,0), c(4,2) ),
  rep(3,6), rep( c(6,2), c(2,4) ),
  rep( c(5,1), c(3,3) ) )

tmp <- dice( 10000, ndice=4 )
ed.out <- sapply(1:4, function(i) ed[[i]][ tmp[[i]] ] )

mean(ed.out[,1] > ed.out[,2])
mean(ed.out[,2] > ed.out[,3])
mean(ed.out[,3] > ed.out[,4])
mean(ed.out[,4] > ed.out[,1])


## 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) )

Description

Takes an integer or vector of integers and returns a vector, list, or matrix of the individual digits (decimal) that make up that number.

Usage

digits(x, n = NULL, simplify = FALSE)

Arguments

Details

This function transforms an integer (or real ignoring the fractional part) into the decimal digits that make of the decimal representation of the number using modular mathematics rather than converting to character, splitting the string, and converting back to numeric.

Value

If x is of length 1 then a vector of the digits is returned.

If x is a vector and simplify is FALSE then a list of vectors is returned, one element for each element of x.

If x is a vector and simplify is TRUE then a matrix with 1 column for each element of x.

Author(s)

Greg Snow [email protected]

See Also

%%, %/%, strsplit

Examples

digits( 12345 )
digits( 567, n=5 )

x <- c(1, 23, 456, 7890)

digits(x)
digits(x, simplify=TRUE)

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

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 [email protected]

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)) )

Description

These functions create a scatterplot of your points and place labels for the points on them. You can then use the mouse to click and drag the labels to new positions with a line stretching between the point and label.

Usage

dynIdentify(x, y, labels = seq_along(x),
  corners = cbind(c(-1, 0, 1, -1, 1, -1, 0, 1),
    c(1, 1, 1, 0, 0, -1, -1, -1)), ...)
TkIdentify(x, y, labels=seq_along(x), hscale=1.75, vscale=1.75,
  corners = cbind( c(-1,0,1,-1,1,-1,0,1), c(1,1,1,0,0,-1,-1,-1) ),...)

Arguments

Details

These functions create a scatterplot of the x and y points with the labels (from the argument above) plotted on top. You can then use the mouse to click and drag the labels to new locations. The Tk version shows the labels being dragged, dynIdentify does not show the labels being dragged, but the label will jump to the new location as soon as you release the mouse button.

The corners argument is a 2 column matrix that gives the allowable points at which the line from the point can attach to the label (so the line does not cover thelabel). The first column represents the x-coordinates and the 2nd column the y-coordinates. A 1 represents the right/top of the label, A -1 is the left/bottom and a 0 is the center. The default values allow attachments at the 4 corners and the centers of the 4 sides of the rectangle bounding the label.

Value

A list of lists with the coordinates of the final positions of the labels and the line ends.

Note

The dynIdentify function only works on windows, TkIdentify should work on any platform with tcltk.

Author(s)

Greg Snow, [email protected]

See Also

identify

Examples

if(interactive()) {
  tmp <- TkIdentify(state.x77[,'Frost'], state.x77[,'Murder'],
  state.abb)
     ### now move the labels

     ### recreate the graph on the current device
  plot( state.x77[,'Frost'], state.x77[,'Murder'],
        xlab='Frost', ylab='Frost')
  text( tmp$labels$x, tmp$labels$y, state.abb )
  segments( state.x77[,'Frost'], state.x77[,'Murder'],
            tmp$lineends$x, tmp$lineends$y )
}

Description

Data from 46 consecutive days on weather variables used to estimate amount of evaporation from the soil.

Usage

data(evap)

Format

A data frame with 46 observations on the following 14 variables.

Obs

Observation number

Month

Month (6-June, 7-July)

day

Day of the month

MaxST

Maximum Soil Temperature

MinST

Minimum Soil Temperature

AvST

Average (integrated) Soil Temperature

MaxAT

Maximum Air Temperature

MinAT

Minimum Air Temperature

AvAT

Average (integrated) Air Temperature

MaxH

Maximum Relative Humidity

MinH

Minimum Relative Humidity

AvH

Average (integrated) Relative Humidity

Wind

Total Wind

Evap

Total evoporation from the soil

Details

The idea of the data is to predict the amount of evaporation given the other variables. Note that the "average" values are scaled differently from the others, this is more an area under the curve measure representing the total/average value.

This dataset was entered by hand from a low quality copy of the paper. If you find any typos, please e-mail them to the package maintainer.

Source

Freund, R.J. (1979) Multicollinearity etc., Some "New" Examples. Proceedings of the Statistical Computing Section, *4*, 111-112.

Examples

data(evap)
pairs(evap[,-c(1,2,3)], panel=panel.smooth)
## maybe str(evap) ; plot(evap) ...

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

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.uni-bielefeld.de/fakultaeten/wirtschaftswissenschaften : S/R - functions : faces

See Also

—

Examples

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")

if(interactive()){
  tke1 <- rep( list(list('slider',from=0,to=1,init=0.5,resolution=0.1)), 15)
  names(tke1) <- c('FaceHeight','FaceWidth','FaceShape','MouthHeight',
	'MouthWidth','SmileCurve','EyesHeight','EyesWidth','HairHeight',
	'HairWidth','HairStyle','NoseHeight','NoseWidth','EarWidth','EarHeight')
  tkfun1 <- function(...){
	tmpmat <- rbind(Min=0,Adjust=unlist(list(...)),Max=1)
	faces(tmpmat, scale=FALSE)
  }

  tkexamp( tkfun1, list(tke1), plotloc='left', hscale=2, vscale=2 )
}

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

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 found on the internet, 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 [email protected]

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')

if(interactive()){
  tke2 <- rep( list(list('slider',from=0,to=1,init=0.5,resolution=0.1)), 18)
  names(tke2) <- c('CenterWidth','TopBottomWidth','FaceHeight','TopWidth',
	'BottomWidth','NoseLength','MouthHeight','MouthCurve','MouthWidth',
	'EyesHeight','EyesBetween','EyeAngle','EyeShape','EyeSize','EyeballPos',
	'EyebrowHeight','EyebrowAngle','EyebrowWidth')
  tkfun2 <- function(...){
	tmpmat <- rbind(Min=0,Adjust=unlist(list(...)),Max=1)
	faces2(tmpmat, scale='none')
  }

  tkexamp( tkfun2, list(tke2), plotloc='left', hscale=2, vscale=2 )
}

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(hscale=1.5, vscale=1.5, wait=FALSE)
plotFagan2(hscale=1.5, vscale=1.5, wait=FALSE)
plotFagan.old()
plotFagan2.old()

Arguments

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

The old functions are run for their side effects and do not return a meaningful value. If wait is FALSE then NULL is returned, if wait is TRUE, then a list with the current values is returned.

Author(s)

Guazzetti Stefano and Greg Snow [email protected]

References

Fagan TJ. Nomogram for Bayes theorem. N Engl J Med 1975;293(5):257-61. https://ebm.bmj.com/content/6/6/164.full

See Also

slider

Examples

fagan.plot(0.8, 2)

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

if(interactive()) {
 plotFagan()

 plotFagan2()
}

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=gpenv$gp)
gp.send(cmd='replot',pipe=gpenv$gp)
gp.plot(x,y,type='p',add=FALSE, title=deparse(substitute(y)),pipe=gpenv$gp)
gp.splot(x,y,z, add=FALSE, title=deparse(substitute(z)), pipe=gpenv$gp,
  datafile=tempfile())

Arguments

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 [email protected]

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')

tmp <- expand.grid(x=1:10, y=1:10)
tmp <- transform(tmp, z=(x-5)*(y-3))
gp.splot(tmp$x, tmp$y, tmp$z)

gp.close()
 
## End(Not run)

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, lower, upper)

Arguments

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 [email protected]

See Also

hdr in the hdrcde package.

Examples

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

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

Description

These functions create a scatterplot then you Hover the mouse pointer over a point in the plot and it will show an id label for that point.

Usage

HWidentify(x, y, label = seq_along(x), lab.col="darkgreen",
pt.col="red", adj=c(0,0), clean=TRUE, xlab = deparse(substitute(x)),
ylab = deparse(substitute(y)), ...)
HTKidentify(x, y, label = seq_along(x), lab.col="darkgreen",
pt.col="red", adj=c(0,0), xlab = deparse(substitute(x)),
ylab = deparse(substitute(y)), ...)

Arguments

Details

This is an alternative to the identify function. The label only shows up for the point currently closest to the mouse pointer. When the mouse pointer moves closer to a different point, the label changes to the one for the new point. The currently labeled point is also highlighted. HWidentify only works on windows, HTKidentify requires the tkrplot package.

Value

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

Author(s)

Greg Snow, [email protected]

See Also

identify

Examples

if( interactive() ){
  tmpx <- runif(25)
  tmpy <- rnorm(25)
  HTKidentify(tmpx,tmpy, LETTERS[1:25], pch=letters)
}

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

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 [email protected]

See Also

xyplot in lattice package

Examples

if(interactive()){
require(stats)
lattice.demo(quakes$long, quakes$lat, quakes$depth)
}

Description

Data on the Growth of The Church of Jesus Christ of Latter-day Saints (commonly known as the Mormon church (https://www.churchofjesuschrist.org/comeuntochrist)).

Usage

data(ldsgrowth)

Format

A data frame with 179 observations on the following 6 variables.

Year

Year from 1830 to 2008

Members

Total number of Members

Wards

Number of Wards and Branches (individual congregations)

Stakes

Number of Stakes (a group of wards/branches)

Missions

Number of Missions

Missionaries

Number of Missionaries called

Details

The data comes from the church records and are as of December 31st of each year.

The church was officially organized on 6 April 1830 (hence the starting year of 1830).

The Missionaries column represents the number of missionaries called each year. Missionaries generally serve for about 2 years.

Source

Deseret News 2010 Church News Almanac

Examples

data(ldsgrowth)
with(ldsgrowth, plot(Year, log(Members)))

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 = FALSE,
  xlim = numeric(0), ylim = numeric(0), verbose = FALSE)

Arguments

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 [email protected]

See Also

loess, locator

Examples

if(interactive()){
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()
}

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

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 [email protected]

See Also

fitdistr in package MASS, mle in package stats4, slider

Examples

if(interactive()){
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
}

Description

These functions/data matricies are examples of what can be passed as the symb argument in the my.symbols function. They are provided both to be used for some common symbols and as examples of what can be passed as the symb argument.

Usage

ms.polygram(n, r=1, adj=pi/2, ...)
ms.polygon(n, r=1, adj=pi/2, ...)
ms.filled.polygon(n, r=1, adj=pi/2, fg=par('fg'), bg=par('fg'), ... )
ms.male
ms.female
ms.arrows(angle, r=1, adj=0.5, length=0.1, ...)
ms.sunflowers(n,r=0.3,adj=pi/2, ...)
ms.image(img, transpose=TRUE, ...)
ms.face(features, ...)

Arguments

Details

These functions/matricies can be passed as the symb argument to the my.symbols function. The represent examples that can be used to create your own symbols or may be used directly.

Value

These functions either return a 2 column matrix of points to be passed to lines or NULL.

Author(s)

Greg Snow [email protected]

See Also

my.symbols, polygon, arrows, lines, faces, also see rasterImage for an alternative to ms.image

Examples

plot(1:10,1:10)
my.symbols(1:10,1:10, ms.polygram, n=1:10, r=seq(0.5,1,length.out=10),
inches=0.3)

my.symbols(1:10,1:10, ms.polygon, n=1:10, add=FALSE, inches=0.3)

my.symbols(1:5, 5:1, ms.filled.polygon, add=FALSE, n=3:7, fg='green',
  bg=c('red','blue','yellow','black','white'), inches=0.3 )

my.symbols( 1:10, 1:10, ms.female, inches=0.3, add=FALSE)
my.symbols( 1:10, 10:1, ms.male, inches=0.3, add=TRUE)

plot(1:10, 1:10)
my.symbols(1:10, 1:10, ms.arrows, angle=runif(10)*2*pi, inches=0.5,
adj=seq(0,1,length.out=10), symb.plots=TRUE)

my.symbols(1:10, 1:10, ms.sunflowers, n=1:10, inches=0.3, add=FALSE)

if( require(png) ) {
  img <- readPNG(system.file("img", "Rlogo.png", package="png"))

  my.symbols( runif(10), runif(10), ms.image, MoreArgs=list(img=img),
                 inches=0.5, symb.plots=TRUE, add=FALSE)
}

tmp.mtcars <- scale(mtcars, center=sapply(mtcars,min),
	scale=sapply(mtcars,function(x) diff(range(x))) )
tmp2.mtcars <- lapply( seq_len(nrow(tmp.mtcars)), function(i) tmp.mtcars[i,] )
my.symbols(mtcars$wt, mtcars$mpg, ms.face, inches=0.3, features=tmp2.mtcars,
	add=FALSE)

Description

This function draws symbols on a plot. It is similar to the builtin symbols function with the difference that it plots symbols defined by the user rather than a prespecified set of symbols.

Usage

my.symbols(x, y=NULL, symb, inches=1, xsize, ysize,
                       add=TRUE,
                       vadj=0.5, hadj=0.5,
                       symb.plots=FALSE,
                       xlab=deparse(substitute(x)),
                       ylab=deparse(substitute(y)), main=NULL,
                       xlim=NULL, ylim=NULL, linesfun=lines,
                       ..., MoreArgs)

Arguments

Details

The symb argument can be a 2 column matrix or a list with components 'x' and 'y' that defines points on the interval [-1,1] that will be connected with lines to draw the symbol. If you want a closed polygon then be sure to replicate the 1st point as the last point.

If any point contains an NA then the line will not be drawn to or from that point. This can be used to create a symbol with disjoint parts that should not be connected.

If symb is a function then it should include a '...' argument along with any arguments to define the symbol. Any unmatched arguments that end up in the '...' argument will be replicated to the same length as 'x' (using the rep function) then the values will be passed one at a time to the symb function. If MoreArgs is specified, the elements of it will also be passed to symb without modification. The symb function can either return a matrix or list with the points that will then be passed to the lines function (see above). Or the function can call the plotting functions itself (set symb.plots to TRUE). High level plotting can be done (plot, hist, and other functions), or low level plotting functions (lines, points, etc) can be used; in this case they should add things to a plot with 'x' and 'y' limits of -1 to 1.

The size of the symbols can be specified by using inches in which case the symbol will be set inside of squares whose sizes are inches size based on the plotting device. The size can also be set using xsize and/or ysize which use the same units as the x and/or y variables. If only one is specified then the box will be square. If both are specified and they do not match the aspect ratio of the plot then the bounding box will not be square and the symbol will be distorted.

Value

This function is run for its side effect of plotting, it returns an invisible NULL.

Note

Since the '...' argument is passed to both lines and symb, the symb function should have a '...' argument so that it will ignore any additional arguments.

Arguments such as 'type' can be passed through the '...' argument if you want the symbol made of something other than lines.

Plotting coordinates and sizes are based on the size of the device at the time the function is called. If you resize the device after plotting, all bets are off.

Currently missing values in x or y are not handled well. It is best if remove all missing values first.

Author(s)

Greg Snow [email protected]

See Also

symbols, subplot, mapply, ms.polygram, lines

Examples

# symb is matrix

my.symbols( 1:10, runif(10), ms.male, add=FALSE, xlab='x',
  ylab='y', inches=0.3, col=c('blue','green'), xlim=c(0,11), ylim=c(-0.1,1.1))
my.symbols( (1:10)+0.5, runif(10), ms.female, add=TRUE, inches=0.3,
  col=c('red','green') )

# symb is function returning matrix

plot(1:10, 1:10)
my.symbols( 1:10, 1:10, ms.polygram, n=1:10, inches=0.3 )


# symb is plotting function
# create a variation on monthplot

fit <- lm( log(co2) ~ time(co2) )
fit.r <- resid(fit)

x <- 1:12
y <- tapply(fit.r, cycle(co2), mean)

tmp.r <- split( fit.r, cycle(co2) )
tmp.r <- lapply( tmp.r, function(x) x-mean(x) )

yl <- do.call('range',tmp.r)

tmpfun <- function(w,data,ylim,...){
  tmp <- data[[w]]
  plot(seq(along=tmp),tmp, type='l', xlab='',ylab='',
       axes=FALSE, ylim=ylim)
  lines(par('usr')[1:2], c(0,0), col='grey')
}

my.symbols(x,y, symb=tmpfun, inches=0.4, add=FALSE, symb.plots=TRUE,
  xlab='Month',ylab='Adjusted CO2', xlim=c(0.5,12.5),
  ylim=c(-0.012,0.012),
  w=1:12, MoreArgs=list(data=tmp.r,ylim=yl) )

# using xsize and ysize
plot( 1:10, (1:10)*100, type='n', xlab='', ylab='' )
my.symbols( 5, 500, ms.polygon, n=250, inches=1.5 )
my.symbols( 5, 500, ms.polygon, n=250, xsize=2, col='blue' )
my.symbols( 5, 500, ms.polygon, n=250, ysize=200, col='green' )
my.symbols( 5, 500, ms.polygon, n=250, xsize=2, ysize=200, col='red' )
abline( v=c(4,6), col='grey' )
abline( h=c(400, 600), col='grey' )


# hand crafted hexagonal grid

x1 <- seq(0, by=2*sqrt(3), length.out=10)
y1 <- seq(0, by=3, length.out=10)

mypoints <- expand.grid(x=x1, y=y1)
mypoints[,1] <- mypoints[,1] + rep( c(0,sqrt(3)), each=10, length.out=100 )

plot(mypoints, asp=1, xlim=c(-2,35))
my.symbols(mypoints, symb=ms.filled.polygon, n=6,
  inches=par('pin')[1]/(diff(par('usr')[1:2]))*4,
  bg=paste('gray',1:100,sep=''), fg='green' )

Description

This dataset is approximately bell shaped, but with some outliers. It is meant to be used for demonstration purposes. If students are tempted to throw out all outliers, then have them work with this data (or use a scaled/centered/shuffled version as errors in a regression problem) and see how many throw away 3/4 of the data before rethinking their strategy.

Usage

data(outliers)

Format

The format is: num [1:100] -1.548 0.172 -0.638 0.233 -0.228 ...

Details

This is simulated data meant to demonstrate "outliers".

Source

Simulated, see the examples section.

Examples

data(outliers)
qqnorm(outliers)
qqline(outliers)
hist(outliers)

o.chuck <- function(x) {  # function to throw away outliers
	qq <- quantile(x, c(1,3)/4, names=FALSE)
	r <- diff(qq) * 1.5
	tst <- x < qq[1] - r | x > qq[2] + r
	if(any(tst)) {
		cat('Removing ', paste(x[tst], collapse=', '), '\n')
		x <- x[!tst]
		out <- Recall(x)
	} else {
		out <- x
	}
	out
}

x <- o.chuck( outliers )
length(x)

if(require(MASS)) {
  char2seed('robust')
  x <- 1:100
  y <- 3 + 2*x + sample(scale(outliers))*10

  plot(x,y)
  fit <- lm(y~x)
  abline(fit, col='red')

  fit.r <- rlm(y~x)
  abline(fit.r, col='blue', lty='dashed')

  rbind(coef(fit), coef(fit.r))
  length(o.chuck(resid(fit)))
}



### The data was generated using code similar to:

char2seed('outlier')

outliers <- rnorm(25)

dir <- 1

while( length(outliers) < 100 ){
	qq <- quantile(c(outliers, dir*Inf), c(1,3)/4)
	outliers <- c(outliers,
		qq[ 1.5 + dir/2 ] + dir*1.55*diff(qq) + dir*abs(rnorm(1)) )
	dir <- -dir
}

Description

This function is similar to the pairs function, but instead of doing all pairwise plots, it takes 2 matricies or data frames and does all combinations of the first on the x-axis with the 2nd on the y-axis. Used with pairs and subsets can spread a scatterplot matrix accross several pages.

Usage

pairs2(x, y, xlabels, ylabels, panel = points, ..., row1attop = TRUE, gap = 1)

Arguments

Details

This functios is similar to the pairs function, but by giving it 2 sets of data it only does the combinations between them. Think of it as giving the upper right or lower left set of plots from pairs. If a regular scatterplot matrix is too small on the page/device then use pairs on subsets of the data to get the diagonal blocks of a scatterplot matrix and this function to get the off diagonal blocks.

Value

This function is run for the side effect of the plot. It does not return anything useful.

Note

Large amounts of the code for this function were blatently borrowed/stolen from the pairs function, the credit for the useful parts should go to the original authors, blame for any problems should go to me. This function is also released under GPL since much of it comes from GPL code.

Author(s)

Greg Snow, [email protected]

See Also

pairs, splom in the lattice package

Examples

pairs2(iris[,1:2], iris[,3:4], col=c('red','green','blue')[iris$Species])

# compare the following plot:
pairs(state.x77, panel=panel.smooth)

# to the following 4 plots

pairs(state.x77[,1:4], panel=panel.smooth)
pairs(state.x77[,5:8], panel=panel.smooth)
pairs2( state.x77[,1:4], state.x77[,5:8], panel=panel.smooth)
pairs2( state.x77[,5:8], state.x77[,1:4], panel=panel.smooth)

Description

This function draws symbols on a lattice plot. It is similar to the builtin symbols function with the difference that it plots symbols defined by the user rather than a prespecified set of symbols.

Usage

panel.my.symbols(x, y, symb, inches=1, polygon=FALSE,
                       ..., symb.plots=FALSE, subscripts, MoreArgs)

Arguments

Details

The symb argument can be a 2 column matrix or a list with components 'x' and 'y' that defines points on the interval [-1,1] that will be connected with lines to draw the symbol. If you want a closed polygon then be sure to replicate the 1st point as the last point or use the polygon option.

If any point contains an NA then the line will not be drawn to or from that point. This can be used to create a symbol with disjoint parts that should not be connected.

If symb is a function then any unmatched arguments that end up in the '...' argument will be replicated to the same length as 'x' (using the rep function) then the values will be passed one at a time to the symb function. If MoreArgs is specified, the elements of it will also be passed to symb without modification. The symb function can either return a matrix or list with the points that will then be passed to the llines function (see above).

Value

This function is run for its side effect of plotting, it returns an invisible NULL.

Note

Plotting coordinates and sizes are based on the size of the device at the time the function is called. If you resize the device after plotting, all bets are off.

Author(s)

Greg Snow [email protected]

See Also

symbols, my.symbols, subplot, mapply, ms.polygram, lines

Examples

if(require(lattice)) {
  tmpdf <- data.frame( x=1:10, y=1:10, g=factor(rep( c("A","B"), each=5 )),
	z=c(1:5,5:1) )

  xyplot( y ~ x, tmpdf, panel=panel.my.symbols, symb=ms.female, inches=0.3 )

  xyplot( y ~ x | g, tmpdf, panel=panel.my.symbols, symb=ms.male, inches=0.3)

  xyplot( y ~ x, tmpdf, panel=panel.superpose, groups=g,
	panel.groups= function(group.number, ...) {
		if(group.number==1) {
			panel.my.symbols(..., symb=ms.male)
		} else {
			panel.my.symbols(..., symb=ms.female)
		} },
		inches=0.3
	)

  xyplot( y ~ x, tmpdf, panel=panel.my.symbols, symb=ms.polygram, n=tmpdf$z,
		inches=0.3)

  xyplot( y ~ x | g, tmpdf, panel=panel.my.symbols, symb=ms.polygram,
		n=tmpdf$z, inches=0.3)

  xyplot( y ~ x, tmpdf, panel=panel.superpose, groups=g,
	panel.groups = panel.my.symbols,
		inches=0.3, symb=ms.polygon, n=tmpdf$z, polygon=TRUE,
		adj=rep(c(0,pi/4),5)
	)
}

Description

This plays the lateral thinking game Petals Around the Rose. This is a game where 5 regular dice are rolled and the players then try to figure out how many petals are around the rose.

Usage

petals(plot = TRUE, txt = TRUE)

Arguments

Details

At least one of the arguments plot and txt needs to be true, otherwise you will be guessing blind (or testing your psychic abilities).

The game is usually played with 5 physical dice, one person who knows the rules (the potentate of the rose, here the computer), and one or more players trying to learn the puzzle. The potentate can only give the players the following 3 rules:

1. The name of the game is "Petals Around the Rose" and the name is significant.

2. The answer is always 0 or an even number.

3. The potentate can tell the answer for any roll after any guesses are made.

The potentate (or other player) then rolls the 5 dice and any players are then allowed to guess. The potentate either confirms a correct guess or tells the correct answer, then the game continues with another roll. Players are not to discuss their reasoning so that each can solve it themselves. When a player thinks they have worked out the reasoning they demonstrate it by getting correct guesses, but not by discussing it with anyone. Generally 6 correct guesses in a row is considered evidence that they have figured out the rules and they are then considered a potentate of the rose.

For this implementation the computer will simulate the role of 5 dice and display the results and ask for a guess of how many petals are around the rose. The player then enters their guess and the computer then either confirms that it is correct or gives the correct answer.

Pressing enter without making a guess ends the game.

Value

This function only returns NULL, it is run for its side effects.

Note

Casual viewing of the function source code is unlikely to reveal the secret (and therefore this could be used as an example of one way to disguise portions of code from casual examination). More on disguising source code is at https://stat.ethz.ch/pipermail/r-devel/2011-October/062236.html.

Some basic debugging can reveal the secret, but that would be cheating and an admission that such a simple game has defeated you, so don't do it, just keep playing until you figure it out.

Author(s)

Greg Snow, [email protected]

References

http://www.borrett.id.au/computing/petals-bg.htm

See Also

dice

Examples

if(interactive()){
  petals()
}

Description

This function attempts to create a script that will recreate the current plot (in the graphics window). You can then edit any parts of the script that you want changed and rerun to get the modified plot.

Usage

plot2script(file='clipboard')

Arguments

Details

This function works with the graphics window and mainly traditional graphics (it may work with lattice or other graphics, but has not really been tested with those).

This function creates a script file (or puts it on the clipboard so that you can past into a script window or text editor) that will recreate the current graph in the current graph window. The script consists of very low level functions (calls to plot.window and axis rather than letting plot handle all this).

If you want the higher level functions that were actually used, then use the history or savehistory commands (this will probably be the better method for most cases).

Some of the low level plotting functions use different arguments to the internal version than the user callable version (box for example), the arguments to these functions may need to be editted before the full script will run correctly.

The lengths of command lines between the creation of the script and what can be run in R do not always match, you may need to manually wrap long lines in the script before it will run properly.

Value

This function is run for its side effects and does not return anything meaningful.

Note

For any serious projects it is best to put your code into a script to begin with and edit the original script rather than using this function.

This function depends on the recordPlot function which can change in any version. Therefore this function should not be considered stable.

Author(s)

Greg Snow [email protected]

See Also

history, savehistory, recordPlot, source

Examples

if(interactive()){

# create a plot
plot(runif(10),rnorm(10))
lines( seq(0,1,length=10), rnorm(10,1,3) )

# create the script
plot2script()

# now paste the script into a script window or text processor.
# edit the ranges in plot.window() and change some colors or
# other options.  Then run the script.
 }

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(hscale=1.5, vscale=1.5, wait=FALSE)
run.power.examp.old()

Arguments

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 returns a list with the parameter settings and the power if wait is TRUE.

run.power.examp.old does not return anything meaningful.

Author(s)

Greg Snow [email protected]

See Also

power.t.test

Examples

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

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

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 [email protected]

See Also

plot, cor

Examples

if(interactive()){
put.points.demo()

x <- rnorm(25, 5, 1)
y <- x + rnorm(25)
put.points.demo(x,y)
}

Description

Simulate and plot p-values from a normal or binomial based test under various conditions. When all the assumptions are true, the p-values should follow an approximate uniform distribution. These functions show that along with how violating the assumptions changes the distribution of the p-values.

Usage

Pvalue.norm.sim(n = 50, mu = 0, mu0 = 0, sigma = 1, sigma0 = sigma,
     test= c("z", "t"), alternative = c("two.sided", "less", "greater", "<>",
           "!=", "<", ">"), alpha = 0.05, B = 10000)
Pvalue.binom.sim(n=100, p=0.5, p0=0.5, test=c('exact','approx'),
                            alternative=c('two.sided', 'less', 'greater',
                            '<>','!=','<','>'),
                            alpha=0.05, B=1000)
run.Pvalue.norm.sim()
run.Pvalue.binom.sim()

Arguments

Details

These functions generate B samples from either a normal or binomial distribution, then compute the P-values for the test of significance on each sample and plot the P-values.

The run.Pvalue.norm.sim and run.Pvalue.binom.sim functions are GUI wrappers for the other 2 functions allowing you to change the parameters and click on "refresh" to run a new set of simulations.

Using NA for sigma0 will result in the sample standard deviations being used (leave blank in the GUI).

When the simulation conditions and the hypothesized values match, the distributions of the p-values will be approximately uniform. Changing the parameter of interest will show the idea of power. Changing the other parameters can show the effects of assumptions not being met.

Value

The P-values are invisibly returned.

Note

Note: the 2-sided p-values for the binomial may not match the results from binom.test and prop.test. The method used here is an approximation for speed.

Author(s)

Greg Snow, [email protected]

References

Murdock, D, Tsai, Y, and Adcock, J (2008) _P-Values are Random Variables_. The American Statistician. (62) 242-245.

See Also

t.test, z.test, binom.test, prop.test, tkexamp

Examples

if(interactive()) {
  run.Pvalue.norm.sim()
  run.Pvalue.binom.sim()
}

Description

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

Usage

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

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

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 [email protected]

See Also

dice, plot.dice, coin.faces, sample

Examples

if(interactive()){
rgl.coin()
flip.rgl.coin()
flip.rgl.coin(1)
flip.rgl.coin(2)

rgl.clear()

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

rgl.clear()

# letters instead of pictures
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()

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) ) )
}

Description

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

Usage

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

Arguments

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 [email protected]

See Also

rgl in package rgl, plot method in package sp

Examples

if(interactive()){

if(require("spData")) {
  data(world)
  rgl.Map(world$geom)
  spheres3d(0,0,0,.999, col='lightblue')
}
}

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

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 [email protected]

See Also

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

Examples

if(interactive()){
roc.demo()
with(CO2,
  roc.demo(uptake[Type=='Mississippi'],
           uptake[Type=='Quebec']       )
)
}

Description

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

Usage

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

Arguments

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 [email protected]

See Also

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

Examples

if(interactive()){
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)
}

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, vscale=1.5, hscale=1.5, wait=FALSE)
run.cor2.examp(n=100, seed, vscale=1.5, hscale=1.5, wait=FALSE)
run.old.cor.examp(n = 100, seed)
run.old.cor2.examp(n = 100, seed)

Arguments

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.

The 2 "old" versions use the default graphics device with a seperate window with the sliders, the versions without "old" in the name include the plot and sliders together in a single tk window.

The size of the plot can be changed by changing the values in the hscale and vscale boxes and clicking on the "Refresh" button.

Value

If wait is TRUE, then the return value is a list with the x and y values of the final plot.

If wait is FALSE (and in the "old" versions) an invisible NULL is returned.

Note

If wait is TRUE then R will wait until you click on the "Exit" button before you can use your R session again. If wait is FALSE then the tk window will appear, but R will regain control so that you can continue to use R as well as interact with the demonstration window.

Author(s)

Greg Snow [email protected]

See Also

cor, tkexamp

Examples

if(interactive()) {
  run.cor2.examp()
}

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

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 [email protected]

See Also

hist, slider

Examples

if(interactive()){
run.hist.demo( rnorm(250, 100, 5) )
}