1 WEIGHTED <- TRUE # weighted analysis?
   2 
   3 d <- read.table('../data_sample/merged/mil.csv', header=F, sep=' ', as.is = T,
   4                 col.names=c('aid', 'byr', 'geno', 'cand', 'test',
   5                   'top', 'off', 'sta', 'nd', 'nh',
   6                   'Cf', 'edc', 'rel', 'ebv',
   7                   'Df', 'edcd', 'reld', 'dpgm',
   8                   'Cr', 'edcr', 'relr', 'ebvr',
   9                   'Gr', 'edcg', 'relg', 'gebv')
  10                 )
  11 d <- d[d$test=='Y', c('edc','dpgm','gebv','ebvr')]
  12 cat('rows=',nrow(d))
  13 attach(d)
  14 if (WEIGHTED) {
  15   reg1 <- lm(dpgm~gebv,weight=edc)
  16   reg2 <- lm(dpgm~ebvr,weight=edc)
  17   reg3 <- lm(dpgm~ebvr+gebv,weight=edc)
  18 } else {
  19   reg1 <- lm(dpgm~gebv)
  20   reg2 <- lm(dpgm~ebvr)
  21   reg3 <- lm(dpgm~ebvr+gebv)
  22 }
  23 print(reg1)
  24 print(reg2)
  25 print(reg3)
  26 # a1 <- reg1$coefficients[1]; b1 <- reg1$coefficients[2]
  27 # a2 <- reg2$coefficients[1]; b2 <- reg2$coefficients[2]
  28 # a3 <- reg3$coefficients[1]; b3 <- reg3$coefficients[2];
  29 # c3 <- reg3$coefficients[3]
  30 R1 <- summary(reg1)$r.squared # R^2 for GEBV
  31 R2 <- summary(reg2)$r.squared # R^2 for EBVr
  32 R3 <- summary(reg3)$r.squared # R^2 for both
  33 cat(sprintf('R2_1=%5.3f  R2_2=%5.3f  R2_3=%5.3f\n',R1,R2,R3))
  34 
  35 ##== WEIGHTED = TRUE ========================================
  36 # [prog]$ R -q --slave <test.R
  37 # rows= 164
  38 #
  39 # lm(formula = dpgm ~ gebv, weights = edc)
  40 # (Intercept)         gebv
  41 #   -548.0635       0.9945
  42 #
  43 # lm(formula = dpgm ~ ebvr, weights = edc)
  44 # (Intercept)         ebvr
  45 #   -357.7728       0.6717
  46 #
  47 # lm(formula = dpgm ~ ebvr + gebv, weights = edc)
  48 # (Intercept)         ebvr         gebv
  49 #   -417.1744      -0.2541       1.1688
  50 #
  51 # R2_1=0.619  R2_2=0.245  R2_3=0.635
  52 #
  53 ##== WEIGHTED = FALSE =======================================
  54 # [prog]$ R -q --slave <test.R
  55 # rows= 164
  56 #
  57 # lm(formula = dpgm ~ gebv)
  58 # (Intercept)         gebv
  59 #    -559.789        1.024
  60 #
  61 # lm(formula = dpgm ~ ebvr)
  62 # (Intercept)         ebvr
  63 #   -445.7595       0.7709
  64 #
  65 # lm(formula = dpgm ~ ebvr + gebv)
  66 # (Intercept)         ebvr         gebv
  67 #   -499.5708      -0.1114       1.0991
  68 #
  69 # R2_1=0.623  R2_2=0.295  R2_3=0.626
  70 #