################################################### ### chunk number 1: ################################################### scores <- c(99,97,72,56,88,80,74,95,66,57,89) # create data set max(scores) min(scores) ################################################### ### chunk number 2: ################################################### mean(scores) ################################################### ### chunk number 3: ################################################### var(scores) ################################################### ### chunk number 4: ################################################### sd(scores) sqrt(var(scores)) ################################################### ### chunk number 5: ################################################### median(scores) ################################################### ### chunk number 6: ################################################### quantile(scores,.25) IQR(scores) fivenum(scores) summary(scores) ################################################### ### chunk number 7: ################################################### boxplot(scores) ################################################### ### chunk number 8: ################################################### hist(scores) ################################################### ### chunk number 9: ################################################### rbinom(1,1,0.5) ################################################### ### chunk number 10: ################################################### rbinom(2,1,0.5) ################################################### ### chunk number 11: ################################################### # Observe right-stone hits for 40 raindrops, but do this # only once (=one trial); assuming ideal conditions, we can # safely suppose that the probability of a right-stone # hit = 0.5 size <- 1 p <- 0.5 fortydrops <- rbinom(40,size,p) ################################################### ### chunk number 12: ################################################### #find out the total number of right-stone hits (or 1's) sum(fortydrops)/40 ################################################### ### chunk number 13: ################################################### op <- par(mfrow=c(2,3),pty = "s") ################################################### ### chunk number 14: ################################################### size <- 1 p <- 0.5 k <- 40 # create a vector "observations" observations <- c(15,25,50,100,500,1000) n <- length(observations) ################################################### ### chunk number 15: ################################################### op <- par(mfrow=c(2,3),pty = "s") for(j in 1:n){ results <- rep(NA,observations[j]) for(i in 1:observations[j]){ results[i] <- sum(rbinom(k,size,p)) } ## end-for # create title on the fly: title <- paste(c("Num Obs.",observations[j],sep=" ")) # plot histogram: hist(results, ylab="Frequency", xlab="No. of R-stone hits", main=title) }#end-for ################################################### ### chunk number 16: ################################################### size <- 1 p <- 0.5 drops <- c(4,40,400,4000) op <- par(mfrow=c(2,2),pty = "s") for(num.drops in drops) { results <- rep(NA,100) for (i in 1:100){ results[i] <- sum(rbinom(num.drops,size,p)) }#end-for maintitle <- paste(num.drops,"drops",sep=" ") hist(results, xlim=range(c(0:num.drops)), xlab="Number of R-stone hits",main=maintitle) }#end-for ################################################### ### chunk number 17: ################################################### drops<-rep(1:400,1) standard.dev <- rep(NA,400) for(j in 1:400) { results <- rep(NA,100) for (i in 1:100){ results[i] <- sum(rbinom(drops[j],size,p)) } standard.dev[j] <- sd(results) } plot(drops,standard.dev, xlim=c(1,400), xlab="Number of drops", ylab="Standard deviation", main=expression("SD "*italic("increases")*" as we increase sample size.")) ################################################### ### chunk number 18: ################################################### size <- 1 p <- 0.5 drops <- c(4,40,400,4000) op <- par(mfrow=c(2,2),pty = "s") for(num.drops in drops) { results <- rep(NA,100) for (i in 1:100){ results[i] <- mean(rbinom(num.drops,size,p)) }#end-for maintitle <- paste(num.drops,"drops",sep=" ") hist(results, xlim=range(c(0:1)), xlab="Proportion of R-stone hits",main=maintitle) }#end-for ################################################### ### chunk number 19: ################################################### drops<-rep(1:400,1) standard.dev <- rep(NA,400) for(j in 1:400) { results <- rep(NA,100) for (i in 1:100){ results[i] <- mean(rbinom(drops[j],size,p)) } standard.dev[j] <- sd(results) } plot(drops,standard.dev, xlim=c(1,400), xlab="Number of drops", ylab="Standard deviation", main=expression("SD now decreases as we increase sample size.")) ################################################### ### chunk number 20: ################################################### size <- 1 p <- 0.5 num.drops <- 4000 results <- rep(NA,100) for (i in 1:100){ results[i] <- sum(rbinom(num.drops,size,p)) } mean.results<-mean(results) n<-floor(mean.results-1) m<-floor(mean.results+1) # candidates we can subtract deviations from: xvalues <- c(1:n,mean.results,m:4000) totaldev <- rep(NA,length(xvalues)) for(i in xvalues){ vectori <- rep(i,100) diffs <- results - vectori sqdeviations <- sum(diffs*diffs) totaldev[i] <- sqdeviations } plot(xvalues,totaldev, xlab="Potential minimizers of squared deviation", ylab="Squared Deviation", main="Squared deviations from the mean versus other numbers") lines(xvalues,totaldev) ################################################### ### chunk number 21: ################################################### #rbinom(40,1,.5) # 40 raindrops observed once binomialprobability <- function(n,p,k){ choose(n,k)*p^k*(1-p)^(n-k) } ################################################### ### chunk number 22: ################################################### binomialprobability(100,2/3,66) ################################################### ### chunk number 23: ################################################### dbinom(66,100,2/3) ################################################### ### chunk number 24: ################################################### binomialprobability(60,2/3,40) ################################################### ### chunk number 25: ################################################### binomialprobability(600,2/3,400) ################################################### ### chunk number 26: ################################################### numballs <- 40 p <- 0.5 probs <- rep(NA,numballs) for(k in 1:numballs){ probs[k] <- binomialprobability(numballs,p,k) } ################################################### ### chunk number 27: ################################################### probs2<-rep(NA,40) for(i in 1:40){ probs2[i]<-dbinom(i,40,0.5) } ################################################### ### chunk number 28: ################################################### head(probs) ################################################### ### chunk number 29: ################################################### (withinone <- sum(probs[19:21])) (withintwo <- sum(probs[18:22])) ################################################### ### chunk number 30: ################################################### # The 95% CI: # Looking at probabilities around a mean value: mean.index <- 20 intervals <- rep(NA,19) for(i in 1:19){ indices <- seq(mean.index-i,mean.index+i,by=1) range <- probs[indices] intervals[i] <- sum(range) } conf.intervals <- data.frame(margin=rep(1:19),probability=intervals) conf.intervals40 <- conf.intervals print(head(conf.intervals)) ################################################### ### chunk number 31: ################################################### plot(conf.intervals$margin,conf.intervals$probability,type="b",xlab="Margins",ylab="Probability",main="Sample size 40") segments(0,.95,5.7,.95,col="gray") segments(5.7,0,5.7,.95,col="gray") ################################################### ### chunk number 32: ################################################### numballs <- 400 p <- 0.5 probs <- rep(NA,numballs) for(k in 1:numballs){ currentk <- binomialprobability(numballs,p,k) probs[k] <- currentk } mean.index <- 200 intervals <- rep(NA,199) for(i in 1:199){ indices <- seq(mean.index-i,mean.index+i,by=1) range <- probs[indices] intervals[i] <- sum(range) } conf.intervals <- data.frame(margin=rep(1:199),probability=intervals) conf.intervals400 <- conf.intervals print(head(conf.intervals)) ################################################### ### chunk number 33: ################################################### plot(conf.intervals$margin,conf.intervals$probability,type="b",xlab="Margins",ylab="Probability",main="Sample size 40") segments(-6,.95,19,.95,col="gray") segments(19,0,19,.95,col="gray") ################################################### ### chunk number 34: ################################################### compute_margins <- function(numballs,p){ probs <- rep(NA,numballs) for(k in 1:numballs){ currentk <- binomialprobability(numballs,p,k) probs[k] <- currentk } mean.index <- numballs*p max.margin <- numballs*p-1 #max.margin <- numballs intervals <- rep(NA,max.margin) for(i in 1:max.margin){ indices <- seq(mean.index-i,mean.index+i,by=1) range <- probs[indices] intervals[i] <- sum(range) } conf.intervals <- data.frame(margin=rep(1:max.margin),probability=intervals) return(conf.intervals) } ################################################### ### chunk number 35: ################################################### plotcis <- function(numballs,p,color="black",margin,maintitle,interval=TRUE){ probs <- rep(NA,numballs+1) for(k in 0:numballs){ currentk <- binomialprobability(numballs,p,k) probs[k+1] <- currentk } proportions <- 0:numballs/numballs plot(proportions,probs,type="l",col="black",xlab="Proportions", ylab="Probability",main=maintitle) if(interval==TRUE){ # draw 95% confidence intervals segments(proportions[(numballs/2+1)-margin],-0.5, proportions[(numballs/2+1)-margin],.06,col=color,lty=1,lwd=2) segments(proportions[(numballs/2+1)+margin],-0.5, proportions[(numballs/2+1)+margin],.06,col=color,lty=1,lwd=2) } } ################################################### ### chunk number 36: ################################################### op <- par(mfrow=c(1,2),pty = "s") plotcis(40,.5,margin=5,maintitle="Sample size 40") plotcis(400,.5,margin=19,maintitle="Sample size 400") ################################################### ### chunk number 37: ################################################### sums<-rep(NA,21) for(i in 0:20){ sums[i+1]<-dbinom(i,40,0.5) } sum(sums) ################################################### ### chunk number 38: ################################################### sum(dbinom(0:20,40,0.5)) ################################################### ### chunk number 39: ################################################### pbinom(20,40,0.5) ################################################### ### chunk number 40: ################################################### newfunction <- function(x,mu,sigma){ 1/(sqrt(2*pi)*sigma)*exp(1)^(-((x - mu)^2/(2*sigma^2)))} ################################################### ### chunk number 41: ################################################### plotcis(40,.5,40,margin=20,maintitle="Comparing the binomial and normal distributions",interval=FALSE) lines(c(1:40)/40,newfunction(c(1:40),20,3),col="black",lty=2) ################################################### ### chunk number 42: ################################################### integrate(dnorm,-1,1) integrate(dnorm,-1.96,1.96) integrate(dnorm,-2,2) ################################################### ### chunk number 43: ################################################### #1000 samples of 40 taken repeatedly: means <- rep(NA,1000) for(i in c(1:1000)){ sample40 <- rnorm(40,mean=60,sd=8) means[i] <- mean(sample40) } ################################################### ### chunk number 44: ################################################### (means40 <- mean(means)) (sd40 <- sd(means)) ################################################### ### chunk number 45: ################################################### hist(means) ################################################### ### chunk number 46: ################################################### samplesize <- 100 # was 40 means <- rep(NA,1000) for(i in c(1:1000)){ sample100 <- rnorm(samplesize,mean=60,sd=8) means[i] <- mean(sample100) } ################################################### ### chunk number 47: ################################################### (means100 <- mean(means)) (sd100 <- sd(means)) ################################################### ### chunk number 48: ################################################### hist(means) ################################################### ### chunk number 49: ################################################### 8/sqrt(40) 8/sqrt(100) ################################################### ### chunk number 50: ################################################### population <- rexp(1000) hist(population) ################################################### ### chunk number 51: ################################################### means <- rep(NA,1000) for(i in c(1:1000)){ samp100 <- rexp(100) means[i] <- mean(samp100) } hist(means) ################################################### ### chunk number 52: ################################################### sample11 <- rnorm(11,mean=60,sd=4) ################################################### ### chunk number 53: ################################################### (estimated.mean <- mean(sample11)) popSD <- 4 sample.size <- length(sample11) (sigma.mu <- popSD/sqrt(sample.size)) ################################################### ### chunk number 54: ################################################### sample44 <- rnorm(44,mean=60,sd=4) estimated.mean <- mean(sample44) sample.size <- length(sample44) (sigma.mu <- 4/sqrt(sample.size)) ################################################### ### chunk number 55: ################################################### sample.sd <- rep(NA,1000) for(i in c(1:1000)){ sample40 <- rnorm(40,mean=60,sd=8) sample.sd[i] <- sd(sample40) } hist(sample.sd) ################################################### ### chunk number 56: ################################################### sample.sd <- rep(NA,1000) for(i in c(1:1000)){ sample40 <- rexp(40) sample.sd[i] <- sd(sample40) } hist(sample.sd) ################################################### ### chunk number 57: ################################################### range <- seq(-4,4,.01) multiplot(2,3) for(i in c(2,5,10,15,20,50)){ plot(range,dnorm(range),lty=1,col="gray") lines(range,dt(range,df=i),lty=2) mtext(paste("df=",i),cex=1.2) } ################################################### ### chunk number 58: ################################################### sample <- rnorm(11,mean=60,sd=4) ################################################### ### chunk number 59: ################################################### print(t.test(sample)$conf.int) ################################################### ### chunk number 60: ################################################### sample <- rnorm(100,mean=60,sd=4) print(t.test(sample)$conf.int) ################################################### ### chunk number 61: ################################################### se <- function(x) { y <- x[!is.na(x)] # remove the missing values, if any sqrt(var(as.vector(y))/length(y)) } ci <- function (scores){ m <- mean(scores,na.rm=TRUE) stderr <- se(scores) len <- length(scores) upper <- m + qt(.975, df=len-1) * stderr lower <- m + qt(.025, df=len-1) * stderr return(data.frame(lower=lower,upper=upper)) } ################################################### ### chunk number 62: ################################################### lower <- rep(NA,100) upper <- rep(NA,100) for(i in 1:100){ sam <- rnorm(100,mean=60,sd=4) lower[i] <- ci(sam)$lower upper[i] <- ci(sam)$upper } cis <- cbind(lower,upper) head(cis) ################################################### ### chunk number 63: ################################################### store <- rep(NA,100) for(i in 1:100){ sam <- rnorm(100,mean=60,sd=4) if(ci(sam)$lower<60 & ci(sam)$upper>60){ store[i] <- TRUE} else { store[i] <- FALSE } } summary(store) ################################################### ### chunk number 64: ################################################### # re-define variance to see whether it is an unbiased estimator newvar <- function(x){ m <- rep(mean(x),length(x)) d <- (x - m)^2 return(sum(d)/length(x)) } correct <- rep(NA,1000) incorrect <- rep(NA,1000) for(i in c(1:1000)){ sample10 <- rnorm(10,mean=0,sd=1) correctvar <- var(sample10) incorrectvar <- newvar(sample10) correct[i] <- correctvar incorrect[i] <- incorrectvar } ################################################### ### chunk number 65: ################################################### op <- par(mfrow=c(1,2)) hist(correct,main=paste("Mean ",round(mean(correct),digits=2),sep=" ")) hist(incorrect,main=paste("Mean ",round(mean(incorrect),digits=2),sep=" ")) ################################################### ### chunk number 66: ################################################### for(i in c(1:1000)){ sample100 <- rnorm(100,mean=0,sd=1) correctvar <- var(sample100) incorrectvar <- newvar(sample100) correct[i] <- correctvar incorrect[i] <- incorrectvar } op <- par(mfrow=c(1,2)) hist(correct,main=paste("Mean",round(mean(correct),digits=2),sep=" ")) hist(incorrect,main=paste("Mean",round(mean(incorrect),digits=2),sep=" ")) ################################################### ### chunk number 67: ################################################### sample <- rnorm(11,mean=60,sd=4) sample.mean <- mean(sample) sigma.mu <- 4/sqrt(11) ################################################### ### chunk number 68: ################################################### range <- seq(55,85,0.01) plot(range,dnorm(range,mean=70, sd=sigma.mu),type="l",ylab="", col="red",main="The null hypothesis") ################################################### ### chunk number 69: ################################################### (z <- (sample.mean - 70)/(4/sqrt(11))) sample <- rnorm(11,mean=60,sd=4) t.test(sample, alternative = "two.sided", mu = 70, conf.level = 0.95) ################################################### ### chunk number 70: ################################################### d <- rep(NA,1000) for(i in c(1:1000)){ sample1 <- rnorm(11,mean=60,sd=4) sample2 <- rnorm(15,mean=62,sd=6) d[i] <- mean(sample1) - mean(sample2) } hist(d) ################################################### ### chunk number 71: ################################################### sample1 <- rnorm(11,mean=60,sd=4) sample2 <- rnorm(15,mean=62,sd=6) t.test(sample1,sample2, mu=0, alternative = "two.sided", conf.level = 0.95,var.equal=FALSE) ################################################### ### chunk number 72: ################################################### # power visualization d <- c() for(i in c(1:1000)){ sample1 <- rnorm(11,mean=60,sd=4) sample2 <- rnorm(15,mean=62,sd=6) currentd <- mean(sample1) - mean(sample2) d <- append(d,currentd) } xvals <- seq(-6,6,.1) plot(xvals,dnorm(xvals,mean=0,sd=1.9633),type="l",lwd=2,ylab="",col="red") arrows(-(2*1.9633),-.05,-(2*1.9633),0.2,angle=0) arrows((2*1.9633),-.05,(2*1.9633),0.2,angle=0) text(-4.5,.008,expression(alpha/2),cex=1.5,col="red") text(4.5,.008,expression(alpha/2),cex=1.5,col="red") text(0,.008,expression(1-alpha),cex=1.5,col="red") ################################################### ### chunk number 73: ################################################### plot(xvals,dnorm(xvals,mean=0,sd=1.9633),type="l",ylab="",col="red") arrows(-(2*1.9633),-.05,-(2*1.9633),0.2,angle=0) arrows((2*1.9633),-.05,(2*1.9633),0.2,angle=0) lines(xvals,dnorm(xvals,mean=mean(d),sd=sd(d)),lwd=2) #lines(density(d),lwd=2) text(-4.5,.008,expression(alpha/2),cex=1.5,col="red") text(4.5,.008,expression(alpha/2),cex=1.5,col="red") text(0,.008,expression(1-alpha),cex=1.5,col="red") ################################################### ### chunk number 74: ################################################### plot(xvals,dnorm(xvals,mean=mean(d),sd=sd(d)),lwd=2,type="l") arrows(-(2*1.9633),-.05,-(2*1.9633),0.2,angle=0) arrows((2*1.9633),-.05,(2*1.9633),0.2,angle=0) text(-2,.08,expression(beta),cex=1.5,col="red") text(-5,.02,expression(1-beta),cex=1.5,col="black") arrows((2*1.9633+1),.04,(2*1.9633+.2),0,angle=45) text((2*1.9633+1),.04,expression(1-beta),cex=1.5,col="black") ################################################### ### chunk number 75: ################################################### numdrops <- 40 p <- .5 n <- c(0:numdrops) num <- numdrops probs <- c() for(k in n){ currentk <- binomialprobability(num,p,k) probs <- append(probs,currentk) } # proportion plot props <- n/num plot(props,probs,type="p",col="limegreen") lines(props,probs,col="limegreen",lwd=2) # draw 95% confidence intervals segments(props[[(numdrops/2+1)-5]],-0.5, props[[(numdrops/2+1)-5]],.06,col="limegreen",lty=1,lwd=2) segments(props[[(numdrops/2+1)+5]],-0.5, props[[(numdrops/2+1)+5]],.06,col="limegreen",lty=1,lwd=2) numdrops <- 400 # was 4 p <- .5 # was .5 n <- c(0:numdrops) num <- numdrops probs <- c() for(k in n){ currentk <- binomialprobability(num,p,k) probs <- append(probs,currentk) } # proportion plot props <- n/num #plot(props,probs,new=TRUE,type="p",col="red") lines(props,probs,col="red",lwd=2,lty=2) # draw 95% confidence intervals segments(props[[(numdrops/2+1)-5]],-0.5, props[[(numdrops/2+1)-5]],.06,col="red",lty=2,lwd=2) segments(props[[(numdrops/2+1)+5]],-0.5, props[[(numdrops/2+1)+5]],.06,col="red",lty=2,lwd=2) leg.txt <- c("Sample of 40","Sample of 400") legend(1,.125, legend = leg.txt, col=c("limegreen","red"), lty=c(1,2), cex=1.2, #magnification of legend lwd=2, xjust=1, yjust=1, merge=TRUE) ################################################### ### chunk number 76: ################################################### TOST <- function(mean1,mean2,theta,n1,n2,sigma){ d <- (mean2 - mean1) t1 <- (d-theta)/(sigma*(sqrt((1/n1)+(1/n2)))) t2 <- (d+theta)/(sigma*(sqrt((1/n1)+(1/n2)))) tcrit <- qt(.95,(n1+n2 - 2)) if((t1 < -tcrit) && (t2 > tcrit)){print(t1) print(t2) print(tcrit) print(c("Equivalent"))} else{print(c("Failed to show equivalence"))}} ################################################### ### chunk number 77: ################################################### d1 <- rep(NA,1000) for(i in c(1:1000)){ sample1 <- rnorm(11,mean=60,sd=4) sample2 <- rnorm(15,mean=63.5,sd=6) d1[i] <- mean(sample1) - mean(sample2) } d2 <- rep(NA,1000) for(i in c(1:1000)){ sample1 <- rnorm(11,mean=60,sd=4) sample2 <- rnorm(15,mean=61,sd=6) d2[i] <- mean(sample1) - mean(sample2) } op <- par(mfrow=c(1,2),pty = "s") plot(density(d1),xlab="",main="Smaller p-value, larger obs. power") arrows(-(2*1.9633),-.05,-(2*1.9633),0.2,angle=0) arrows((2*1.9633),-.05,(2*1.9633),0.2,angle=0) plot(density(d2),xlab="",main="Larger p-value, smaller obs. power") arrows(-(2*1.9633),-.05,-(2*1.9633),0.2,angle=0) arrows((2*1.9633),-.05,(2*1.9633),0.2,angle=0) ################################################### ### chunk number 78: ################################################### power.t.test(n=NULL, delta=200, sd=200, sig.level=0.05, power=.80, type=c("two.sample"), alternative=c("two.sided")) ################################################### ### chunk number 79: ################################################### error <- rep(NA,1000) for(i in c(1:1000)){ sample <- rnorm(10,mean=0,sd=1) currentmean <- mean(sample) error[i] <- currentmean - 0 } plot(density(error)) ################################################### ### chunk number 80: ################################################### x <- seq(c(1:100,by=.005)) plot(density(rf(10000,2,6)), xlim=range(0,5),xlab="",main="An F-distribution with parameters 2 and 6.") ################################################### ### chunk number 81: ################################################### scores <- c(9,1,2,10,2,6,1,5,0) subj <- paste("s",rep(c(1:9),1),sep="") group <- paste("g",rep(c(1:3),1,each=3),sep="") data1 <- data.frame(scores,group,subj) ################################################### ### chunk number 82: ################################################### g1data1 <- subset(data1,group=="g1")$scores g2data1 <- subset(data1,group=="g2")$scores g3data1 <- subset(data1,group=="g3")$scores SSwithin <- sum((mean(g1data1)-g1data1)^2)+sum((mean(g2data1)-g2data1)^2)+sum((mean(g3data1)-g3data1)^2) Dfwithin <- 9 - 3 (MSwithin <- SSwithin/Dfwithin) grandmean <- mean(data1$scores) SSbetween <- 3*(mean(g1data1)-grandmean)^2+3*(mean(g2data1)-grandmean)^2+3*(mean(g3data1)-grandmean)^2 Dfbetween <- 3 - 1 (MSbetween <- SSbetween/Dfbetween) ################################################### ### chunk number 83: ################################################### (MSbetween/MSwithin) ################################################### ### chunk number 84: ################################################### aov.fm <- aov(scores ~ group+Error(subj),data1) summary(aov.fm) ################################################### ### chunk number 85: ################################################### (aov.fm) ################################################### ### chunk number 86: ################################################### lm.fm <- lm(scores~group,data=data1) (mm.fm <- model.matrix(lm.fm)) (cf.fm <- coefficients(lm.fm)) (res.fm <- residuals(lm.fm)) (cbind(mm.fm,res.fm)) ################################################### ### chunk number 87: ################################################### (data1$scores) ################################################### ### chunk number 88: ################################################### (mm.fm[,1]) ################################################### ### chunk number 89: ################################################### (mm.fm[,2:3]) ################################################### ### chunk number 90: ################################################### (mm.fm[,1]) ################################################### ### chunk number 91: ################################################### (cf.fm[2]) ################################################### ### chunk number 92: ################################################### (cf.fm[3]) ################################################### ### chunk number 93: ################################################### (mm.fm)[,2:3] ################################################### ### chunk number 94: ################################################### (anova(lm.fm)) summary(aov.fm) ################################################### ### chunk number 95: ################################################### scores <- c(3,4,5,7,6,5,1,2,3) subj <- paste("s",rep(c(1:9),1),sep="") group <- paste("g",rep(c(1:3),1,each=3),sep="") data2 <- data.frame(scores,group,subj) ################################################### ### chunk number 96: ################################################### pop1 <- rnorm(1000,mean=60,sd=4) pop2 <- rnorm(1000,mean=62,sd=4) pop3 <- rnorm(1000,mean=64,sd=4) variances <- rep(NA,1000) for(i in c(1:1000)){ s1 <- sample(pop1,11) s2 <- sample(pop2,11) s3 <- sample(pop3,11) means1 <- mean(s1) means2 <- mean(s2) means3 <- mean(s3) variances[i] <- var(c(means1,means2,means3)) } meanvar <- mean(variances) plot(density(variances),main="",xlab="") ################################################### ### chunk number 97: ################################################### pop1 <- rnorm(1000,mean=60,sd=4) pop2 <- rnorm(1000,mean=60,sd=4) pop3 <- rnorm(1000,mean=60,sd=4) ss <- function(sample){ m <- rep(mean(sample),length(sample)) m2 <- (sample - m)^2 result <- sum(m2) result } mswithin <- function(s1,s2,s3){ N <- sum(length(s1),length(s2),length(s3)) DF <- N - 3 msw <- (ss(s1)+ss(s2)+ss(s3))/DF msw } ################################################### ### chunk number 98: ################################################### mswithins <- rep(NA,1000) for(i in c(1:1000)){ s1 <- sample(pop1,11) s2 <- sample(pop2,15) s3 <- sample(pop3,20) mswithins[i] <- mswithin(s1,s2,s3) } plot(density(mswithins)) ################################################### ### chunk number 99: ################################################### msbetween <- function(s1,s2,s3){ gm <- mean(c(s1,s2,s3)) m1 <- mean(s1) m2 <- mean(s2) m3 <- mean(s3) msb <- (length(s1)*(m1 - gm)^2 + length(s2)*(m2 - gm)^2 + length(s3)*(m3 - gm)^2)/2 return(msb) } msbetweens <- rep(NA,1000) for(i in c(1:1000)){ s1 <- sample(pop1,11) s2 <- sample(pop2,15) s3 <- sample(pop3,20) msbetweens[i] <- msbetween(s1,s2,s3) } plot(density(msbetweens)) ################################################### ### chunk number 100: ################################################### pop1 <- rnorm(1000,mean=60,sd=4) pop2 <- rnorm(1000,mean=62,sd=4) pop3 <- rnorm(1000,mean=64,sd=4) mswithins <- rep(NA,1000) for(i in c(1:1000)){ s1 <- sample(pop1,11) s2 <- sample(pop2,15) s3 <- sample(pop3,20) mswithins[i] <- mswithin(s1,s2,s3) } plot(density(mswithins)) ################################################### ### chunk number 101: ################################################### pop1 <- rnorm(1000,mean=60,sd=4) pop2 <- rnorm(1000,mean=62,sd=4) pop3 <- rnorm(1000,mean=64,sd=4) msbetweens <- rep(NA,1000) for(i in c(1:1000)){ s1 <- sample(pop1,11) s2 <- sample(pop2,15) s3 <- sample(pop3,20) msbetweens[i] <- msbetween(s1,s2,s3) } plot(density(msbetweens)) ################################################### ### chunk number 102: ################################################### Fratio <- function(msbetween,mswithin){ Fvalue <- msbetween/mswithin return(Fvalue) } pop1 <- rnorm(1000,mean=60,sd=4) pop2 <- rnorm(1000,mean=60,sd=4) pop3 <- rnorm(1000,mean=60,sd=4) Fs <- rep(NA,1000) for(i in c(1:1000)){ s1 <- sample(pop1,11) s2 <- sample(pop2,15) s3 <- sample(pop3,20) Fs[i] <- msbetween(s1,s2,s3)/mswithin(s1,s2,s3) } plot(density(Fs),xlim=range(0,8)) ################################################### ### chunk number 103: ################################################### pop1 <- rnorm(1000,mean=60,sd=2) pop2 <- rnorm(1000,mean=60,sd=5) pop3 <- rnorm(1000,mean=60,sd=10) Fs <- rep(NA,1000) for(i in c(1:1000)){ s1 <- sample(pop1,11) s2 <- sample(pop2,15) s3 <- sample(pop3,20) Fs[i] <- msbetween(s1,s2,s3)/mswithin(s1,s2,s3) } plot(density(Fs)) ################################################### ### chunk number 104: ################################################### library(faraway) data(stat500) attach(stat500) x <- mean(midterm) sdx <- sd(midterm) y <- mean(final) sdy <- sd(final) # three ways to see normality: #createPS("diststat500.ps") op <- par(mfrow=c(3,2)) hist(final) hist(midterm) #dev.off() #createPS("boxplotstat500.ps") #op <- par(mfrow=c(2,1)) boxplot(final) boxplot(midterm) #dev.off() #createPS("qqstat500.ps") #op <- par(mfrow=c(2,1)) qqnorm(final) qqnorm(midterm) ################################################### ### chunk number 105: ################################################### # bimodality and what it means: sample1 <- rnorm(1000,mean=1,sd=2) sample10 <- rnorm(1000,mean=10,sd=2) op <- par(mfrow=c(3,1)) #createPS("bimodalhist.ps") hist(append(sample1,sample10)) #dev.off() #createPS("bimodalboxplot.ps") boxplot(append(sample1,sample10)) #dev.off() #createPS("bimodalqq.ps") qqnorm(append(sample1,sample10)) #dev.off() ################################################### ### chunk number 106: ################################################### sample50 <- rnorm(1000,mean=50,sd=2) sample20 <- rnorm(1000,mean=20,sd=2) #createPS("1-10-20.ps") op <- par(mfrow=c(2,1)) hist(append(append(sample1,sample10),sample20)) qqnorm(append(append(sample1,sample10),sample20)) #dev.off() ################################################### ### chunk number 107: ################################################### op <- par(mfrow=c(2,1)) hist(append(append(sample1,sample10),sample50)) qqnorm(append(append(sample1,sample10),sample50)) ################################################### ### chunk number 108: ################################################### # a "perfect" relationship: final2 <- final #createPS("trivial.ps") plot(final ~ final2,xlab="final") #dev.off() ################################################### ### chunk number 109: ################################################### final2 <- final #createPS("trivialwithline.ps") plot(final ~ final2,xlab="final") abline(0,1,col="red") #dev.off() ################################################### ### chunk number 110: ################################################### #createPS("finalmidterm.ps") plot(final ~ midterm) abline(0,1,col="red") #dev.off() ################################################### ### chunk number 111: ################################################### #createPS("centredraw.ps") plot(final ~ midterm) arrows(x,min(final), x,max(final),code=0) arrows(min(midterm),y, max(midterm),y,code=0) #dev.off() detach(stat500) ################################################### ### chunk number 112: ################################################### #createPS("centredz2.ps") scaledstat500 <- data.frame(scale(stat500)) attach(scaledstat500) plot(final ~ midterm) arrows(mean(midterm),min(final), mean(midterm),max(final),code=0) arrows(min(midterm),mean(final), max(midterm),mean(final),code=0) text(1,2,labels=expression(x[i]%*%y[i]),cex=1.2,col="green") text(1.5,2,labels=c("= +ve"),cex=1.2,col="green") text(-1,-2,labels=expression(x[i]%*%y[i]),cex=1.2,col="green") text(-.5,-2,labels=c("= +ve"),cex=1.2,col="green") text(1,-2,labels=expression(x[i]%*%y[i]),cex=1.2,col="red") text(1.5,-2,labels=c("= -ve"),cex=1.2,col="red") text(-1,-2,labels=expression(x[i]%*%y[i]),cex=1.2,col="green") text(-.5,-2,labels=c("= +ve"),cex=1.2,col="green") text(-1,2,labels=expression(x[i]%*%y[i]),cex=1.2,col="red") text(-.5,2,labels=c("= -ve"),cex=1.2,col="red") #dev.off() ################################################### ### chunk number 113: ################################################### # home-made correlation r: sum(final*midterm)/(length(final)-1) # r-provided correlation: cor(midterm,final) ################################################### ### chunk number 114: ################################################### library(UsingR) data(galton) attach(galton) gx <- mean(parent) gsdx <- sd(parent) gy <- mean(child) gsdy <- sd(child) #createPS("galtondata.ps") plot(child ~ parent) arrows(gx,min(child), gx,max(child),code=0) arrows(min(parent),gy, max(parent),gy,code=0) #dev.off() detach(galton) ################################################### ### chunk number 115: ################################################### # more subtle question: #createPS("oneSDabove.ps") plot(final ~ midterm) arrows(mean(midterm),min(final), mean(midterm),max(final),code=0) arrows(min(midterm),mean(final), max(midterm),mean(final),code=0) arrows(1-.5/sdx,min(final), 1-.5/sdx,max(final),code=0) arrows(1+.5/sdx,min(final), 1+.5/sdx,max(final),code=0) #dev.off() ################################################### ### chunk number 116: ################################################### oneSDsubsampleabove <- subset(subset(scaledstat500,(1 -.5/sdx) < midterm), (1+.5/sdx) > midterm) yoneSDabove <- mean(oneSDsubsampleabove$final) ################################################### ### chunk number 117: ################################################### #createPS("oneSDbelow.ps") plot(final ~ midterm) arrows(mean(midterm),min(final), mean(midterm),max(final),code=0) arrows(min(midterm),mean(final), max(midterm),mean(final),code=0) arrows(-1-.5/sdx,min(final), -1-.5/sdx,max(final),code=0) arrows(-1+.5/sdx,min(final), -1+.5/sdx,max(final),code=0) #dev.off() ################################################### ### chunk number 118: ################################################### oneSDsubsamplebelow <- subset( subset(scaledstat500, (-1 -.5/sdx) < midterm), (-1+.5/sdx) > midterm) oneSDsubsamplebelow yoneSDbelow <- mean(oneSDsubsamplebelow$final) yoneSDbelow ################################################### ### chunk number 119: ################################################### #createPS("ratios.ps") plot(final ~ midterm) arrows(mean(midterm),min(final), mean(midterm),max(final),code=0) arrows(min(midterm),mean(final), max(midterm),mean(final),code=0) abline(0,1,col="red") abline(0,.5452,col="green") text(1.5,2,labels=c("1:1 ratio of change"),col="red") text(1.45,.3,labels=c("0.5:1 ratio of change"),col="green") #dev.off() ################################################### ### chunk number 120: ################################################### detach(scaledstat500) ################################################### ### chunk number 121: ################################################### attach(stat500) # fitting by least squares lm.stat500 <- lm(final ~ midterm) #createPS("finalmidtermreg.ps") plot(final ~ midterm) abline(lm.stat500,col="red") text(15,24,labels=c("y = 15.0462 + 0.5633x"),cex=1.2,col="red") #dev.off() ################################################### ### chunk number 122: ################################################### anova(lm.stat500) ################################################### ### chunk number 123: ################################################### summary(aov(final~midterm,stat500)) ################################################### ### chunk number 124: ################################################### summary(lm.stat500) ################################################### ### chunk number 125: ################################################### data1<-c( 49,47,46,47,48,47,41,46,43,47,46,45, 48,46,47,45,49,44,44,45,42,45,45,40, 49,46,47,45,49,45,41,43,44,46,45,40, 45,43,44,45,48,46,40,45,40,45,47,40) # across subjects then conditions Hays.mul.df <- as.data.frame(matrix(data1, ncol= 4, dimnames = list(paste("subj", 1:12), c("Shape1.Color1", "Shape2.Color1", "Shape1.Color2", "Shape2.Color2")))) Hays.df <- data.frame(rt = data1, subj = factor(rep(paste("subj", 1:12, sep=""), 4)), shape = factor(rep(rep(c("shape1", "shape2"), c(12, 12)), 2)), color = factor(rep(c("color1", "color2"), c(24, 24)))) ################################################### ### chunk number 126: ################################################### #write(t(Hays.df),file="hays.txt",ncolumns=4) anova.fm <- aov(rt ~ shape * color + Error(subj/(shape * color)), data=Hays.df) summary(anova.fm) coefficients(anova.fm) ################################################### ### chunk number 127: ################################################### # The error strata: # Regular SS-within: c1 <- Hays.mul.df$Shape1.Color1 c2 <- Hays.mul.df$Shape1.Color2 c3 <- Hays.mul.df$Shape2.Color1 c4 <- Hays.mul.df$Shape2.Color2 # sum((mean(c1)-c1)^2)+sum((mean(c2)-c2)^2)+sum((mean(c3)-c3)^2)+sum((mean(c4)-c4)^2)# 284 ################################################### ### chunk number 128: ################################################### #(colmeans-Hays.mul.df)^2 # Shape2 - Shape1 difference within subjects Shape1 <- (c1+c2)*.5 Shape2 <- (c3+c4)*.5 DShape <- Shape2 - Shape1 SumSqShape <- sum((mean(DShape) - DShape)^2) # 17.5 sdShape <- sd(DShape) n <- 12 barD <- mean(DShape) (F <- (n*(barD^2))/(sdShape^2)) # 7.542857 ################################################### ### chunk number 129: ################################################### # Color2 - Color1 difference within subjects Color1 <- (c1+c3)*.5 Color2 <- (c2+c4)*.5 DColor <- Color2 - Color1 sum((mean(DColor) - DColor)^2) # 9.5 sdColor <- sd(DColor) n <- 12 barD <- mean(DColor) (F <- (n*(barD^2))/(sdColor^2)) # 13.89474 ################################################### ### chunk number 130: ################################################### # Color-Shape interaction by subject Shapes <- (c1-c2)*.5 Colors <- (c3-c4)*.5 DSC <- (Shapes - Colors) sum((mean(DSC)-DSC)^2) # 30.5 sdDSC <- sd(DSC) n <- 12 barD <- mean(DSC) (F <- (n*(barD^2))/(sdDSC^2)) # 0 # Subjects SS: #SubjectRTs <- (c1 + c2 + c3 + c4)*.5 #sum((mean(SubjectRTs) - SubjectRTs)^2) # 226.5 ################################################### ### chunk number 131: ################################################### # Create a [row x col] multiplot multiplot <- function(row,col){ op <- par(mfrow=c(row,col),pty="s") } #Lattice plots: some sensible defaults for scales: scalelist <- list(x=list(alternating=1), y=list(alternating=1), tck=c(.5)) #function for plotting the regression lines: drawfittedline <- function(x,y){ panel.xyplot(x,y) panel.lmline(x,y,type="l",lwd=1,col="white")} ################################################### ### chunk number 132: ################################################### library(lme4) MathAchieve<-read.table("mathachieve.txt") colnames(MathAchieve) <- c("School","Minority","Sex","SES","MathAch","MEANSES") head(MathAchieve) ################################################### ### chunk number 133: ################################################### lm0 <- lm(MathAch ~ SES, data=MathAchieve) summary(lm0) coefficients(lm0) ################################################### ### chunk number 134: ################################################### population <- rexp(1000) plot(density(population)) ################################################### ### chunk number 135: ################################################### sample.size <- 100 n.sim <- 1000 means <- rep(NA,n.sim) for(i in c(1:n.sim)){ sample100 <- sample(population,sample.size) means[i] <- mean(sample100) } plot(density(means)) ################################################### ### chunk number 136: ################################################### lm1 <- lm(MathAch ~ SES, data=subset(MathAchieve,School==1224)) summary(lm1) (lm1$coefficients) ################################################### ### chunk number 137: ################################################### lm2 <- lm(MathAch ~ SES, data=subset(MathAchieve,School==1288)) summary(lm2) (lm2$coefficients) ################################################### ### chunk number 138: ################################################### multiplot(1,2) plot(MathAch ~ SES,data=subset(MathAchieve,School==1224),main="School 1224") abline(lm1$coefficients) plot(MathAch ~ SES,data=subset(MathAchieve,School==1288),main="School 1228") abline(lm2$coefficients) ################################################### ### chunk number 139: ################################################### library(lattice) (print(xyplot(MathAch~SES|factor(School), MathAchieve, xlab="Student SES", ylab="Math Achievement", panel=drawfittedline, scales=scalelist))) ################################################### ### chunk number 140: ################################################### lme1 <- lmList(MathAch ~ SES|School, MathAchieve) ## equivalent way (i.e., the intercept is implicit in the command above): lme1 <- lmList(MathAch ~ 1+SES|School, MathAchieve) ## Aside: you can drop the intercept as follows ## lme1 <- lmList(MathAch ~ SES-1|School, MathAchieve) lme1$"1224" lme1$"1288" ################################################### ### chunk number 141: ################################################### t.test(coef(lme1)[1]) t.test(coef(lme1)[2]) ################################################### ### chunk number 142: ################################################### lme1 <- lmList(MathAch ~ SES|School, MathAchieve) intercepts <- coef(lme1)[1] slopes <- coef(lme1)[2] (cov(intercepts,slopes)) (cov(intercepts,slopes)/sqrt(var(intercepts)*var(slopes))) (cor(intercepts,slopes)) ################################################### ### chunk number 143: ################################################### intslopes<-data.frame(intercepts,slopes) colnames(intslopes) <- c("Intercepts","Slopes") plot(Intercepts~Slopes,intslopes) lm.intslopes <- lm(Intercepts~Slopes,data=intslopes) abline(coefficients(lm.intslopes)) ################################################### ### chunk number 144: ################################################### MathAchSchool <- read.table("mathachschool.txt") colnames(MathAchSchool) <- c("School", "Size","Sector", "PRACAD","DISCLIM","HIMINTY","MEANSES") MathScores <- merge(MathAchieve,MathAchSchool,by="School") ################################################### ### chunk number 145: ################################################### (contrasts(MathScores$Sector)) ################################################### ### chunk number 146: ################################################### #lmepub <- lmList(MathAch ~ SES|School, subset(MathScores,Sector=="Public")) #slopespub <- coef(lmepub)[2] #mean(slopespub) ################################################### ### chunk number 147: ################################################### #lmecath <- lmList(MathAch ~ SES|School, subset(MathScores,Sector=="Catholic")) #slopescath <- coef(lmecath)[2] #mean(slopescath) - mean(slopespub) ################################################### ### chunk number 148: ################################################### #library(lme4) #lme2 <- lmList(MathAch ~ SES|School, subset(MathScores,Sector=="Public")) #intercepts <- coef(lme2)[1] #slopes <- coef(lme2)[2] #intslopes<-data.frame(intercepts,slopes) #colnames(intslopes) <- c("Intercepts","Slopes") #plot(Intercepts~Slopes,intslopes) #lm.intslopes <- lm(Intercepts~Slopes,data=intslopes) #abline(coefficients(lm.intslopes)) ################################################### ### chunk number 149: ################################################### #lme2 <- lmList(MathAch ~ SES|School, subset(MathScores,Sector=="Catholic")) #intercepts <- coef(lme2)[1] #slopes <- coef(lme2)[2] #intslopes<-data.frame(intercepts,slopes) #colnames(intslopes) <- c("Intercepts","Slopes") #plot(Intercepts~Slopes,intslopes) #lm.intslopes <- lm(Intercepts~Slopes,data=intslopes) #abline(coefficients(lm.intslopes)) ################################################### ### chunk number 150: ################################################### lme1.fm <- lmer(MathAch~SES+Sector+(1+SES|School),MathScores) summary(lme1.fm) #library(lme4) #lmer1.fm <- lmer(MathAch~SES+Sector+(SES|School),MathScores) #summary(lmer1.fm) ################################################### ### chunk number 151: ################################################### fixef(lme1.fm) ################################################### ### chunk number 152: ################################################### females <- rbinom(40,1,.4) females sum(females) ################################################### ### chunk number 153: ################################################### populationsize <- 1000 samplesize <- 40 p <- .4 # probability of females compute95CIpopulation <- function(populationsize,samplesize,p){ females <- rbinom(samplesize,1,.4) samplesumfemales <- sum(females) sdsample <- sqrt(samplesize*p*(1-p)) sample95CI <- c(samplesumfemales - 2*sdsample, samplesumfemales + 2*sdsample) population95CI <- populationsize/samplesize * sample95CI print(population95CI)} ################################################### ### chunk number 154: ################################################### compute95CIpopulation(1000,40,.4) compute95CIpopulation(1000,40,.4) compute95CIpopulation(1000,40,.4) compute95CIpopulation(1000,40,.4) compute95CIpopulation(1000,40,.4) compute95CIpopulation(1000,40,.4) compute95CIpopulation(1000,40,.4) compute95CIpopulation(1000,40,.4) compute95CIpopulation(1000,40,.4) compute95CIpopulation(1000,40,.4) compute95CIpopulation(1000,40,.4) compute95CIpopulation(1000,40,.4) compute95CIpopulation(1000,40,.4) compute95CIpopulation(1000,40,.4)