help with ols and contrast functions in Design library

Dear All,

I'm trying to use the ols function in the Design library (version  
2.1.1) of R to estimate parameters of a linear model, and then use the  
contrast function in the same library to test various contrasts.

As a simple example, suppose I have three factors:  feature (3  
levels), group (2 levels), and patient (3 levels).  Patient is coded  
as a non-unique identifier and is therefore nested within group.

response <- rnorm(length(example$LOG_ABUNDANCE), mean = 12)
feature <- rep(c(1,2,3), 6)
group <- c(rep(c(1,2),each=9))
patient <- rep(rep(c(1,2,3), each=3),2)

myData <- data.frame(patient=factor(patient), group=factor(group),
feature=factor(feature), response=response)

I use the ols command to fit the linear model, but I receive the  
following error.

fit <- ols(response ~ feature*group + group/patient, myData)

 > fit <- ols(response ~ feature*group + group/patient, myData)

Error in if (!length(fname) || !any(fname == zname)) { :
   missing value where TRUE/FALSE needed

Because of this, I tried using a unique identifier for patient using  
the following command.

myData$group.patient <- with(myData, group:patient)[drop=TRUE]

Running the same model with this factor will correct the error, but  
leaves me with an 'NA' for one of the estimated model parameters.

 > fit2 <- ols(response ~ feature*group + group.patient, myData)
 > fit2

Linear Regression Model

ols(formula = response ~ feature * group + group.patient, data = myData)

          n Model L.R.       d.f.         R2      Sigma
         18      4.659         10     0.2281      1.122

Residuals:
     Min      1Q  Median      3Q     Max
-1.1466 -0.5854 -0.2545  0.6834  1.4900

Coefficients:
                       Value Std. Error          t  Pr(>|t|)
Intercept           12.7116  9.442e-01  1.346e+01 2.928e-06
feature=2           -0.4795  1.370e+00 -3.500e-01 7.367e-01
feature=3           -0.0948  1.389e+00 -6.828e-02 9.475e-01
group=2             -0.7218  3.586e+15 -2.013e-16 1.000e+00
group.patient=1:2   -1.1455  1.120e+00 -1.023e+00 3.405e-01
group.patient=1:3   -0.5619  9.894e-01 -5.679e-01 5.879e-01
group.patient=2:1   -0.1402  3.586e+15 -3.909e-17 1.000e+00
group.patient=2:2   -0.1699  3.586e+15 -4.738e-17 1.000e+00
group.patient=2:3        NA  1.438e+00         NA        NA
feature=2 * group=2  0.1224  1.669e+00  7.330e-02 9.436e-01
feature=3 * group=2 -0.1970  3.586e+15 -5.494e-17 1.000e+00


When I try to test a contrast based on this fit, the 'NA' apparently  
prevents the estimation of the contrast.

 > contrast(fit2, list(group='1', feature=levels(myData$feature),  

group.patient=levels(myData$group.patient)), list(group='2',  
feature=levels(myData$feature), group.patient=levels(myData 
$group.patient)), type="average")
   Contrast         S.E. Lower Upper  t Pr(>|t|)
1       NA 2.390489e+15    NA    NA NA       NA

Error d.f.= 8

Any suggestions?

Sincerely,
Tim Clough



	[[alternative HTML version deleted]]