## Chapter 3 

## Sec 3.1: Regularities 
## ss 3.1.1: Deterministic models 
## Footnote code 
## Simplified version of graph 
curve(0.5*9.8*x^2, from=0, to=5, xaxs="i", yaxs="i",  
      xlab = "Time (sec)", ylab = "Distance fallen (m)") 
text(0.25, 105, expression("Distance" == frac(1,2)*phantom(i)*g*t^2), adj=0) 

## ss 3.1.2: Models that include a random component 
## Footnote code 
## Plot depression vs weight: data frame roller (DAAG) 
plot(depression ~ weight, data = roller, xlim=c(0,1.04*max(weight)), 
     ylim=c(0,1.04*max(depression)),  
     xaxs="i", yaxs="i",   # "i"=inner: Fit axes exactly to the limits 
     xlab = "Weight of roller (t)", ylab = "Depression(mm)", pch = 16)  
abline(0, 2.25)            # A slope of 2.25 looks about right 

##                      Generalizing from models 
##                        Which model is best? 
## ss 3.1.3: Fitting models -- the model formula 
## Fit line - by default, this fits intercept & slope. 
## requires data frame roller (DAAG) 
roller.lm <- lm(depression ~ weight, data=roller) 
## Compare with the code used to plot the data 
plot(depression ~ weight, data=roller) 
## Add a the fitted line to the plot 
abline(roller.lm) 
## Footnote code 
## For a model that omits the intercept term, specify 
lm(depression ~ -1 + weight, data=roller)   

##                    Fitted values and residuals 
##                           Model objects 
names(roller.lm)     # Get names of list elements 
coef(roller.lm) 
roller.lm$coefficients 

## Sec 3.2: Distributions: Models for the Random Component 

## ss 3.2.1: Discrete distributions 
##                       Bernoulli distribution 
##                       Binomial distribution 
## Footnote code 
## To get the labeling (0, 1, 2) as in the text, specify: 
probs <- dbinom(0:2, size=2, prob=0.5) 
names(probs) <- 0:2 
probs 
pbinom(q=2, size=4, prob=0.5) 
pbinom(q=4, size=50, prob=0.2) 
qbinom(p = 0.65, size = 4, prob = 0.5) 
## Check result 
sum(dbinom(0:2, size=4, prob=.5)) 

##                        Poisson distribution 

##              Means, variances and standard deviations 
## ss 3.2.2: Continuous distributions 
##                        Normal distribution 
## Plot the normal density, in the range -3 to 3 
z <- pretty(c(-3,3), 30)   # Find ~30 equally spaced points 
ht <- dnorm(z)             # By default: mean=0, variance=1 
plot(z, ht, type="l", xlab="Normal deviate", ylab="Density", yaxs="i") 
 # yaxs="i" locates the axes at the limits of the data 
## Footnote code 
## Plot the normal density, in the range -3.25 to 3.25 
z <- pretty(c(-3.25,3.25), 30) # Find ~30 equally spaced points 
ht <- dnorm(z)                 # By default: mean=0, variance=1 
plot(z, ht, type="l", xlab="Normal deviate", ylab="Ordinate", yaxs="i") 
polygon(c(z[z <= 1.0], 1.0), c(dnorm(z[z <= 1.0]), 0), col="grey") 
 # Around 84.1% of the total area is to the left of the vertical line. 
pnorm(1.0)         # by default, mean=0 and SD=1 
## Footnote code 
## Additional examples 
pnorm(0)             # .5    (exactly half the area is to the left of the mean) 
pnorm(-1.96)         # .025 
pnorm(1.96)          # .975 
pnorm(1.96, mean=2)  # .484  (a normal distribution with mean 2 and SD 1) 
pnorm(1.96, sd=2)    # .836  (sd = standard deviation) 
qnorm(.9)          # 90th percentile; mean=0 and SD=1 
## Footnote code 
## Additional examples 
qnorm(0.841)         # 1.0 
qnorm(0.5)           # 0 
qnorm(0.975)         # 1.96 
qnorm(c(.1,.2,.3))   # -1.282 -0.842 -0.524  (10th, 20th and 30th percentiles) 
qnorm(.1, mean=100, sd=10)  # 87.2 (10th percentile, mean=100, SD=10) 

##                   Other continuous distributions 
##              Different ways to describe distributions 

## Sec 3.3: The Uses of Random Numbers 
## ss 3.3.1: Simulation 
rbinom(10, size=1, p=.5)   # 10 Bernoulli trials, prob=0.5 
rbinom(25, size=4, prob=0.5)   
set.seed(9388) 
rpois(20, 3) 
options(digits=2)  # Suggest number of digits to display 
rnorm(10)          # 10 random values from the normal distribution 
runif(n = 20, min=0, max=1) # 20 numbers, uniform distn on (0, 1) 
rexp(n=10, rate=3)          # 10 numbers, exponential, mean 1/3. 

## ss 3.3.2: Sampling from populations 
## For the sequence below, precede with set.seed(3676) 
sample(1:9384, 15, replace=FALSE) 
## For the sequence below, precede with set.seed(366) 
split(sample(seq(1:10)), rep(c("Control","Treatment"), 5)) 
 # sample(1:10) gives a random re-arrangement (permutation)  
 # of 1, 2, ..., 10 

##                          Cluster sampling 
##                         Resampling methods 

## Sec 3.4: Model Assumptions 
## ss 3.4.1: Random sampling assumptions -- independence 
## ss 3.4.2: Checks for normality 
##                     Examination of histograms 
## Footnote code 
## The following gives a rough equivalent of the figure: 
set.seed (21) # Use to reproduce the data in the figure 
par(mfrow=c(2,3)) 
x <- pretty(c(6.5,13.5), 40) 
for(i in 1:5){ 
    y <- rnorm(50, mean=10, sd=1) 
    hist(y, prob=TRUE, xlim=c(6.5,13.5), ylim=c(0,0.5), main="") 
    lines(x, dnorm(x,10,1)) 
    } 
par(mfrow=c(1,1)) 
rm(x, y) 

##                    The normal probability plot 
## Use qreference() (DAAG) 
## With seed=21, the random numbers are as in the previous figure 
qreference(m=50, seed=21, nrep=5, nrows=1)  # 50 values per panel 
## Footnote code 
set.seed(21) # Use the same setting as for the previous figure 
library(lattice) 
qqmath(~rnorm(50*5)|rep(1:5,rep(50,5)), layout=c(5,1), aspect=1) 

##    The sample plot, set alongside plots for random normal data 
## Compare normal probability plot for normal-ambient difference 
## with simulated normal values: data frame pair65 (DAAG) 
qreference(pair65$heated - pair65$ambient, nrep=8) 

##             Formal statistical testing for normality? 
## ss 3.4.3: Checking other model assumptions 
## ss 3.4.4: Are non-parametric methods the answer? 
## ss 3.4.5: Why models matter -- adding across contingency tables 

## Sec 3.5: Recap 

## Sec 3.6: Further Reading 
x <- rpois(7, 78.3) 
mean(x); var(x) 
Markov <- function (N=100, initial.value=1, P) 
{ 
    X <- numeric(N) 
    X[1] <- initial.value + 1 
    n <- nrow(P) 
    for (i in 2:N){ 
    X[i] <- sample(1:n, size=1, prob=P[X[i-1],])} 
    X - 1 
}