Statistics Biological I Biometry 3010 / Natural Resources 3130 / Statistical Science 2200 Lab 3 Jennifer Chung Section 401 9/26/11 Lab 3 will focus on using R scripts to organize and apply commands. Scripts are files that contain commands that can be executed in part via highlighting or as a whole in much the same way as any computer program would function. For part of this lab we will need the Class2011.txt data, so read in this file so that it is in the R environment. class = read.table("Class2011.txt",header=T,sep=",") names(Class) 1. (5 points) To create a script file in Windows to to File/New Script; on the Mac go to New Documents. A blank script file should appear. Copy and paste the following commands into the Script Program Window:
stem(class$head) mean(class$head) sd(class$head) hist(class$head)
Now highlight just the first line in the Script Window, and if you are using Windows version of R click on Ctrl+R; in Mac Cmd+Return. You should see the command and the solution (i.e. the stem plot or the mean head size) show up in the R Console. Try this with the other commands. (Another approach is to just place the active cursor at the end of the line and do Ctrl+R or Cmd+Return to run the line.) In order to run all the commands at once, we typically highlight the whole thing and click on Ctrl+R (Cmd Return on the Mac). Copy and paste the output solutions into the box below. Make a histogram of the data, but without the mis-measured head sizes (if there are any). To do this, do something like ClassHead = Class$Head[Class$Head > ??] # What should ?? be? hist(ClassHead)
Make the plot relatively small to save space and paper. Calculate the mean plus or minus 2 standard deviations for this data. Indicate where these are relative to the data in the histogram by placing a mark on or below the figure where the two points occur. Approximately what percentage of the data lies between these marks? class = read.table("C:\\Users\\B30A\\Desktop\\Class2011.txt",header=T,sep=",") > names(class) [1] "head" "foot" "sex" "class" > stem(class$head) The decimal point is at the | 51 | 0 52 | 00 53 | 00 54 | 000000 55 | 00000000000 56 | 00000000000000000000000000 57 | 000000000000000000 58 | 00000000000000000000000 59 | 000000000000000 60 | 0000000000000 61 | 000 62 | 000 63 | 64 | 65 | 0 66 | 67 | 68 | 0 > mean(class$head) [1] 57.384 > sd(class$head) [1] 2.428938 > hist(class$head)
Histogram of class$head
Frequency
0 50
10
20
30
40
55
60 class$head
65
Histogram of classhead
Frequency
0
5
10
15
20
25
52
54
56 classhead
58
60
62
classhead = class$head[class$head < 63] > hist(classhead) mean(class$head)+(2*sd(class$head)) [1] 62.24188
> mean(class$head)-(2*sd(class$head)) [1] 52.52612 Approximately 95% of the data lies between these marks. Note that this script can be saved for future reference by using the File/Save command. This typically saves the file with a .R postscript, but you may have to add that yourself, so check on that. 2. (5 points) Create a script file to generate a key for identifying elements of a point plot, namely plot character type (pch), color (col), and character size (cex). Copy and paste the following code into a script file and execute using the run arrow button. What does the type="n" do? Eliminate it and see.
plot(seq(20),seq(20),type='n') for(i in seq(20)) points(i,i,pch=i,col=i,cex=i/2)
plot(seq(20),seq(20),type='n') > for(i in seq(20)) points(i,i,pch=i,col=i,cex=i/2)
seq(20)
5
10
15
20
5
10 seq(20)
15
20
plot(seq(20),seq(20)) > for(i in seq(20)) points(i,i,pch=i,col=i,cex=i/2)
seq(20)
5
10
15
20
5
10 seq(20)
15
20
3. (5 points) Add to the existing script file commands that will generate a key for identifying of elements a line plot, namely location of a horizontal line (h), color (col), and line width (lwd). Copy and paste the following code into the existing script file and execute using the run arrow button.
plot(seq(20),seq(20),type='n') for(i in seq(20)) abline(h=i,col=i,lwd=i/2)
plot(seq(20),seq(20),type='n') > for(i in seq(20)) abline(h=i,col=i,lwd=i/2)
seq(20)
5
10
15
20
5
10 seq(20)
15
20
Save the script containing the commands for creating these two keys in a file in your local directory called Key.R. This may be useful to have in the future when you are making figures for future assignments or in your own research.
4. (5 points) Generate probability distributions using the R functions Binomial (dbinom) and Poisson (dpois). For the parameters n and p this gives us the probabilities of getting the values 0 through 50 when the mean value is n*p=10 for each distribution. If you have time, try varying the parameters n and p a bit to see what alternate distributions you can produce. What parameter is usually used to represent the mean of the Poisson distribution? How closely does the Poisson approximate the binomial?
n = 50 p = 0.2 x = seq(0,n,by=1) y.binom = dbinom(x,n,p) y.pois = dpois(x,n*p) plot( x, y.binom) lines(x, y.pois, col=2)
n = 50 > p = 0.2 > x = seq(0,n,by=1) > y.binom = dbinom(x,n,p) > y.pois = dpois(x,n*p) > > plot( x, y.binom) > lines(x, y.pois, col=2)
y.binom
0.00 0
0.02
0.04
0.06
0.08
0.10
0.12
0.14
10
20 x
30
40
50
The poisson distribution approximates the binomial distribution fairly closely, with only a few exceptions around the highest values from when x is between 9 and 13. Lambda is used to represent the mean of the Poisson distribution.
5. (5 points) Generate 1000 new "random" values using the binomial and Poisson probability distribution functions and the parameters n and p. This is stochastic simulation. Apply a histogram to each data set generated and compare the results. How do the binomial and Poisson distributions compare for this set of parameter values? How do they compare to the plot of the theoretical values you made above. Make a table comparing the empirical mean and variance from each of your stochastic simulations with the theoretical mean and variance based on the parameters n and p for each distribution. (Make this a 2 x 4 table, mean and variance for the rows, theoretical and empirical of the binomial and Poisson for the columns.) How do the results compare? How do the variances compare between the binomial and Poisson?
n = 50 p = 0.2 r.binom = rbinom(1000,size=n,prob=p) r.pois = rpois(1000,lambda=n*p) # This part will make the xscale the same r.breaks = hist(c(r.binom,r.pois),plot=F)$breaks par(mfrow=c(2,1)) hist(r.binom,breaks=r.breaks) hist(r.pois,breaks=r.breaks) par(mfrow=c(1,1))
n = 50 > p = 0.2 > > r.binom = rbinom(1000,size=n,prob=p) > r.pois = rpois(1000,lambda=n*p) > r.breaks = hist(c(r.binom,r.pois),plot=F)$breaks > par(mfrow=c(2,1)) > hist(r.binom,breaks=r.breaks) hist(r.pois,breaks=r.breaks)
Histogram of r.binom
250 Frequency 0 50 150
5
10 r.binom
15
20
Histogram of r.pois
250 Frequency 0 50 150
5
10 r.pois
15
20
The binomial and poisson distributions are almost identical, with a few differences in the frequency of values greater than 14. However, they are both fairly symmetric and unimodal. Compared to the theoretical plot above, the histograms have the same frequencies, even with a greater n value.
Save the commands for Questions 4 and 5 into a script labeled BinomPois.R. Doing so will allow you to access it easier later. Make sure you spend some time trying to understand what the scripts are doing. Make substitutions for the values in the scripts, and add or eliminate commands to get a better feel for what is going on.