In the vcdExtra package on R-Forge, I have functions and generic methods
for calculating log odds ratios
for R x C x strata tables. I'd like to define methods for fitting
weighted lm()s to the resulting loddsratio objects,
but I'm having problems figuring out how to do this generally.
# install.packages("vcdExtra", repos="http://R-Forge.R-Project.org";)
library(vcdExtra)
> fung.lor <- loddsratio(Fungicide)
> fung.lor
log odds ratios for group and outcome by sex, strain
strain
sex 1 2
M -1.596015 -0.8266786
F -1.386294 -0.6317782
>
> fung.lor.df <- as.data.frame(fung.lor)
> fung.lor.df
group outcome sex strain LOR ASE
1 Control:Treated Tumor:NoTumor M 1 -1.5960149 0.7394909
2 Control:Treated Tumor:NoTumor F 1 -1.3862944 0.9574271
3 Control:Treated Tumor:NoTumor M 2 -0.8266786 0.6587325
4 Control:Treated Tumor:NoTumor F 2 -0.6317782 1.1905545
>
Now, I want to test whether the odds ratios differ by sex or strain, so
I do a weighted lm()
> fung.mod <- lm(LOR ~ sex + strain, data=fung.lor.df, weights=1/ASE^2)
> anova(fung.mod)
Analysis of Variance Table
Response: LOR
Df Sum Sq Mean Sq F value Pr(>F)
sex 1 0.00744 0.00744 112.3 0.05990 .
strain 1 0.84732 0.84732 12788.1 0.00563 **
Residuals 1 0.00007 0.00007
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
I tried to write a generic formula method to do this, but I keep running
into errors:
lor <- function(x, ...)
UseMethod("lor")
lor.formula <- function(formula, data, subset, weights,
model = TRUE, x = FALSE, y = FALSE,
contrasts = NULL, ...) {
data <- as.data.frame(data)
if (missing(weights)) {
if (! "ASE" %in% names(data)) stop("data does not contain an ASE column")
data$weights <- 1/data$ASE^2
}
lm(formula, data, subset, weights=weights,
model = model, x = x, y = y,
contrasts = contrasts, ...)
}
> lor(LOR ~ strain+sex, fung.lor)
Error in xj[i] : invalid subscript type 'closure'
> lor(LOR ~ strain+sex, fung.lor.df)
Error in xj[i] : invalid subscript type 'closure'
>
> traceback()
8: `[.data.frame`(list(LOR = c(-1.59601489210196, -1.38629436111989,
-0.826678573184468, -0.631778234183653), strain = c(1L, 1L, 2L,
2L), sex = c(1L, 2L, 1L, 2L), `(weights)` = c(1.82866556836903,
1.09090909090909, 2.30452674897119, 0.705507123112907)), function (x,
...)
UseMethod("subset"), , FALSE)
7: model.frame.default(formula = formula, data = data, subset = subset,
weights = weights, drop.unused.levels = TRUE)
6: model.frame(formula = formula, data = data, subset = subset,
weights = weights, drop.unused.levels = TRUE)
5: eval(expr, envir, enclos)
4: eval(mf, parent.frame())
3: lm(formula, data, subset, weights = weights, model = model, x = x,
y = y, contrasts = contrasts, ...)
2: lor.formula(LOR ~ strain + sex, fung.lor.df)
1: lor(LOR ~ strain + sex, fung.lor.df)
>
How can I make this work?
--
Michael Friendly Email: friendly AT yorku DOT ca
Professor, Psychology Dept.
York University Voice: 416 736-5115 x66249 Fax: 416 736-5814
4700 Keele Street Web: http://www.datavis.ca
Toronto, ONT M3J 1P3 CANADA
______________________________________________
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.
