#######################################################################
# Mean comparison and ANCOVA estimates
#######################################################################

samp  <- read.table("Diet0001.dat", row.names="ID",
   col.names=c("ID", "DISTR.1", "BLACK", "NBHISP", "GRADE", "SLFHLTH",
   "SLFWGHT", "WORKHARD", "GOODQUAL", "PHYSFIT", "PROUD",
   "LIKESLF", "ACCEPTED", "FEELLOVD", "DIET", "Y.0", "Y.1", "DISTR.2",
   "PI.TRUE"))

#######################################################################
# Mean comparison:

fit <- lm( DISTR.2 ~ DIET, data=samp)

# > print(summary(fit))
# 
# Call:
# lm(formula = DISTR.2 ~ DIET, data = samp)
# 
# Residuals:
#      Min       1Q   Median       3Q      Max 
# -0.70314 -0.36314 -0.09353  0.25647  1.91647 
# 
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)    
# (Intercept) 0.643525   0.006733  95.581  < 2e-16 ***
# DIET        0.059614   0.014931   3.993 6.61e-05 ***
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
# 
# Residual standard error: 0.4655 on 5998 degrees of freedom
# Multiple R-Squared: 0.002651,   Adjusted R-squared: 0.002484 
# F-statistic: 15.94 on 1 and 5998 DF,  p-value: 6.612e-05 


#######################################################################
# ANCOVA adjusting for DISTR.1:

fit <- lm( DISTR.2 ~ DISTR.1 + DIET, data=samp)

# > print(summary(fit))
# 
# Call:
# lm(formula = DISTR.2 ~ DISTR.1 + DIET, data = samp)
# 
# Residuals:
#      Min       1Q   Median       3Q      Max 
# -1.86469 -0.23638 -0.07043  0.18786  1.88682 
# 
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)    
# (Intercept) 0.266385   0.009193  28.978   <2e-16 ***
# DISTR.1     0.611411   0.011820  51.728   <2e-16 ***
# DIET        0.003039   0.012465   0.244    0.807    
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
# 
# Residual standard error: 0.3871 on 5997 degrees of freedom
# Multiple R-Squared: 0.3104,     Adjusted R-squared: 0.3101 
# F-statistic:  1349 on 2 and 5997 DF,  p-value: < 2.2e-16 



#######################################################################
# ANCOVA adjusting for all covariates:

fit <- lm( DISTR.2 ~ DISTR.1 + BLACK + NBHISP + GRADE + 
   SLFHLTH + SLFWGHT + WORKHARD + GOODQUAL + PHYSFIT + PROUD + 
   LIKESLF + ACCEPTED + FEELLOVD + DIET, data=samp)

# > print(summary(fit))
# 
# Call:
# lm(formula = DISTR.2 ~ DISTR.1 + BLACK + NBHISP + GRADE + SLFHLTH + 
#     SLFWGHT + WORKHARD + GOODQUAL + PHYSFIT + PROUD + LIKESLF + 
#     ACCEPTED + FEELLOVD + DIET, data = samp)
# 
# Residuals:
#      Min       1Q   Median       3Q      Max 
# -1.69102 -0.22781 -0.05734  0.18025  1.81213 
# 
# Coefficients:
#               Estimate Std. Error t value Pr(>|t|)    
# (Intercept) -0.0052087  0.0422214  -0.123 0.901821    
# DISTR.1      0.5178247  0.0125785  41.168  < 2e-16 ***
# BLACK        0.0745764  0.0120532   6.187 6.53e-10 ***
# NBHISP       0.0289905  0.0141752   2.045 0.040884 *  
# GRADE        0.0019352  0.0035672   0.543 0.587490    
# SLFHLTH      0.0196974  0.0059809   3.293 0.000996 ***
# SLFWGHT      0.0038613  0.0069990   0.552 0.581183    
# WORKHARD    -0.0121711  0.0056656  -2.148 0.031734 *  
# GOODQUAL     0.0209810  0.0098513   2.130 0.033231 *  
# PHYSFIT     -0.0005642  0.0069283  -0.081 0.935099    
# PROUD        0.0376379  0.0101757   3.699 0.000219 ***
# LIKESLF      0.0242944  0.0064040   3.794 0.000150 ***
# ACCEPTED     0.0167152  0.0068032   2.457 0.014040 *  
# FEELLOVD     0.0389178  0.0085254   4.565 5.10e-06 ***
# DIET        -0.0136658  0.0129321  -1.057 0.290675    
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
# 
# Residual standard error: 0.3758 on 5985 degrees of freedom
# Multiple R-Squared: 0.3513,     Adjusted R-squared: 0.3498 
# F-statistic: 231.5 on 14 and 5985 DF,  p-value: < 2.2e-16 



