set.seed(1357) nba <- read.csv("http://stat.ufl.edu/~winner/data/nba_ht_wt.csv",header=T) attach(nba); names(nba) nba.samp <- sample(1:length(Height),12,replace=F) nba.sample <- nba[nba.samp,c(1,3,4)] nba.sample detach(nba) nhl <- read.csv("http://stat.ufl.edu/~winner/data/nhl_ht_wt.csv",header=T) attach(nhl); names(nhl) nhl.samp <- sample(1:length(Height),12,replace=F) nhl.sample <- nhl[nhl.samp,c(2,4,5)] nhl.sample detach(nhl) epl <- read.csv("http://stat.ufl.edu/~winner/data/epl_2015_ht_wt.csv",header=T) attach(epl); names(epl) epl.samp <- sample(1:length(Height),12,replace=F) epl.sample <- epl[epl.samp,c(2,6,7)] epl.sample all.sample <- rbind(nba.sample,nhl.sample,epl.sample) all.sample plot(Height,Weight) NBA <- c(rep(1,12),rep(0,24)) NHL <- c(rep(0,12),rep(1,12),rep(0,12)) League <- c(rep(1,12),rep(2,12),rep(3,12)) League <- factor(League) all.sample1 <- data.frame(all.sample,NBA,NHL,League) wtht.1 <- lm(Weight ~ Height + NBA + NHL, data=all.sample1) summary(wtht.1) anova(wtht.1) drop1(wtht.1,test="F") wtht.2 <- lm(Weight ~ Height, data=all.sample1) summary(wtht.2) anova(wtht.2) anova(wtht.2,wtht.1) ht.x <- seq(62,82,0.1) yhat.nba <- coef(wtht.1)[1] + ht.x*coef(wtht.1)[2] + coef(wtht.1)[3] yhat.nhl <- coef(wtht.1)[1] + ht.x*coef(wtht.1)[2] + coef(wtht.1)[4] yhat.epl <- coef(wtht.1)[1] + ht.x*coef(wtht.1)[2] plot(Height,Weight,pch=as.numeric(League),ylim=c(120,240)) lines(ht.x,yhat.nba,lty=1) lines(ht.x,yhat.nhl,lty=2) lines(ht.x,yhat.epl,lty=5) legend(65,240,c("NBA","NHL","EPL"),pch=c(1,2,3),lty=c(1,2,5)) wtht.3 <- lm(Weight ~ Height + NBA + NHL + I(Height*NBA) + I(Height*NHL), data=all.sample1) summary(wtht.3) anova(wtht.3) drop1(wtht.3, test="F") anova(wtht.1,wtht.3) ht.x <- seq(62,82,0.1) yhat.nba3 <- coef(wtht.3)[1] + ht.x*coef(wtht.3)[2] + coef(wtht.3)[3] + ht.x*coef(wtht.3)[5] yhat.nhl3 <- coef(wtht.3)[1] + ht.x*coef(wtht.3)[2] + coef(wtht.3)[4] + ht.x*coef(wtht.3)[6] yhat.epl3 <- coef(wtht.3)[1] + ht.x*coef(wtht.3)[2] plot(Height,Weight,pch=as.numeric(League),ylim=c(120,240)) lines(ht.x,yhat.nba3,lty=1) lines(ht.x,yhat.nhl3,lty=2) lines(ht.x,yhat.epl3,lty=5) legend(65,240,c("NBA","NHL","EPL"),pch=c(1,2,3),lty=c(1,2,5)) wtht.nba <- lm(Weight[League==1] ~ Height[League==1],data=all.sample1) summary(wtht.nba) wtht.nhl <- lm(Weight[League==2] ~ Height[League==2],data=all.sample1) summary(wtht.nhl) wtht.epl <- lm(Weight[League==3] ~ Height[League==3],data=all.sample1) summary(wtht.epl) (df.nba <- df.residual(wtht.nba)); (s2.nba <- deviance(wtht.nba)/df.nba) (df.nhl <- df.residual(wtht.nhl)); (s2.nhl <- deviance(wtht.nhl)/df.nhl) (df.epl <- df.residual(wtht.epl)); (s2.epl <- deviance(wtht.epl)/df.epl) (df.all <- df.residual(wtht.3)); (s2.all <- deviance(wtht.3)/df.all) (s2.chk <- (df.nba*s2.nba+df.nhl*s2.nhl+df.epl*s2.epl)/df.all) Bart1 <- 1/(1+(1/df.nba+1/df.nhl+1/df.epl-1/df.all)/(3*(3-1))) Bart2 <- df.all*log(s2.all) - (df.nba*log(s2.nba) + df.nhl*log(s2.nhl) + df.epl*log(s2.epl)) Bart.X2 <- Bart1*Bart2 X2.05 <- qchisq(.95,3-1) X2.p <- 1-pchisq(Bart.X2,3-1) Bart.out <- cbind(Bart.X2,3-1,X2.05,X2.p) colnames(Bart.out) <- c("X2-Stat", "DF", "X2(.05)", "P-value") round(Bart.out,4) all.sample1