Timothy Clough
2009-Nov-05 19:54 UTC
[R] 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]]
Frank E Harrell Jr
2009-Nov-05 19:59 UTC
[R] help with ols and contrast functions in Design library
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 neededSorry, 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