#######################################################################
# ANCOVA with all baseline-by-treatment interactions

# center the covariates at overall means
sampc <- samp
sampc$DISTR.1 <- sampc$DISTR.1 - mean(sampc$DISTR.1)
sampc$BLACK <- sampc$BLACK - mean(sampc$BLACK)
sampc$NBHISP <- sampc$NBHISP - mean(sampc$NBHISP)
sampc$GRADE <- sampc$GRADE - mean(sampc$GRADE)
sampc$SLFHLTH <- sampc$SLFHLTH - mean(sampc$SLFHLTH)
sampc$SLFWGHT <- sampc$SLFWGHT - mean(sampc$SLFWGHT)
sampc$WORKHARD <- sampc$WORKHARD - mean(sampc$WORKHARD)
sampc$GOODQUAL <- sampc$GOODQUAL - mean(sampc$GOODQUAL)
sampc$PHYSFIT <- sampc$PHYSFIT - mean(sampc$PHYSFIT)
sampc$PROUD <- sampc$PROUD - mean(sampc$PROUD)
sampc$LIKESLF <- sampc$LIKESLF - mean(sampc$LIKESLF)
sampc$ACCEPTED <- sampc$ACCEPTED - mean(sampc$ACCEPTED)
sampc$FEELLOVD <- sampc$FEELLOVD - mean(sampc$FEELLOVD)

# fit the model to estimate ACE
fit <- glm( DISTR.2 ~  DIET + 
   DISTR.1 + BLACK + NBHISP + GRADE + 
   SLFHLTH + SLFWGHT + WORKHARD + GOODQUAL + PHYSFIT + PROUD + 
   LIKESLF + ACCEPTED + FEELLOVD +
   DISTR.1:DIET + BLACK:DIET + NBHISP:DIET + GRADE:DIET + 
   SLFHLTH:DIET + SLFWGHT:DIET + WORKHARD:DIET + GOODQUAL:DIET +
   PHYSFIT:DIET + PROUD:DIET + 
   LIKESLF:DIET + ACCEPTED:DIET + FEELLOVD:DIET,
   data=sampc)

