dat <- matrix(c( 0, 0, 0, 1, 3, 4, 15, 7, 0, 0, 0, 1, 5, 13, 11, 0, 0, 1, 1, 2, 25, 15, 6, 0, 0, 1, 3, 7, 14, 7, 4, 2, 0, 1, 7, 15, 28, 8, 2, 0, 0, 1, 7, 18, 15, 6, 0, 0, 0, 4, 10, 12, 8, 2, 0, 0, 0, 5, 11, 2, 3, 0, 0, 0, 9, 12, 10, 3, 1, 0, 0, 0), nrow=9,ncol=8,byrow=T) rownames(dat) <- 71:63 colnames(dat) <- seq(16.25,19.75,0.5) rowtots <- apply(dat,1,sum) coltots <- apply(dat,2,sum) grandtot <- sum(dat) print(dat) freq <- as.vector(dat) sta <- as.numeric(rep(rownames(dat),ncol(dat))) cub <- as.numeric(rep(colnames(dat),rep(nrow(dat),ncol(dat)))) stature <- rep(sta,freq) cubit <- rep(cub,freq) mean.stature <- sum(sta*freq)/grandtot mean.cubit <- sum(cub*freq)/grandtot var.stature <- sum((sta-mean.stature)^2*freq)/(grandtot-1) var.cubit <- sum((cub-mean.cubit)^2*freq)/(grandtot-1) cov.stature.cubit <- sum((sta-mean.stature)*(cub-mean.cubit)*freq)/(grandtot-1) cor.stature.cubit <- cor(stature,cubit) cat("Correlation between cubit and stature is",cor(stature,cubit),"\n") cat("Note while column totals are correctly given as\n") print(coltots) cat("some row totals are wrong - right values are\n") print(rowtots) cat("\nNeglecting extreme values data reduces to\n") print(dat[2:8,2:7]) freq.co <- freq freq.co[(sta<64)|(sta>70)|(cub>19.5)] <- 0 grandtot.co <- sum(freq.co) stature.co <- rep(sta,freq.co) cubit.co <- rep(cub,freq.co) mean.stature.co <- sum(sta*freq.co)/grandtot.co mean.cubit.co <- sum(cub*freq.co)/grandtot.co var.stature.co <- sum((sta-mean.stature.co)^2*freq.co)/(grandtot.co-1) var.cubit.co <- sum((cub-mean.cubit.co)^2*freq.co)/(grandtot.co-1) cov.stature.cubit.co <- sum((sta-mean.stature.co)*(cub-mean.cubit.co)*freq.co)/(grandtot.co-1) cor.stature.cubit.co <- cor(stature.co,cubit.co) cat("Corrected correlation between cubit and stature is", cor(stature.co,cubit.co),"\n") d <- as.vector(t(dat)) cub.factor <- gl(8,1,72) stat.factor <- gl(9,8,72) cub.length <- 15.25 + as.numeric(cub.factor) stat.length <- 70.5 - as.numeric(stat.factor) cub.reps <- rep(cub.length,d) stat.reps <- rep(stat.length,d) plot(cub.length,stat.length,type="n", xlab="Left cubit", ylab="Stature") text(cub.length,stat.length,d)