No one replied to this, so I'll try again, with a simple example. I calculate a set of log odds ratios, and turn them into a data frame as follows:

> library(vcdExtra)
> (lor.CM <- loddsratio(CoalMiners))
log odds ratios for Wheeze and Breathlessness by Age

   25-29    30-34    35-39    40-44    45-49    50-54    55-59    60-64
3.695261 3.398339 3.140658 3.014687 2.782049 2.926395 2.440571 2.637954
>
> (lor.CM.df <- as.data.frame(lor.CM))
  Wheeze Breathlessness   Age      LOR        ASE
1  W:NoW          B:NoB 25-29 3.695261 0.16471778
2  W:NoW          B:NoB 30-34 3.398339 0.07733658
3  W:NoW          B:NoB 35-39 3.140658 0.03341311
4  W:NoW          B:NoB 40-44 3.014687 0.02866111
5  W:NoW          B:NoB 45-49 2.782049 0.01875164
6  W:NoW          B:NoB 50-54 2.926395 0.01585918
7  W:NoW          B:NoB 55-59 2.440571 0.01452057
8  W:NoW          B:NoB 60-64 2.637954 0.02159903

Now I want to fit a linear model by WLS, LOR ~ Age, which can do directly as

> lm(LOR ~ as.numeric(Age), weights=1/ASE, data=lor.CM.df)

Call:
lm(formula = LOR ~ as.numeric(Age), data = lor.CM.df, weights = 1/ASE)

Coefficients:
    (Intercept)  as.numeric(Age)
         3.5850          -0.1376

But, I want to do the fitting in my own function, the simplest version is

my.lm <- function(formula, data, subset, weights) {
        lm(formula, data, subset, weights)
}

But there is obviously some magic about formula objects and evaluation environments, because I don't understand why this doesn't work.

> my.lm(LOR ~ as.numeric(Age), weights=1/ASE, data=lor.CM.df)
Error in model.frame.default(formula = formula, data = data, subset = subset, :
  invalid type (closure) for variable '(weights)'
>

A second question: Age is a factor, and as.numeric(Age) gives me 1:8.

What simple expression on lor.CM.df$Age would give me either the lower
limits (here: seq(25, 60, by = 5)) or midpoints of these Age intervals
(here: seq(27, 62, by = 5))?

best,
-Michael


On 12/16/2010 9:43 AM, Michael Friendly wrote:
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.

Reply via email to