[R] help with ols and contrast functions in Design library

[Frank E Harrell Jr]

Timothy Clough wrote:
> 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

Sorry, Design, and its replacement rms, do not support nested effects. 
Also, any model that results in an NA as a parameter estimate will not 
work properly in Design/rms.

Frank

> 
> 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]]
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 


-- 
Frank E Harrell Jr   Professor and Chair           School of Medicine
                      Department of Biostatistics   Vanderbilt University