天道酬勤,学无止境

R - how to pass formula to a with(df, glm(y ~ x)) construction inside a function

I'm using the mice package in R to multiply-impute some missing data. I need to be able to specify a formula that is passed to a with(df, glm(y ~ x)) construction inside of a function. This with() construction is the format used by the mice package to fit the regression model separately within each of the imputed datasets.

However, I cannot figure out the scoping problems preventing me from successfully passing the formula as an argument. Here is a reproducible example:

library(mice)

data(mtcars)
mtcars[5, 5] <- NA # introduce a missing value to be imputed

mtcars.imp = mice(mtcars, m = 5)

# works correctly outside of function
with(mtcars.imp, glm(mpg ~ cyl))

fit_model_mi = function(formula) {
  with(mtcars.imp, glm(formula))
}

# doesn't work when trying to pass formula into function   
fit_model_mi("mpg ~ cyl")

Also see here for the same question being asked on R help, although it does not receive an answer.

标签

评论

Try wrapping the formula in as.formula

fit_model_mi = function(formula) {
    with(mtcars.imp, glm(as.formula(formula)) )
}

Seems to work:

> fit_model_mi("mpg ~ cyl")
call :
with.mids(data = mtcars.imp, expr = glm(as.formula(formula)))

call1 :
mice(data = mtcars, m = 5)

nmis :
 mpg  cyl disp   hp drat   wt qsec   vs   am gear carb 
   0    0    0    0    1    0    0    0    0    0    0 

analyses :
[[1]]

Call:  glm(formula = as.formula(formula))

Coefficients:
(Intercept)          cyl  
     37.885       -2.876  

Degrees of Freedom: 31 Total (i.e. Null);  30 Residual
Null Deviance:      1126 
Residual Deviance: 308.3    AIC: 169.3

You can also attach your data by

attach(mtcars)

Result shown

fit_model_mi("mpg ~ cyl")
call :
with.mids(data = mtcars.imp, expr = glm(formula))

call1 :
mice(data = mtcars, m = 5)

nmis :
 mpg  cyl disp   hp drat   wt qsec   vs   am gear carb 
   0    0    0    0    1    0    0    0    0    0    0 

analyses :
[[1]]

Call:  glm(formula = formula)

Coefficients:
(Intercept)          cyl  
     37.885       -2.876  

Degrees of Freedom: 31 Total (i.e. Null);  30 Residual
Null Deviance:      1126 
Residual Deviance: 308.3    AIC: 169.3

受限制的 HTML

  • 允许的HTML标签:<a href hreflang> <em> <strong> <cite> <blockquote cite> <code> <ul type> <ol start type> <li> <dl> <dt> <dd> <h2 id> <h3 id> <h4 id> <h5 id> <h6 id>
  • 自动断行和分段。
  • 网页和电子邮件地址自动转换为链接。

