Nov <- c(23.9,43.3,36.3,40.6,57.0,52.5,46.1,142.0,112.6,23.7)
Dec <- c(41.0,52.0,18.7,55.0,40.0,29.2,51.0,17.6,46.6,57.0)
n1 <- length(Nov)
x1 <- Nov
y1 <- Dec
x1bar <- mean(x1)
y1bar <- mean(y1)
Sxx1 <- sum((x1-x1bar)^2)
Syy1 <- sum((y1-y1bar)^2)
Sxy1 <- sum((x1-x1bar)*(y1-y1bar))
See1 <- Syy1 - Sxy1^2/Sxx1
cat("n1 =",n1,"x1bar =",x1bar,"y1bar =",y1bar,"\n")
cat("Sxx1 =",Sxx1,"Syy1 =",Syy1,"Sxy1 =",Sxy1,"See1 =",See1,"\n")
a1 <- y1bar
b1 <- Sxy1/Sxx1
r1 <- cor(x1,y1)
cat("a1 =",a1,"b1 =",b1,"r1 =",r1,"\n")
NewNov <- c(34.1,62.0,106.9,34.1,68.3,81.0)
NewDec <- c(12.3,90.4,28.8,106.2,62.3,50.5)
n2 <- length(NewNov)
x2 <- NewNov
y2 <- NewDec
x2bar <- mean(x2)
y2bar <- mean(y2)
Sxx2 <- sum((x2-x2bar)^2)
Syy2 <- sum((y2-y2bar)^2)
Sxy2 <- sum((x2-x2bar)*(y2-y2bar))
See2 <- Syy2 - Sxy2^2/Sxx2
cat("n2 =",n2,"x2bar =",x2bar,"y2bar =",y2bar,"\n")
cat("Sxx2 =",Sxx2,"Syy2 =",Syy2,"Sxy2 =",Sxy2,"See2 =",See2,"\n")
a2 <- y2bar
b2 <- Sxy2/Sxx2
r2 <- cor(x2,y2)
cat("a2 =",a2,"b2 =",b2,"r2 =",r2,"\n")
x0 <- c(x1,x2)
y0 <- c(y1,y2)
n0 <- length(x0)
x0bar <- mean(x0)
y0bar <- mean(y0)
Sxx0 <- sum((x0-x0bar)^2)
Syy0 <- sum((y0-y0bar)^2)
Sxy0 <- sum((x0-x0bar)*(y0-y0bar))
See0 <- Syy0 - Sxy0^2/Sxx0
cat("n0 =",n0,"x0bar =",x0bar,"y0bar =",y0bar,"\n")
cat("Sxx0 =",Sxx0,"Syy0 =",Syy0,"Sxy0 =",Sxy0,"See0 =",See0,"\n")
a0 <- y0bar
b0 <- Sxy0/Sxx0
r0 <- cor(x0,y0)
cat("a0 =",a0,"b0 =",b0,"r0 =",r0,"\n")
n <- n1 + n2
xbar <- (n1*x1bar + n2*x2bar)/n
ybar <- (n1*y1bar + n2*y2bar)/n
nh <- (1/n1 + 1/n2)^(-1)
Sxxc <- nh*(x1bar-x2bar)^2
Sxyc <- nh*(x1bar-x2bar)*(y1bar-y2bar)
bc <- Sxyc/Sxxc
Sxx <- Sxx1 + Sxx2 + Sxxc
a <- ybar
b <- (b1*Sxx1 + b2*Sxx2 + bc*Sxxc)/Sxx
See <- See1+See2+
       ((b1-b2)^2*Sxx1*Sxx2+(b2-bc)^2*Sxx2*Sxxc+(bc-b1)^2*Sxxc*Sxx1)/Sxx
See <- See1+See2+
       ((b1-b2)^2*Sxx1*Sxx2+(b2-bc)^2*Sxx2*Sxxc+(bc-b1)^2*Sxxc*Sxx1)/Sxx
cat("\nNote calculated values n, xbar, ybar, Sxx, See, a and b\n")
cat("are same as values n0, x0bar, y0bar, Sxx0, See0, a0 and b0\n")
cat("found from combined sample\n")
cat("n  =",n,"xbar  =",xbar,"ybar  =",ybar,"\n")
cat("Sxx  =",Sxx,"See  =",See,"a  =",a,"b  =",b,"\n")
cat("n0 =",n0,"x0bar =",x0bar,"y0bar =",y0bar,"\n")
cat("Sxx0 =",Sxx0,"See0 =",See0,"a0 =",a0,"b0 =",b0,"\n")
mod1 <- lm(y1 ~ x1)
mod2 <- lm(y2 ~ x2)
mod0 <- lm(y0 ~ x0)
print(summary(mod1))
print(summary(mod2))
print(summary(mod0))