# ---- A small sample of using R ---- # 
a <- 3;
b <- -1.5;
n1 <- 100;
n2=10;
windows(); opar <- par(mfrow = c(2,2), oma = c(0, 0, 2.7, 0)); 
# 4 plots in one graphic window ---   can you have, say, 6 plots in one window?  #

x <- seq(-5, 5, len=n1); # n1 equally spaced points in [-5, 5] #
x;
z <- a+b*x;          # a: intercept, b: slope #
plot(x,z,main=paste("sample size=", n1));
sd1 <-2
y<-z+rnorm(n1)*sd1;   # rnorm(M) gives you M random numbers from N(0,1) #
plot(x,y, main=paste("sample size=",n1, ",s.d.=", sd1));

hist(z,freq=FALSE, main=paste("sample size=",n1, ",s.d.=", 0))   # Histogram of z's #
hist(y,freq=FALSE, main=paste("sample size=",n1, ",s.d.=", sd1))   # Do you see any difference in these two histograms? #


windows(); opar <- par(mfrow = c(2,2), oma = c(0, 0, 2.7, 0))
y<-z+rnorm(n1)*sd1;
plot(x,y, main=paste("sample size=",n1, ",s.d.=", sd1));
fm <- lm(y~x);       # Simple regression; y regress on x # 
summary(fm);         # Check if estimates for intercept and slope good #
abline(fm, lty=3,col="red");           # Add the estimated reg. line to the current plot #
abline(a, b, lty=1,col=1);  # The true regression line #


sd2 <- 16;
t<-z+rnorm(n1)*sd2;
plot(x,t,main=paste("sample size=",n1, ",s.d.=", sd2));
abline(a, b, lty=1, col=1);
fm1 <- lm(t~x);
summary(fm1);
abline(fm1, lty=3,col="blue");

hist(y,freq=FALSE, main=paste("sample size=",n1, ",s.d.=", sd1))
hist(t,freq=FALSE, main=paste("sample size=",n1, ",s.d.=", sd2))


windows(); opar <- par(mfrow = c(2,2), oma = c(0, 0, 2.7, 0))

x1 <- seq(-5, 5, len=n2)
z1 <- a+b*x1;
#plot(x1,z1,main=paste("sample size=",n2));#

y1<-z1+rnorm(n2)*sd1;
plot(x1,y1,main=paste("sample size=",n2, ",s.d.=", sd1));
fm2 <- lm(y1~x1);
summary(fm2);
abline(fm2,lty=3,col="red")
abline(a, b, lty=1,col=1);


t1<-z1+rnorm(n2)*sd2;
plot(x1,t1,main=paste("sample size=",n2, ",s.d.=", sd2));
fm3 <- lm(t1~x1);
abline(fm3,lty=3,col="blue")
abline(a, b, lty=1,col=1);
summary(fm3);
hist(y1,freq=FALSE, main=paste("sample size=",n2, ",s.d.=", sd1))
hist(t1,freq=FALSE, main=paste("sample size=",n2, ",s.d.=", sd2))


windows(); opar <- par(mfrow = c(2,2), oma = c(0, 0, 2.7, 0))
y<-z+rnorm(n1)*sd1;
plot(x,y, main=paste("sample size=",n1, ",s.d.=", sd1));
fm <- lm(y~x);       # Simple regression; y regress on x # 
summary(fm);         # Check if estimates for intercept and slope good #
abline(fm, lty=3,col="red");           # Add the estimated reg. line to the current plot #
abline(a, b, lty=1,col=1);  # The true regression line #
anova(fm);              # Obtain the ANOVA table    #
sd2 <- 16;
t<-z+rnorm(n1)*sd2;
plot(x,t,main=paste("sample size=",n1, ",s.d.=", sd2));
abline(a, b, lty=1, col=1);
fm1 <- lm(t~x);
summary(fm1);
abline(fm1, lty=3,col="blue");
anova(fm1);             #the ANOVA table    #

y1<-z1+rnorm(n2)*sd1;
plot(x1,y1,main=paste("sample size=",n2, ",s.d.=", sd1));
fm2 <- lm(y1~x1);
summary(fm2);
abline(fm2, lty=3,col="red")
abline(a, b, lty=1,col=1);
anova(fm2);                #the ANOVA table    #

t1<-z1+rnorm(n2)*sd2;
plot(x1,t1,main=paste("sample size=",n2, ",s.d.=", sd2));
fm3 <- lm(t1~x1);
abline(fm3, lty=3,col="blue")
abline(a, b, lty=1,col=1);
summary(fm3);
anova(fm3);                #the ANOVA table    #
# Can you see the effect of sample size and standard deviation on the fitted line? #
# Play with R as much as possible this week, you will need to use R 
# for your homework shortly.    

# Always remove all objects before ending your R session # 
rm(list=ls(all=TRUE));

q() # To quit R #

Come to my office if you have any question!    ^_^  Yes.