相关推荐
  • 如何调试“对比只能应用于具有 2 个或更多级别的因素”错误?(How to debug "contrasts can be applied only to factors with 2 or more levels" error?)
    问题 以下是我正在使用的所有变量: str(ad.train) $ Date : Factor w/ 427 levels "2012-03-24","2012-03-29",..: 4 7 12 14 19 21 24 29 31 34 ... $ Team : Factor w/ 18 levels "Adelaide","Brisbane Lions",..: 1 1 1 1 1 1 1 1 1 1 ... $ Season : int 2012 2012 2012 2012 2012 2012 2012 2012 2012 2012 ... $ Round : Factor w/ 28 levels "EF","GF","PF",..: 5 16 21 22 23 24 25 26 27 6 ... $ Score : int 137 82 84 96 110 99 122 124 49 111 ... $ Margin : int 69 18 -56 46 19 5 50 69 -26 29 ... $ WinLoss : Factor w/ 2 levels "0","1": 2 2 1 2 2 2 2 2 1 2 ... $ Opposition : Factor w/ 18 levels "Adelaide","Brisbane Lions",..: 8 18 10
  • 如何调试“对比只能应用于具有 2 个或更多级别的因素”错误?(How to debug "contrasts can be applied only to factors with 2 or more levels" error?)
    问题 以下是我正在使用的所有变量: str(ad.train) $ Date : Factor w/ 427 levels "2012-03-24","2012-03-29",..: 4 7 12 14 19 21 24 29 31 34 ... $ Team : Factor w/ 18 levels "Adelaide","Brisbane Lions",..: 1 1 1 1 1 1 1 1 1 1 ... $ Season : int 2012 2012 2012 2012 2012 2012 2012 2012 2012 2012 ... $ Round : Factor w/ 28 levels "EF","GF","PF",..: 5 16 21 22 23 24 25 26 27 6 ... $ Score : int 137 82 84 96 110 99 122 124 49 111 ... $ Margin : int 69 18 -56 46 19 5 50 69 -26 29 ... $ WinLoss : Factor w/ 2 levels "0","1": 2 2 1 2 2 2 2 2 1 2 ... $ Opposition : Factor w/ 18 levels "Adelaide","Brisbane Lions",..: 8 18 10
  • 如何调试“对比只能应用于具有 2 个或更多级别的因素”错误?(How to debug "contrasts can be applied only to factors with 2 or more levels" error?)
    问题 以下是我正在使用的所有变量: str(ad.train) $ Date : Factor w/ 427 levels "2012-03-24","2012-03-29",..: 4 7 12 14 19 21 24 29 31 34 ... $ Team : Factor w/ 18 levels "Adelaide","Brisbane Lions",..: 1 1 1 1 1 1 1 1 1 1 ... $ Season : int 2012 2012 2012 2012 2012 2012 2012 2012 2012 2012 ... $ Round : Factor w/ 28 levels "EF","GF","PF",..: 5 16 21 22 23 24 25 26 27 6 ... $ Score : int 137 82 84 96 110 99 122 124 49 111 ... $ Margin : int 69 18 -56 46 19 5 50 69 -26 29 ... $ WinLoss : Factor w/ 2 levels "0","1": 2 2 1 2 2 2 2 2 1 2 ... $ Opposition : Factor w/ 18 levels "Adelaide","Brisbane Lions",..: 8 18 10
  • 如何调试“对比只能应用于具有 2 个或更多级别的因素”错误?(How to debug "contrasts can be applied only to factors with 2 or more levels" error?)
    问题 以下是我正在使用的所有变量: str(ad.train) $ Date : Factor w/ 427 levels "2012-03-24","2012-03-29",..: 4 7 12 14 19 21 24 29 31 34 ... $ Team : Factor w/ 18 levels "Adelaide","Brisbane Lions",..: 1 1 1 1 1 1 1 1 1 1 ... $ Season : int 2012 2012 2012 2012 2012 2012 2012 2012 2012 2012 ... $ Round : Factor w/ 28 levels "EF","GF","PF",..: 5 16 21 22 23 24 25 26 27 6 ... $ Score : int 137 82 84 96 110 99 122 124 49 111 ... $ Margin : int 69 18 -56 46 19 5 50 69 -26 29 ... $ WinLoss : Factor w/ 2 levels "0","1": 2 2 1 2 2 2 2 2 1 2 ... $ Opposition : Factor w/ 18 levels "Adelaide","Brisbane Lions",..: 8 18 10
  • R:将参数传递给R函数中的glm(R : Pass argument to glm inside an R function)
    问题 我想习惯于解决R中的问题。我想在函数内部调用函数glm()但它不起作用,显然是出于范围限制的原因,我没有设法通过函数assign()或eval()进行修复。 eval() 。 这是一个简化的版本: ao <- function (y, x, phi = seq (0,1,0.1), dataset, weights) { logLikvector <- rep(0,length(phi)) # vector of zeros to be replaced thereafter for (i in 1:length(phi)) { # loop to use glm() fit <- glm (y ~ x, data = dataset, family = binomial, weights = weights) logLikvector[i] <- logLik(fit) # get log likelihood } logLikvector } 现在我想在数据集中使用函数ao() ao (y = Prop, x = Age, dataset = mydata, weights = Total) 这不起作用,但是以下起作用: ao (y = mydata$Prop, x = mydata$Age, dataset = mydata, weights = mydata
  • 如何使 group_by 和 lm 快速?(How to make group_by and lm fast?)
    问题 这是一个样本。 df <- tibble( subject = rep(letters[1:7], c(5, 6, 7, 5, 2, 5, 2)), day = c(3:7, 2:7, 1:7, 3:7, 6:7, 3:7, 6:7), x1 = runif(32), x2 = rpois(32, 3), x3 = rnorm(32), x4 = rnorm(32, 1, 5)) df %>% group_by(subject) %>% summarise( coef_x1 = lm(x1 ~ day)$coefficients[2], coef_x2 = lm(x2 ~ day)$coefficients[2], coef_x3 = lm(x3 ~ day)$coefficients[2], coef_x4 = lm(x4 ~ day)$coefficients[2]) 这个数据很小,所以性能没有问题。 但是我的数据太大了,大约有 1,000,000 行和 200,000 个主题,而且这段代码很慢。 我认为原因不是lm的速度,而是很多主题和子集。 回答1 理论上 首先,您可以拟合具有多个 LHS 的线性模型。 其次,显式数据拆分并不是 group-by 回归的唯一方法(或推荐的方法)。 请参阅 R 回归分析:分析特定种族的数据和 R:为每个类别构建单独的模型。 因此
  • 如何在R中添加不同的趋势线?(How do I add different trend lines in R?)
    问题 我知道如何使用lm和abline函数添加线性趋势线,但是如何添加其他趋势线,例如对数,指数和幂趋势线? 回答1 这是我之前准备的: # set the margins tmpmar <- par("mar") tmpmar[3] <- 0.5 par(mar=tmpmar) # get underlying plot x <- 1:10 y <- jitter(x^2) plot(x, y, pch=20) # basic straight line of fit fit <- glm(y~x) co <- coef(fit) abline(fit, col="blue", lwd=2) # exponential f <- function(x,a,b) {a * exp(b * x)} fit <- nls(y ~ f(x,a,b), start = c(a=1, b=1)) co <- coef(fit) curve(f(x, a=co[1], b=co[2]), add = TRUE, col="green", lwd=2) # logarithmic f <- function(x,a,b) {a * log(x) + b} fit <- nls(y ~ f(x,a,b), start = c(a=1, b=1)) co <- coef(fit) curve(f(x
  • R错误,提示“模型未全部适合相同大小的数据集”(R error which says “Models were not all fitted to the same size of dataset”)
    问题 我创建了两个广义线性模型,如下所示: glm1 <-glm(Y ~ X1 + X2 + X3, family=binomial(link=logit)) glm2 <-glm(Y ~ X1 + X2, family=binomial(link=logit)) 然后,我使用anova函数: anova(glm2,glm1) 但收到错误消息: “ anova.glmlist(c(list(object),dotargs)中的错误,色散=色散,: 模型并非都适合于相同大小的数据集” 这是什么意思,我该如何解决? 我已经在代码的开头attach了数据集,因此两个模型都在同一个数据集上工作。 回答1 该错误的主要原因是当一个或多个预测变量中缺少值时。 在最新版本的R中,默认操作是忽略所有缺少任何值的行(以前的默认设置是产生错误)。 因此,例如,如果数据框有100行,并且X3中缺少一个值,那么您的模型glm1将适合99行数据(将X3丢失的行删除),但是glm2对象将适合全部100行数据(由于它不使用X3,因此不需要删除任何行)。 因此,由于两个模型适合不同的数据集(以及如何计算自由度等),因此anova函数会给您带来错误。 一种解决方案是创建一个仅包含将在至少一个模型中使用的列的新数据框,并删除所有缺少值的所有行(使用na.omit或na.exclude函数可以轻松实现这一点)
  • 在 R 中使用 sample.split 不正确的数据拆分和逻辑回归问题(Incorrect splitting of data using sample.split in R and issue with logistic regression)
    问题 我有2个问题。 当我尝试将我的数据拆分为测试集和训练集时,使用如下所示的sample.split ,采样完成得相当不清楚。 我的意思是数据 d 的长度为 392,因此,4:1 除法应显示 0.8*392=313.6 即测试集中的 313 或 314 行,但显示的长度为 304。有什么我可以失踪? require(caTools) set.seed(101) samplev = sample.split(d[,], SplitRatio= 0.80) train = subset(d, samplev == TRUE) test = subset(d, samplev == FALSE) 我正在尝试将拆分数据如下用于 R 中的逻辑回归任务,如下所示- #Training m <- glm(mpg01~ . -name, data= train, family = binomial(link = 'logit')) out2 <- predict.glm(m, test, type = "response") class2 <- vector() for (i in 1:length(out2)) { if(out2[i] >= 0.5) { class2[i] <- 1 } else { class2[i] <- 0 } } r2 <- table(class2, test
  • 如何在一个循环中拟合多个交互模型?(How to fit multiple interaction models in a loop?)
    问题 假设我有 3 个响应变量 A、C 和 M,我想为所有可能的模型拟合一个模型,即拟合 Y ~ A、Y ~ C、Y ~ M、Y ~ A * C、Y ~ A * M、Y ~ C * M 等。有没有一种快速的方法来做到这一点,而无需每次都手动指定交互? 我不想写 M1 = glm(Y ~ A , data = subs, family = "poisson") M2 = glm(Y ~ C , data = subs, family = "poisson") M3 = glm(Y ~ M , data = subs, family = "poisson") M4 = glm(Y ~ A*C , data = subs, family = "poisson") ... 实际上,我有 3 个以上的变量,并且想要某种循环,这是否可能。 谢谢 回答1 这应该有效: glmulti::glmulti( Y = "Y", xr = c("A", "C", "M"), data = subs, filename = "my_results", family = "poisson" ) 它将创建一个文本文件my_results.txt其中包含有关每个模型的信息。 您也可以将其与其他包、 leaps 、 bestglm其他包作为单行包一起使用。 回答2 这是一种函数式编程方法。 您创建数据
  • update() 具有局部协变量的函数内的模型(update() a model inside a function with local covariate)
    问题 我需要从函数内部更新回归模型。 理想情况下,该函数应该适用于任何类型的模型( lm 、 glm 、 multinom 、 clm )。 更准确地说,我需要添加一个或多个在函数内部定义的协变量。 这是一个例子。 MyUpdate <- function(model){ randData <- data.frame(var1=rnorm(length(model$residuals))) model2 <- update(model, ".~.+randData$var1") return(model2) } 这是一个使用示例 data(iris) model1 <- lm(Sepal.Length~Species, data=iris) model2 <- MyUpdate(model1) eval(expr,envir,enclos)中的错误:找不到对象“randData” 这是 glm 的另一个例子 model1 <- glm(Sepal.Length>5~Species, data=iris, family=binomial) model2 <- MyUpdate(model1) 任何想法? 回答1 问题是var1在数据框和模型的环境中MyUpdate但不在MyUpdate的环境中。 1)为避免此问题,不仅使用修订后的公式更新模型,还使用包含var1的修订数据框更新模型:
  • R - 在 data.table 中使用 glm(R - using glm inside a data.table)
    问题 我正在尝试在 data.table 中做一些 glm 以生成按关键因素分割的建模结果。 我一直在成功地做到这一点: 高水平glm glm(modellingDF,formula=Outcome~IntCol + DecCol,family=binomial(link=logit)) 具有单列的作用域 glm 建模DF[,列表(结果, 拟合=glm(x,公式=结果~IntCol ,family=binomial(link=logit))$fitted ), by=variable] 具有两个整数列的作用域 glm 建模DF[,列表(结果, 拟合=glm(x,公式=结果~IntCol + IntCol2 ,family=binomial(link=logit))$fitted ), by=variable] 但是,当我尝试使用我的十进制列在范围内执行高级 glm 时,它会产生此错误 Error in model.frame.default(formula = Outcome ~ IntCol + DecCol, data = x, : variable lengths differ (found for 'DecCol') 我想可能是由于分区的长度可变,所以我用一个可重现的例子进行了测试: library("data.table") testing<-data.table
  • 如何在R中循环/重复线性回归(How to Loop/Repeat a Linear Regression in R)
    问题 我已经弄清楚了如何用4个变量在R中创建一个表,该表用于多个线性回归。 每次回归的因变量(Lung)取自22,000列的csv表的一列。 自变量之一(血液)取自相似表的相应列。 每列代表一个特定基因的水平,这就是为什么它们如此之多的原因。 还有两个附加变量(每个患者的年龄和性别)。 当我输入线性回归方程式时,我使用lm(Lung [,1]〜Blood [,1] + Age + Gender),它适用于一个基因。 我正在寻找一种输入此方程式的方法,并让R计算肺和血液的所有其余列,并希望将系数输出到表格中。 任何帮助,将不胜感激! 回答1 您想运行22,000个线性回归并提取系数吗? 从编码的角度来看,这很容易做到。 set.seed(1) # number of columns in the Lung and Blood data.frames. 22,000 for you? n <- 5 # dummy data obs <- 50 # observations Lung <- data.frame(matrix(rnorm(obs*n), ncol=n)) Blood <- data.frame(matrix(rnorm(obs*n), ncol=n)) Age <- sample(20:80, obs) Gender <- factor(rbinom(obs, 1, .5
  • ggplot2:如何绘制正交回归线?(ggplot2: How to plot an orthogonal regression line?)
    问题 我已经在两个不同的视觉感知测试中测试了大量参与者——现在,我想看看这两个测试的表现在多大程度上相关。 为了可视化相关性,我使用ggplot()在 R 中绘制了一个散点图,并拟合了一条回归线(使用stat_smooth() )。 但是,由于我的x和y变量都是性能度量,因此在拟合回归线时我需要将它们都考虑在内——因此,我不能使用简单的线性回归(使用stat_smooth(method="lm") ),但是而是需要拟合正交回归(或总最小二乘法)。 我该怎么做呢? 我知道我可以在stat_smooth()指定formula ,但我不知道要使用什么公式。 据我了解,没有任何预设方法( lm, glm, gam, loess, rlm )适用。 回答1 事实证明,您可以从 (x,y) 的主成分分析中提取斜率和截距,如下所示。 这只是稍微简单一点,在基础 R 中运行,并给出与在MethComp使用Deming(...)相同的结果。 # same `x and `y` as @user20650's answer df <- data.frame(y, x) pca <- prcomp(~x+y, df) slp <- with(pca, rotation[2,1] / rotation[1,1]) int <- with(pca, center[2] - slp*center[1])
  • R - model.frame() 和非标准评估(R - model.frame() and non-standard evaluation)
    问题 我对我试图编写的函数的行为感到困惑。 我的例子来自survival包,但我认为这个问题比这更普遍。 基本上,以下代码 library(survival) data(bladder) ## this will load "bladder", "bladder1" and "bladder2" mod_init <- coxph(Surv(start, stop, event) ~ rx + number, data = bladder2, method = "breslow") survfit(mod_init) 将产生一个我感兴趣的对象。但是,当我在函数中编写它时, my_function <- function(formula, data) { mod_init <- coxph(formula = formula, data = data, method = "breslow") survfit(mod_init) } my_function(Surv(start, stop, event) ~ rx + number, data = bladder2) 该函数将在最后一行返回错误: Error in eval(predvars, data, env) : invalid 'envir' argument of type 'closure' 10 eval
  • 如何从R中的glm模型中检索相关矩阵(How to retrieve correlation matrix from glm models in R)
    问题 我正在使用 nlme 包中的 gls 函数。 您可以复制并粘贴以下代码以重现我的分析。 library(nlme) # Needed for gls function # Read in wide format tlc = read.table("http://www.hsph.harvard.edu/fitzmaur/ala2e/tlc.dat",header=FALSE) names(tlc) = c("id","trt","y0","y1","y4","y6") tlc$trt = factor(tlc$trt, levels=c("P","A"), labels=c("Placebo","Succimer")) # Convert to long format tlc.long = reshape(tlc, idvar="id", varying=c("y0","y1","y4","y6"), v.names="y", timevar="time", direction="long") # Create week numerical variable tlc.long$week = tlc.long$time-1 tlc.long$week[tlc.long$week==2] = 4 tlc.long$week[tlc.long$week==3] = 6 tlc
  • R 或 Python - 循环测试数据 - 未来 24 小时的预测验证(每天 96 个值)(R or Python - loop the test data - Prediction validation next 24 hours (96 values each day))
    问题 我有一个大数据集,低于训练和测试数据集 train_data 是从 2016-01-29 到 2017-12-31 head(train_data) date Date_time Temp Ptot JFK AEH ART CS CP 1 2016-01-29 2016-01-29 00:00:00 30.3 1443.888 52.87707 49.36879 28.96548 6.239999 49.61212 2 2016-01-29 2016-01-29 00:15:00 30.3 1410.522 49.50248 49.58356 26.37977 5.024000 49.19649 3 2016-01-29 2016-01-29 00:30:00 30.3 1403.191 50.79809 49.04253 26.15317 5.055999 47.48126 4 2016-01-29 2016-01-29 00:45:00 30.3 1384.337 48.88359 49.14100 24.52135 5.088000 46.19261 5 2016-01-29 2016-01-29 01:00:00 30.1 1356.690 46.61842 48.80624 24.28208 5.024000 43.00352 6 2016-01-29 2016
  • 将 statsmodel 估计与 scikit-learn 交叉验证一起使用,这可能吗?(Using statsmodel estimations with scikit-learn cross validation, is it possible?)
    问题 我将此问题发布到 Cross Validated 论坛,后来意识到这可能会在 stackoverlfow 中找到合适的受众。 我正在寻找一种方法,我可以使用fit蟒蛇statsmodel ontained饲料为目标(结果) cross_val_score scikit学习cross_validation方法? 附加的链接表明它可能是可能的,但我没有成功。 我收到以下错误 estimator 应该是一个在 0x7fa6e801c590 处实现 'fit' 方法的估计器 statsmodels.discrete.discrete_model.BinaryResultsWrapper 对象被传递 参考这个链接 回答1 实际上,您不能直接在statsmodels对象上使用cross_val_score ,因为接口不同:在 statsmodels 中 训练数据直接传入构造函数一个单独的对象包含模型估计的结果 但是,您可以编写一个简单的包装器,使statsmodels对象看起来像sklearn估计器: import statsmodels.api as sm from sklearn.base import BaseEstimator, RegressorMixin class SMWrapper(BaseEstimator, RegressorMixin): """ A
  • 对来自小鼠 (R) 的部分估算数据运行 glm.mids(Run glm.mids on a subset of imputed data from mice (R))
    问题 当我尝试在mids插补对象的子集上运行glm.mids时出现错误: library(mice) imp2 = mice(nhanes) glm.mids( (hyp==2)~bmi+chl, data=imp2, subset=(age==1) ) 给出神秘的错误信息 "Error in eval(expr, envir, enclos) : ..1 used in an incorrect context, no ... to look in" 即使语法适用于原始数据集上的常规glm : glm( (hyp==2)~bmi+chl, data=nhanes, subset=(age==1) ) 文档?glm.mids没有专门解决subset但说您可以将其他参数传递给glm 。 如果我不能将subset与glm.mids一起glm.mids ,是否有一种直接对mids列表对象进行子集化的好方法? 回答1 我冒昧地重写了glm.mids 。 它有点笨拙。 这个问题似乎源于将属性传递给 glm 的隐含性质。 另请参阅这些帖子: https://stat.ethz.ch/pipermail/r-help/2003-November/041537.html http://r.789695.n4.nabble.com/Question-on-passing-the-subset
  • 为 geom_smooth gam 中的每一行获取调整后的 r 平方值(Getting adjusted r-squared value for each line in a geom_smooth gam)
    问题 我使用 ggplot2 生成了下图。 PlotEchi = ggplot(data=Echinoidea, aes(x=Year, y=mean, group = aspect, linetype = aspect, shape=aspect)) + geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=.025, position=pd) + geom_point(position=pd, size=2) + geom_smooth(method = "gam", formula = y~s(x, k=3), se=F, size = 0.5,colour="black") + xlab("") + ylab("Abundance (mean +/- SE)") + facet_wrap(~ species, scales = "free", ncol=1) + scale_y_continuous(limits=c(min(y=0), max(Echinoidea$mean+Echinoidea$se))) + scale_x_continuous(limits=c(min(Echinoidea$Year-0.125), max(Echinoidea$Year+0.125)))