WEIGHTED <- TRUE # weighted analysis?

d <- read.table('../data_sample/merged/mil.csv', header=F, sep=' ', as.is = T,
                col.names=c('aid', 'byr', 'geno', 'cand', 'test',
                  'top', 'off', 'sta', 'nd', 'nh',
                  'Cf', 'edc', 'rel', 'ebv',
                  'Df', 'edcd', 'reld', 'dpgm',
                  'Cr', 'edcr', 'relr', 'ebvr',
                  'Gr', 'edcg', 'relg', 'gebv')
                )
d <- d[d$test=='Y', c('edc','dpgm','gebv','ebvr')]
cat('rows=',nrow(d))
attach(d)
if (WEIGHTED) {
  reg1 <- lm(dpgm~gebv,weight=edc)
  reg2 <- lm(dpgm~ebvr,weight=edc)
  reg3 <- lm(dpgm~ebvr+gebv,weight=edc)
} else {
  reg1 <- lm(dpgm~gebv)
  reg2 <- lm(dpgm~ebvr)
  reg3 <- lm(dpgm~ebvr+gebv)
}
print(reg1)
print(reg2)
print(reg3)
# a1 <- reg1$coefficients[1]; b1 <- reg1$coefficients[2]
# a2 <- reg2$coefficients[1]; b2 <- reg2$coefficients[2]
# a3 <- reg3$coefficients[1]; b3 <- reg3$coefficients[2];
# c3 <- reg3$coefficients[3]
R1 <- summary(reg1)$r.squared # R^2 for GEBV
R2 <- summary(reg2)$r.squared # R^2 for EBVr
R3 <- summary(reg3)$r.squared # R^2 for both
cat(sprintf('R2_1=%5.3f  R2_2=%5.3f  R2_3=%5.3f\n',R1,R2,R3))

##== WEIGHTED = TRUE ========================================
# [prog]$ R -q --slave <test.R
# rows= 164
#
# lm(formula = dpgm ~ gebv, weights = edc)
# (Intercept)         gebv
#   -548.0635       0.9945
#
# lm(formula = dpgm ~ ebvr, weights = edc)
# (Intercept)         ebvr
#   -357.7728       0.6717
#
# lm(formula = dpgm ~ ebvr + gebv, weights = edc)
# (Intercept)         ebvr         gebv
#   -417.1744      -0.2541       1.1688
#
# R2_1=0.619  R2_2=0.245  R2_3=0.635
#
##== WEIGHTED = FALSE =======================================
# [prog]$ R -q --slave <test.R
# rows= 164
#
# lm(formula = dpgm ~ gebv)
# (Intercept)         gebv
#    -559.789        1.024
#
# lm(formula = dpgm ~ ebvr)
# (Intercept)         ebvr
#   -445.7595       0.7709
#
# lm(formula = dpgm ~ ebvr + gebv)
# (Intercept)         ebvr         gebv
#   -499.5708      -0.1114       1.0991
#
# R2_1=0.623  R2_2=0.295  R2_3=0.626
#