## com12a of october 2023
#  with stk4011, small babies etc.

data <- matrix( scan("smallbabies23_data", skip=14), ncol=11, byrow=TRUE )

idnumber <- data[ ,1]
yy0 <- data[ ,2]
age <- data[ ,3]
lwt <- data[, 4] * 0.45349 # from pounds to kg 
race <- data[ ,5]
smoke <- data[ ,6]
ptl <- data[ ,7]
ht <- data[ ,8]
ui <- data[ ,9]
ftv <- data[ ,10]
bwt <- data[ ,11] / 1000 # in kg, not gram

yy <- bwt
nn <- length(yy) 

x1 <- age
x2 <- lwt
x3 <- smoke 

hello123 <- lm(yy ~ x1 + x2 + x3)
summary(hello123)

resi123 <- hello123$res

# Coefficients:
#              Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  2.362720   0.300687   7.858 3.11e-13 ***
# x1           0.007093   0.009925   0.715   0.4757    age 
# x2           0.008862   0.003792   2.337   0.0205 *  weight before pregnant 
# x3          -0.267213   0.105802  -2.526   0.0124 *  smoking

hello1 <- lm(yy ~ x1)
summary(hello1)
sigma1 <- sqrt( sum( (hello1$res)^2) / (nn-2) ) 

hello2 <- lm(yy ~ x2)
summary(hello2)
sigma2 <- sqrt( sum( (hello2$res)^2) / (nn-2) ) 

hello3 <- lm(yy ~ x3)
summary(hello3)
sigma3 <- sqrt( sum( (hello3$res)^2) / (nn-2) ) 

print(round(c(sigma1,sigma2,sigma3),4))
#  0.7280 0.7182 0.7178

Sigma <- cov(cbind(x1,x2,x3))

cov2cor(Sigma)

## for x3, smoking: 
#  y = a + b x + noise = 3.0549 - 0.2817 x + noise 
#  as much as 0.28 kg lower with smoking mothers 
#  c.i.: -0.2817 plus-minus 1.96*0.1069
lower3 <- hello3$coef[2] - 1.69*0.1069
upper3 <- hello3$coef[2] + 1.69*0.1069
print(c(lower3,upper3)) # (-0.462, -0.101) 

# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  3.05496    0.06693  45.642  < 2e-16 ***
# x3          -0.28171    0.10697  -2.634  0.00916 **