# > print(summary(fit))
# 
# Call:
# glm(formula = DISTR.2 ~ DIET + DISTR.1 + BLACK + NBHISP + GRADE + 
#     SLFHLTH + SLFWGHT + WORKHARD + GOODQUAL + PHYSFIT + PROUD + 
#     LIKESLF + ACCEPTED + FEELLOVD + DISTR.1:DIET + BLACK:DIET + 
#     NBHISP:DIET + GRADE:DIET + SLFHLTH:DIET + SLFWGHT:DIET + 
#     WORKHARD:DIET + GOODQUAL:DIET + PHYSFIT:DIET + PROUD:DIET + 
#     LIKESLF:DIET + ACCEPTED:DIET + FEELLOVD:DIET, data = sampc)
# 
# Deviance Residuals: 
#      Min        1Q    Median        3Q       Max  
# -1.72288  -0.22801  -0.05723   0.18175   1.82660  
# 
# Coefficients:
#                 Estimate Std. Error t value Pr(>|t|)    
# (Intercept)    0.6587668  0.0055344 119.031  < 2e-16 ***
# DIET          -0.0058574  0.0149019  -0.393 0.694286    
# DISTR.1        0.5210525  0.0143816  36.230  < 2e-16 ***
# BLACK          0.0742012  0.0131869   5.627 1.92e-08 ***
# NBHISP         0.0241334  0.0159814   1.510 0.131071    
# GRADE          0.0021585  0.0039652   0.544 0.586210    
# SLFHLTH        0.0257745  0.0067034   3.845 0.000122 ***
# SLFWGHT        0.0026900  0.0077575   0.347 0.728781    
# WORKHARD      -0.0194426  0.0062698  -3.101 0.001938 ** 
# GOODQUAL       0.0154691  0.0110400   1.401 0.161212    
# PHYSFIT        0.0054135  0.0077812   0.696 0.486635    
# PROUD          0.0329261  0.0113575   2.899 0.003756 ** 
# LIKESLF        0.0225898  0.0072863   3.100 0.001942 ** 
# ACCEPTED       0.0159821  0.0077195   2.070 0.038462 *  
# FEELLOVD       0.0412613  0.0095958   4.300 1.74e-05 ***
# DIET:DISTR.1  -0.0106838  0.0296963  -0.360 0.719033    
# DIET:BLACK    -0.0003361  0.0327604  -0.010 0.991814    
# DIET:NBHISP    0.0238344  0.0346707   0.687 0.491826    
# DIET:GRADE     0.0023927  0.0091386   0.262 0.793465    
# DIET:SLFHLTH  -0.0297672  0.0150042  -1.984 0.047310 *  
# DIET:SLFWGHT  -0.0011506  0.0181310  -0.063 0.949404    
# DIET:WORKHARD  0.0392602  0.0147076   2.669 0.007620 ** 
# DIET:GOODQUAL  0.0260125  0.0246058   1.057 0.290477    
# DIET:PHYSFIT  -0.0266321  0.0171818  -1.550 0.121190    
# DIET:PROUD     0.0240841  0.0255893   0.941 0.346650    
# DIET:LIKESLF   0.0086665  0.0153017   0.566 0.571161    
# DIET:ACCEPTED -0.0015125  0.0164325  -0.092 0.926664    
# DIET:FEELLOVD -0.0125126  0.0209234  -0.598 0.549850    
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
# 
# (Dispersion parameter for gaussian family taken to be 0.1410809)
# 
#     Null deviance: 1303.09  on 5999  degrees of freedom
# Residual deviance:  842.54  on 5972  degrees of freedom
# AIC: 5306.7
# 
# Number of Fisher Scoring iterations: 2



# center the covariates at means among the treated
sampt <- samp
sampt$DISTR.1 <- sampt$DISTR.1 - mean(sampt$DISTR.1[sampt$DIET==1])
sampt$BLACK <- sampt$BLACK - mean(sampt$BLACK[sampt$DIET==1])
sampt$NBHISP <- sampt$NBHISP - mean(sampt$NBHISP[sampt$DIET==1])
sampt$GRADE <- sampt$GRADE - mean(sampt$GRADE[sampt$DIET==1])
sampt$SLFHLTH <- sampt$SLFHLTH - mean(sampt$SLFHLTH[sampt$DIET==1])
sampt$SLFWGHT <- sampt$SLFWGHT - mean(sampt$SLFWGHT[sampt$DIET==1])
sampt$WORKHARD <- sampt$WORKHARD - mean(sampt$WORKHARD[sampt$DIET==1])
sampt$GOODQUAL <- sampt$GOODQUAL - mean(sampt$GOODQUAL[sampt$DIET==1])
sampt$PHYSFIT <- sampt$PHYSFIT - mean(sampt$PHYSFIT[sampt$DIET==1])
sampt$PROUD <- sampt$PROUD - mean(sampt$PROUD[sampt$DIET==1])
sampt$LIKESLF <- sampt$LIKESLF - mean(sampt$LIKESLF[sampt$DIET==1])
sampt$ACCEPTED <- sampt$ACCEPTED - mean(sampt$ACCEPTED[sampt$DIET==1])
sampt$FEELLOVD <- sampt$FEELLOVD - mean(sampt$FEELLOV[sampt$DIET==1])

