> data_y group age y 1 1 31 17.05 2 1 23 4.96 3 1 27 10.40 4 1 28 11.05 5 1 22 0.26 6 1 24 2.51 7 2 23 -0.87 8 2 22 -10.74 9 2 22 -3.27 10 2 25 -1.97 11 2 27 7.50 12 2 20 -7.25 > > age_c <- age-mean(age) > num_grp <- max(group) > group <- factor(group) > > grp_rep <- rep(0,num_grp) > grp_mean_age <- rep(0,num_grp) > grp_se_age <- rep(0,num_grp) > grp_mean_y <- rep(0,num_grp) > grp_se_y <- rep(0,num_grp) > cum_rep <- 0 > E_XX <- 0 > > for (i1 in 1:num_grp) { + grp_rep[i1] <- length(y[group==i1]) + grp_mean_age[i1] <- mean(age[group==i1]) + grp_se_age[i1] <- sd(age[group==i1])/sqrt(grp_rep[i1]) + grp_mean_y[i1] <- mean(y[group==i1]) + grp_se_y[i1] <- sd(y[group==i1])/sqrt(grp_rep[i1]) + for (i2 in (cum_rep+1):(cum_rep+grp_rep[i1])) E_XX <- E_XX + (age[i2]-grp_mean_age[i1])^2 + + cum_rep <- cum_rep + grp_rep[i1] + } > > print (cbind(grp_mean_age,grp_se_age)) grp_mean_age grp_se_age [1,] 25.83333 1.400397 [2,] 23.16667 1.013794 > print (cbind(grp_mean_y,grp_se_y)) grp_mean_y grp_se_y [1,] 7.705000 2.554291 [2,] -2.766667 2.540106 > > drop1(ancv_y.mod0 <- lm(y ~ age_c + group, test="F")) Single term deletions Model: y ~ age_c + group Df Sum of Sq RSS AIC 70.39 27.230 age_c 1 318.91 389.30 45.753 group 1 71.79 142.18 33.666 Warning message: In lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) : extra arguments test are just disregarded. > > ancv_y.mod1 <- lm(y ~ -1 + age_c + group) > anova(ancv_y.mod1) Analysis of Variance Table Response: y Df Sum Sq Mean Sq F value Pr(>F) age_c 1 576.09 576.09 73.6594 1.257e-05 *** group 2 144.95 72.47 9.2666 0.006527 ** Residuals 9 70.39 7.82 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > summary.lm(ancv_y.mod1) Call: lm(formula = y ~ -1 + age_c + group) Residuals: Min 1Q Median 3Q Max -5.7731 -0.9902 0.1395 1.8254 3.0374 Coefficients: Estimate Std. Error t value Pr(>|t|) age_c 1.8859 0.2953 6.386 0.000127 *** group1 5.1905 1.2077 4.298 0.001997 ** group2 -0.2521 1.2077 -0.209 0.839270 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2.797 on 9 degrees of freedom Multiple R-squared: 0.9111, Adjusted R-squared: 0.8814 F-statistic: 30.73 on 3 and 9 DF, p-value: 4.648e-05 > > ancv_y.mod2 <- lm(y ~ -1 + group) > anova(ancv_y.mod2) Analysis of Variance Table Response: y Df Sum Sq Mean Sq F value Pr(>F) group 2 402.13 201.06 5.1648 0.0288 * Residuals 10 389.30 38.93 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > summary.lm(ancv_y.mod2) Call: lm(formula = y ~ -1 + group) Residuals: Min 1Q Median 3Q Max -7.9733 -4.6613 0.1467 2.8575 10.2667 Coefficients: Estimate Std. Error t value Pr(>|t|) group1 7.705 2.547 3.025 0.0128 * group2 -2.767 2.547 -1.086 0.3029 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 6.239 on 10 degrees of freedom Multiple R-squared: 0.5081, Adjusted R-squared: 0.4097 F-statistic: 5.165 on 2 and 10 DF, p-value: 0.0288 > > ancv_y.mod3 <- lm(y ~ age_c) > anova(ancv_y.mod3) Analysis of Variance Table Response: y Df Sum Sq Mean Sq F value Pr(>F) age_c 1 576.09 576.09 40.519 8.187e-05 *** Residuals 10 142.18 14.22 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > summary.lm(ancv_y.mod3) Call: lm(formula = y ~ age_c) Residuals: Min 1Q Median 3Q Max -7.5138 -0.3365 0.3053 1.4438 5.9081 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.4692 1.0885 2.268 0.0467 * age_c 2.2782 0.3579 6.365 8.19e-05 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 3.771 on 10 degrees of freedom Multiple R-squared: 0.8021, Adjusted R-squared: 0.7823 F-statistic: 40.52 on 1 and 10 DF, p-value: 8.187e-05 > > coef(ancv_y.mod1) age_c group1 group2 1.8858922 5.1904771 -0.2521437 > mse <- deviance(ancv_y.mod1)/df.residual(ancv_y.mod1) > mu_hat <- coef(ancv_y.mod1)[2:(num_grp+1)] > beta_hat <- coef(ancv_y.mod1)[1] > > > yb_adj <- mu_hat > se_yb_adj <- sqrt(mse*((1/grp_rep)+((grp_mean_age-mean(age))^2)/E_XX)) > > print(cbind(yb_adj,se_yb_adj)) yb_adj se_yb_adj group1 5.1904771 1.207708 group2 -0.2521437 1.207708 > > diff_adj_mean <- rep(0,(num_grp*(num_grp-1)/2)) > se_diff_adj_mean <- rep(0,(num_grp*(num_grp-1)/2)) > trt_id1 <- rep(0,(num_grp*(num_grp-1)/2)) > trt_id2 <- rep(0,(num_grp*(num_grp-1)/2)) > t_diff <- rep(0,(num_grp*(num_grp-1)/2)) > p_diff <- rep(0,(num_grp*(num_grp-1)/2)) > p_diff_bon <- rep(0,(num_grp*(num_grp-1)/2)) > pair_id <- 0 > > for (i1 in 1:(num_grp-1)) { + for (i2 in (i1+1):num_grp) { + pair_id <- pair_id+1 + + diff_adj_mean[pair_id] <- yb_adj[i1]-yb_adj[i2] + + se_diff_adj_mean[pair_id] <- + sqrt(mse*((1/grp_rep[i1] + 1/grp_rep[i2])+(((grp_mean_age[i1]-grp_mean_age[i2])^2)/E_XX))) + + t_diff[pair_id] <- diff_adj_mean[pair_id]/se_diff_adj_mean[pair_id] + p_diff[pair_id] <- 1-pt(abs(t_diff[pair_id]),df.residual(ancv_y.mod1)) + p_diff_bon[pair_id] <- min((((num_grp*(num_grp-1)/2))*p_diff[pair_id]),1) + + trt_id1[pair_id] <- i1 + trt_id2[pair_id] <- i2 + }} > > > > print(cbind(trt_id1,trt_id2,diff_adj_mean,se_diff_adj_mean, + t_diff,p_diff,p_diff_bon)) trt_id1 trt_id2 diff_adj_mean se_diff_adj_mean t_diff p_diff [1,] 1 2 5.442621 1.796453 3.029649 0.007127396 p_diff_bon [1,] 0.007127396 > > >