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@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
______________________________________________
R-help@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.
