假设我有一个响应变量和一个包含三个协变量的数据(作为玩具示例):
y = c(1,4,6)
d = data.frame(x1 = c(4,-1,3), x2 = c(3,9,8), x3 = c(4,-4,-2))
我想对数据进行线性回归:
fit = lm(y ~ d$x1 + d$x2 + d$y2)
有没有办法编写公式,这样我就不必写出每个协变量？例如,像
fit = lm(y ~ d)
(我希望数据框中的每个变量都是协变量.)我问,因为我的数据框中实际上有50个变量,所以我想避免写出来x1 + x2 + x3 + etc.
在R predict.lm中,根据线性回归的结果计算预测,并提供计算这些预测的置信区间.根据手册,这些间隔是基于拟合的误差方差,而不是基于系数的误差间隔.
另一方面,基于逻辑和泊松回归计算预测的predict.glm(在其他几个中)没有置信区间的选项.我甚至很难想象如何计算这种置信区间,以便为泊松和逻辑回归提供有意义的见解.
是否存在为此类预测提供置信区间有意义的情况？他们怎么解释？这些案例中的假设是什么？
我正在使用该twang包来创建倾向得分,这些得分在二项式glm中用作权重survey::svyglm.代码看起来像这样:
pscore <- ps(ppci ~ var1+var2+.........., data=dt....)
dt$w <- get.weights(pscore, stop.method="es.mean")
design.ps <- svydesign(ids=~1, weights=~w, data=dt,)
glm1 <- svyglm(m30 ~ ppci, design=design.ps,family=binomial)
这会产生以下警告:
Warning message:
In eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
有谁知道我做错了什么？
我不确定这个消息在stats.SE上是否会更好,但总的来说我想我会先在这里尝试一下.
我想估计一个响应的协变量效应,其值取[0,1]中的值.也就是说,响应变量的值存在于0-1(包括)之间.我想使用Papke和Wooldridge(1996)描述的分数logit模型,见下文:
http://faculty.smu.edu/millimet/classes/eco6375/papers/papke%20wooldridge%201996.pdf
是否有R函数(或库)来促进分数logit模型的估计？我能glm()以某种方式修改吗？
我很欣赏@Jibler的评论 - 这可以从分数logit模型得到估计的beta值.但是,正如@Ben指出的那样,鉴于此规范,将无法正确估计SE.
我认为这是一个更受欢迎的经济学模型,因此STATA期刊撰稿人对此进行了充分讨论: http://fmwww.bc.edu/EC-C/S2013/823/EC823.S2013.nn06.slides.pdf http:/ /www.stata.com/meeting/germany10/germany10_buis.pdf
我能够从Papke和Wooldridge 401k计划示例中获得数据(见下文).在我看来,至少在分数logit模型中的稳健性是通过方差的三明治估计 - Papke和Wooldridge的方程(9)得到的.也就是说,等式(10)继续展示如何通过将估计的vcov矩阵与标准glm(...,family=binomial(link=logit))拟合预先乘以Pearson残差的估计来获得稳健性.
Buis的幻灯片似乎sandwich()使用参数vce(robust)实现了小数logit估计器的一种形式.这些与sandwich()R 中的函数的应用完全一致,与标准二项式GLM一致.我假设,但我不确定,因为我不是STATA,这与Baum的论点一样简单robust吗？如果有人拥有STATA并且可以检查这将有所帮助.family=quasibinomialGLM给出的模型给出了略微不同的SE估计值.但它似乎也是分数logit模型的均值/方差参数的合理估计.
下面是一些R代码,它复制了上面Buis文章中给出的数据拟合(它还显示了准二项式模型如何给出略微不同的SE估计):
##
## Replicate what some STATA Journal editors call "fractional logit"
## get data from: "http://fmwww.bc.edu/repec/bocode/k/k401.dta"
##
library(sandwich)
library(foreign)
X <- read.dta("F:/ProportionsDepVar/k401.dta")
class(X)
names(X)
dim(X)
X$totemp1 <- X$totemp/10000
glmfit <- glm(prate ~ mrate + totemp1 + age + sole, family=binomial(link=logit), data=X)
summary(glmfit)
##
## And the SE's …我正在运行很多回归,我只对一个特定变量的系数和p值的影响感兴趣.因此,在我的脚本中,我希望能够从glm摘要中提取p值(获得系数本身很容易).我知道查看p值的唯一方法是使用summary(myReg).还有其他方法吗？
例如:
fit <- glm(y ~ x1 + x2, myData)
x1Coeff <- fit$coefficients[2] # only returns coefficient, of course
x1pValue <- ???
我曾尝试将其fit$coefficients作为矩阵处理,但我仍然无法简单地提取p值.
是否有可能做到这一点？
谢谢!
我正在尝试使用10倍交叉验证来估计逻辑回归.
#import libraries
library(car); library(caret); library(e1071); library(verification)
#data import and preparation
data(Chile)
chile <- na.omit(Chile) #remove "na's"
chile <- chile[chile$vote == "Y" | chile$vote == "N" , ] #only "Y" and "N" required
chile$vote <- factor(chile$vote) #required to remove unwanted levels
chile$income <- factor(chile$income) # treat income as a factor
目标是估计一个glm模型,该模型预测投票"Y"或"N"的结果取决于相关的解释变量,并且基于最终模型,计算混淆矩阵和ROC曲线以掌握不同阈值水平的模型行为.
型号选择导致:
res.chileIII <- glm(vote ~
sex +
education +
statusquo ,
family = binomial(),
data = chile)
#prediction
chile.pred <- predict.glm(res.chileIII, type = "response")
产生:
> head(chile.pred)
1 2 …以下是我正在使用的所有变量:
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 …我想强制特定的变量进入glm回归而不完全指定每个变量.我的真实数据集有大约200个变量.到目前为止,我还没有在我的在线搜索中找到这样的样本.
例如(只有3个变量):
n=200
set.seed(39)
samp = data.frame(W1 = runif(n, min = 0, max = 1), W2=runif(n, min = 0, max = 5))
samp = transform(samp, # add A
A = rbinom(n, 1, 1/(1+exp(-(W1^2-4*W1+1)))))
samp = transform(samp, # add Y
Y = rbinom(n, 1,1/(1+exp(-(A-sin(W1^2)+sin(W2^2)*A+10*log(W1)*A+15*log(W2)-1+rnorm(1,mean=0,sd=.25))))))
如果我想包含所有主要术语,这有一个简单的快捷方式:
glm(Y~., family=binomial, data=samp)
但是说我想要包括所有主要术语(W1,W2和A)加上W2 ^ 2:
glm(Y~A+W1+W2+I(W2^2), family=binomial, data=samp)
这有什么捷径吗？
[在发布之前编辑自己:]这个有效! glm(formula = Y ~ . + I(W2^2), family = binomial, data = samp)
好的,那么这个呢!
我想省略一个主要的术语变量,只包括两个主要术语(A,W2)和W2 ^ 2和W2 ^ 2:A:
glm(Y~A+W2+A*I(W2^2), family=binomial, data=samp)
显然只有几个变量不需要快捷方式,但我使用高维数据.当前数据集"仅"有200个变量,但其他一些变量有成千上万个.
当我训练只是使用时glm,一切正常,我甚至没有接近耗尽记忆.但是当我跑步时train(..., method='glm'),我的内存耗尽.
这是因为train为交叉验证的每次迭代存储了大量数据(或者无论trControl程序是什么)？我在看,trainControl我找不到如何防止这种...任何提示？我只关心性能总结和预测的响应.
(我知道它与存储参数调整网格搜索的每次迭代中的数据无关,因为我认为没有glm的网格.)
我有50个变量.这就是我在glm中使用它们的方式.
var = glm(Stuff ~ ., data=mydata, family=binomial)
但我想排除其中的两个.那么如何排除2具体？我希望会有这样的事情:
var = glm(Stuff ~ . # notthisstuff, data=mydata, family=binomial)
想法？