# Load library
library(nlme)

# Response variable
response <- c(rnorm(5, 1, 3), rnorm(5, 6, 1), rnorm(5, 10, 5))

# Dataframe with two categorical variables
foo <- data.frame(response = response, X = rep(letters[1:3], each = 5), Y = rep(LETTERS[1:3], 
    each = 5))

# gls model allowing for different variance per level of Y
m1 <- gls(response ~ X, weights = varIdent(form = ~1 | Y), data = foo)

# anova command indicating the factor 'X' is significant
anova(m1)

## Denom. DF: 12 
##             numDF F-value p-value
## (Intercept)     1   478.2  <.0001
## X               2    10.9   0.002

# summary command indicating levels b and c are greater than level a of
# factor 'X'.  WHAT TEST DO I CITE IN A MANUSCRIPT FOR THESE CONTRASTS?
summary(m1)

## Generalized least squares fit by REML
##   Model: response ~ X 
##   Data: foo 
##     AIC   BIC logLik
##   70.55 73.46 -29.27
## 
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | Y 
##  Parameter estimates:
##      A      B      C 
## 1.0000 0.1762 1.7840 
## 
## Coefficients:
##              Value Std.Error t-value p-value
## (Intercept) -0.550     1.493  -0.368  0.7190
## Xb           6.324     1.516   4.173  0.0013
## Xc          11.707     3.053   3.835  0.0024
## 
##  Correlation: 
##    (Intr) Xb    
## Xb -0.985       
## Xc -0.489  0.482
## 
## Standardized residuals:
##     Min      Q1     Med      Q3     Max 
## -1.2975 -0.7330 -0.1790  0.8634  1.6068 
## 
## Residual standard error: 3.338 
## Degrees of freedom: 15 total; 12 residual