# fit the model to estimate ACE.1
fit <- glm( DISTR.2 ~  DIET + 
   DISTR.1 + BLACK + NBHISP + GRADE + 
   SLFHLTH + SLFWGHT + WORKHARD + GOODQUAL + PHYSFIT + PROUD + 
   LIKESLF + ACCEPTED + FEELLOVD +
   DISTR.1:DIET + BLACK:DIET + NBHISP:DIET + GRADE:DIET + 
   SLFHLTH:DIET + SLFWGHT:DIET + WORKHARD:DIET + GOODQUAL:DIET +
   PHYSFIT:DIET + PROUD:DIET + 
   LIKESLF:DIET + ACCEPTED:DIET + FEELLOVD:DIET,
   data=sampt)

# > print(summary(fit))
# 
# Call:
# glm(formula = DISTR.2 ~ DIET + DISTR.1 + BLACK + NBHISP + GRADE + 
#     SLFHLTH + SLFWGHT + WORKHARD + GOODQUAL + PHYSFIT + PROUD + 
#     LIKESLF + ACCEPTED + FEELLOVD + DISTR.1:DIET + BLACK:DIET + 
#     NBHISP:DIET + GRADE:DIET + SLFHLTH:DIET + SLFWGHT:DIET + 
#     WORKHARD:DIET + GOODQUAL:DIET + PHYSFIT:DIET + PROUD:DIET + 
#     LIKESLF:DIET + ACCEPTED:DIET + FEELLOVD:DIET, data = sampt)
# 
# Deviance Residuals: 
#      Min        1Q    Median        3Q       Max  
# -1.72288  -0.22801  -0.05723   0.18175   1.82660  
# 
# Coefficients:
#                 Estimate Std. Error t value Pr(>|t|)    
# (Intercept)    0.7184844  0.0075152  95.604  < 2e-16 ***
# DIET          -0.0153450  0.0131194  -1.170 0.242192    
# DISTR.1        0.5210525  0.0143816  36.230  < 2e-16 ***
# BLACK          0.0742012  0.0131869   5.627 1.92e-08 ***
# NBHISP         0.0241334  0.0159814   1.510 0.131071    
# GRADE          0.0021585  0.0039652   0.544 0.586210    
# SLFHLTH        0.0257745  0.0067034   3.845 0.000122 ***
# SLFWGHT        0.0026900  0.0077575   0.347 0.728781    
# WORKHARD      -0.0194426  0.0062698  -3.101 0.001938 ** 
# GOODQUAL       0.0154691  0.0110400   1.401 0.161212    
# PHYSFIT        0.0054135  0.0077812   0.696 0.486635    
# PROUD          0.0329261  0.0113575   2.899 0.003756 ** 
# LIKESLF        0.0225898  0.0072863   3.100 0.001942 ** 
# ACCEPTED       0.0159821  0.0077195   2.070 0.038462 *  
# FEELLOVD       0.0412613  0.0095958   4.300 1.74e-05 ***
# DIET:DISTR.1  -0.0106838  0.0296963  -0.360 0.719033    
# DIET:BLACK    -0.0003361  0.0327604  -0.010 0.991814    
# DIET:NBHISP    0.0238344  0.0346707   0.687 0.491826    
# DIET:GRADE     0.0023927  0.0091386   0.262 0.793465    
# DIET:SLFHLTH  -0.0297672  0.0150042  -1.984 0.047310 *  
# DIET:SLFWGHT  -0.0011506  0.0181310  -0.063 0.949404    
# DIET:WORKHARD  0.0392602  0.0147076   2.669 0.007620 ** 
# DIET:GOODQUAL  0.0260125  0.0246058   1.057 0.290477    
# DIET:PHYSFIT  -0.0266321  0.0171818  -1.550 0.121190    
# DIET:PROUD     0.0240841  0.0255893   0.941 0.346650    
# DIET:LIKESLF   0.0086665  0.0153017   0.566 0.571161    
# DIET:ACCEPTED -0.0015125  0.0164325  -0.092 0.926664    
# DIET:FEELLOVD -0.0125126  0.0209234  -0.598 0.549850    
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
# 
# (Dispersion parameter for gaussian family taken to be 0.1410809)
# 
#     Null deviance: 1303.09  on 5999  degrees of freedom
# Residual deviance:  842.54  on 5972  degrees of freedom
# AIC: 5306.7
# 
# Number of Fisher Scoring iterations: 2