0

I am trying to write a custom function. I have the working code:

m <- lm(age ~ religion * degree * country + gender, WVS)
e_df <- data.frame(e <- effect("religion:degree:country", m))
e_df$predictor <- paste(e_df[,3], e_df[,2], e_df[,1], sep="\n\n") # merge first three columns
ifelse(names(e_df)=="fit", "age", names(e_df))->names(e_df) # rename column 1
print(subset(e_df, select = c(predictor, age, se)), digits = 3)

I tried to write it into the custom function: 

threewaycovar<- function(outcome, V1, V2, V3, covar, d) 
{require(effects)
m <- lm(outcome ~ V1 * V2 * V3 + covar, d)
e <- effect("V1:V2:V3", m)
e_df <- data.frame(e)
e_df$predictor <- paste(e_df[,3], e_df[,2], e_df[,1], sep="\n\n") # merge     first three columns
ifelse(names(e_df)=="fit", "outcome", names(e_df))->names(e_df) # rename column 1
op <-(subset(e_df, select = c(predictor, outcome, se)))
return(op)}

However, the custom function keeps giving me the error message

" Error in eval(expr, envir, enclos) : object 'outcome' not found".

I am not sure why I am getting this error message. I tried to run it using 

threewaycovar(WVS$age, WVS$religion, WVS$degree, WVS$country, WVS$gender, WVS)

I saw this post Error in eval(expr, envir, enclos) - contradiction? and took out a reference to column numbers that relates to outcome but my function still is not working. Can anyone tell me why my function is not working?

Edit in response to comments: WVS is a default package included in R so my data is reproducible. My traceback is 

14 eval(expr, envir, enclos) 
13 eval(predvars, data, env) 
12 model.frame.default(outcome ~ V1 * V2 * V3 + covar + (V1 + V2 + 
V3 + covar), data = structure(list(poverty = structure(c(1L, 2L, 1L, 3L, 1L, 2L, 3L, 1L, 1L, 1L, 1L, 2L, 1L, 3L, 1L, 3L, 3L, 1L, 2L, 2L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 2L, 1L, 3L, 1L,  ... 
11 model.frame(outcome ~ V1 * V2 * V3 + covar + (V1 + V2 + V3 + 
covar), data = structure(list(poverty = structure(c(1L, 2L, 1L, 3L, 1L, 2L, 3L, 1L, 1L, 1L, 1L, 2L, 1L, 3L, 1L, 3L, 3L, 1L, 2L, 2L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 2L, 1L, 3L, 1L, 1L,  ... 
10 eval(expr, envir, enclos) 
9 eval(call("model.frame", ff, data = data, subset = subset, na.action = naa), 
envir) 
8 expand.model.frame(mod, all.predictors) 
7 na.omit(expand.model.frame(mod, all.predictors)) 
6 Analyze.model(focal.predictors, mod, xlevels, default.levels, 
formula.rhs, partial.residuals = partial.residuals, quantiles = quantiles, 
x.var = x.var, data = data, typical = typical) 
5 Effect.lm(predictors, mod, vcov. = vcov., ...) 
4 Effect(predictors, mod, vcov. = vcov., ...) 
3 effect.default("V1:V2:V3", m) 
2 effect("V1:V2:V3", m) 
1 threewaycovar(WVS$age, WVS$religion, WVS$degree, WVS$country, 
WVS$gender, WVS)

I am trying to make a function that will give me the estimated means and standard errors for a three way interaction controlling for at least one co-variate. Above is what I have tried thus far.

M--
  • 25,431
  • 8
  • 61
  • 93
Jim
  • 1
  • 1
  • 4
  • 1
    You should provide some sample input data to make your problem [reproducible](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) and make it easier to help you. The output of `traceback()` would be useful to see which function is throwing the error. It would also be better to describe what you are trying to do rather than just showing how you attempted it. There may be better ways, – MrFlick May 17 '16 at 23:13
  • 1
    `ifelse(names(e_df)=="fit", "outcome", names(e_df))->names(e_df)` – HubertL May 17 '16 at 23:20
  • @HubertL thank you for your comment but I am still getting the same error message. – Jim May 17 '16 at 23:34
  • @MrFlick I used one of the defualt packages to make my problem reproducible and I added the traceback among other edits above. Do you have a better way to get the estimated means and standard errors for a three way interaction controlling for covariates? summary SE does not work in this case. – Jim May 17 '16 at 23:34
  • Where does `WVS` come from? It's not part of the base or stats packages. – MrFlick May 17 '16 at 23:38
  • 1
    @Jim `WVS` is a data set, not a package. It's included in the `effects` package, not a default package. – Gregor Thomas May 17 '16 at 23:42
  • 1
    As to your problem, most functions in R where a user can specify a model follow one of two possibilities: (1) the model is specified by a `formula` and a data set is passed in, e.g., `lm`, `glm`, `lmer`, `polr`, ... (2) the response is passed in as a single vector and a design matrix is passed in for the predictors, e.g., `xgboost`, `glmnet`, ... Your function is trying to go a third way. Don't. Just do the same thing that has been done and works for everyone else. Pick way (1), it will be easier on the user than way (2) with all those interactions. – Gregor Thomas May 17 '16 at 23:50

1 Answers1

4

The main problem you are having is you can't just throw variables into formulas or expressions. Here's a re-write of the function where you pass columns names as strings rather than vector values. 

library(effects)
threewaycovar <- function(outcome, effects, covar, d) {
    stopifnot(all(c(outcome, effects, covar) %in% names(d)))
    form <- as.formula(paste(outcome,"~", paste(effects,collapse="*"), "+", paste(covar, collapse="+")))
    m <- lm(form, d)
    e <- effect(paste(effects, collapse=":"), m)
    e_df <- data.frame(e)
    e_df$predictor <- apply(e_df[,rev(seq_along(effects))], 1, paste, collapse="\n\n") # merge     first three columns
    ifelse(names(e_df)=="fit", "outcome", names(e_df))->names(e_df) # rename column 1
    op <-e_df[,c("predictor", "outcome", "se")]
    return(op)
}

threewaycovar("age", c("religion","degree","country"), "gender", WVS)

Here we build our formula dynamically as a string, preserving the original names from the data set.

MrFlick
  • 195,160
  • 17
  • 277
  • 295