# R examples for Statistics 572 Lecture
# January 25, 2008
# Bret larget

# All text after a # is a comment that R ignores.

# All of these R commands are described in the file
# lecture2.pdf

# You can have R execute these commands by typing
#
#   source("lecture2.R")
#
# at the R prompt.
#
# For this to work, you need to change the working directory
# to the directory where this file is stored
# and this file would need to be named lecture2.R.

# Note also that when R reads and executes a list of
# commands in a file, it will only print to screen
# or create graphics plots when there is an explicit
# call to print or plot.
#
# So, for example, the command
#
#   dbinom(7,13,0.6)
#
# will compute this number, but it will neither save nor print the result.
# We could save it as a new object called a with this command
#
#   a = dbinom(7,13,0.6)
#
# and use it later. If we just want it to print, do this.
#
#   print( dbinom(7,13,0.6) )
#

#
# Examples of probability calculations
#

# Binomial probability of 7 sucesses, n=13, p=0.6
dbinom(7,13,0.6)

# Poisson probabilities of 0:2, mu=5
dpois(0:2,5)

# Normal curve area to the left of 93 with mu=100, sigma=15
pnorm(93,100,15)

# 0.1 and 0.9 quantiles from normal curve with mu=100, sigma=15
qnorm(c(0.1,0.9),100,15)

# a random sample of 12 independent binomial random variables with n=10, p=0.3
rbinom(12,10,0.3)

# demonstration of the set.seed function.
set.seed(3451)
x1 = rnorm(5)
x2 = rnorm(5)
set.seed(3451)
x3 = rnorm(5)
x1
x2
x3

#
# Reading in data to a data frame.
#

fevdata = read.table("fevdata.txt",header=T)
str(fevdata)

# Calculations of a variable in the data frame.
with(fevdata, mean(age))
with(fevdata, sd(age))
with( fevdata, sapply(split(ht,sex),mean) )

#
# Loading library for lattice graphics
#

library(lattice)

# Histogram of height
plot( histogram( ~ ht, data=fevdata ) )

# Same histogram with a different color and number of bins.
plot( histogram( ~ ht, col="red", nint=24, data=fevdata ) )

# A function to be used as the panel argument in the histogram function
# to overlay a normal curve over the histogram.

overlayNormal = function(x,ncol="black",nlwd=2,...) {
  panel.histogram(x,...)
  panel.mathdensity(dmath = dnorm, col=ncol, args = list(mean=mean(x),sd=sd(x)),lwd=nlwd)
}

# The actual call to create the plot.
plot( histogram( ~ ht, type="density", panel=overlayNormal, data=fevdata ) )

# Two separate histograms, one for each sex
plot( histogram( ~ ht | sex, data=fevdata ) )

# Side-by-side box-and-whisker plots
plot( bwplot( ht ~ sex, data=fevdata ) )

# A scatterplot with a key
plot( xyplot( ht ~ age, group=sex, auto.key=T,data=fevdata ) )