我总共有95个潜在的预测变量,我想将这个数字减少到具有更强预测能力的变量.到目前为止,我的计划是编写一些代码:

• 在循环内选择6个随机预测变量并对它们执行逐步回归(方向=两者).
• 此循环将继续100,000次迭代,以确保输入每个可能的组合.
• 预测器(来自摘要命令)的重要性将基于p值.其中显着值<0.05被编码为'1'并且> 0.05对于6个预测变量(或更少)来说是'0'.预测变量名称保留在循环输出表中.

我似乎无法使用95列创建单个输出表,并使用为循环的每次迭代生成的6列来写入每个单独的行.

那么有没有办法添加到使用以下创建的数组:

results <- array(NA,c(100000,95)) 

列名由:

colnames(results)<-c(<inputdata>)

我正在使用auto.arima包中的函数拟合模型forecast.例如,我得到的是AR(1)的模型.然后我从这个模型中提取残差.这如何产生与原始矢量相同数量的残差？如果这是AR(1)模型,那么残差的数量应该比原始时间序列的维数小1.我错过了什么？

例:

require(forecast)
arprocess = as.numeric(arima.sim(model = list(ar=.5), n=100))
#auto.arima(arprocess, d=0, D=0, ic="bic", stationary=T)
#  Series: arprocess 
#  ARIMA(1,0,0) with zero mean     

#  Coefficients:
#          ar1
#       0.5198
# s.e.  0.0867

# sigma^2 estimated as 1.403:  log likelihood=-158.99
# AIC=321.97   AICc=322.1   BIC=327.18
r = resid(auto.arima(arprocess, d=0, D=0, ic="bic", stationary=T))
> length(r)
  [1] 100

更新:深入研究代码auto.arima,我看到它使用Arima了反过来使用的代码stats:::arima.因此,问题实际上是如何stats:::arima计算第一次观察的残差？

我想知道,如果SAS可以将所有数据集变量都包含在回归模型中而无需全部输入.我以前用过R,我想要一样的东西:

model <- lm(y ~ ., data = d)

但在SAS.提前致谢!

我有一个分位数回归模型,我有兴趣估计.25,.5和.875分位数的效果.我的模型中的系数彼此不同,其方式符合我的模型的实质性实质理论.

下一步是测试一个分位数的特定解释变量的系数是否与另一个分位数的估计系数显着不同.我该如何测试？此外,我还想测试给定分位数的该变量的系数是否与OLS模型中的estimnate显着不同.我怎么做？

我对任何答案感兴趣,虽然我更喜欢一个涉及R的答案.这里有一些测试代码:(注意:这不是我的实际模型或数据,但是一个简单的例子,因为数据在R安装中可用)

data(airquality)
library(quantreg)
summary(rq(Ozone ~ Solar.R + Wind + Temp, tau = c(.25, .5, .75), data = airquality, method = "br"), se = "nid")

tau: [1] 0.25

Coefficients:
            Value     Std. Error t value   Pr(>|t|) 
(Intercept) -69.92874  12.18362   -5.73957   0.00000
Solar.R       0.06220   0.00917    6.77995   0.00000
Wind         -2.63528   0.59364   -4.43918   0.00002
Temp          1.43521   0.14363    9.99260   0.00000

Call: rq(formula = Ozone ~ Solar.R + Wind + Temp, tau = c(0.25, 0.5, 
    0.75), data = airquality, method = "br")

tau: [1] 0.5 …

我正在使用具有许多NA值的数据框.虽然我能够创建一个线性模型,但由于缺少值和缺少指标列,我随后无法将模型的拟合值与原始数据对齐.

这是一个可重复的例子:

library(MASS)
dat <- Aids2
# Add NA's 
dat[floor(runif(100, min = 1, max = nrow(dat))),3] <- NA
# Create a model
model <- lm(death ~ diag + age, data = dat)
# Different Values
length(fitted.values(model))
# 2745
nrow(dat)
# 2843

我正在从Stata转换到R.在Stata中,如果我将因子级别(比如-0和1)标记为(M和F),则0和1将保持不变.此外,在大多数软件(包括Excel和SPSS)中,这是虚拟变量线性回归所必需的.

但是,我注意到R默认因子水平为1,2而不是0,1.我不知道为什么R会这样做,尽管内部(并且正确地)回归假定0和1作为因子变量.我将不胜感激任何帮助.

这是我做的:

尝试#1:

sex<-c(0,1,0,1,1)
sex<-factor(sex,levels = c(1,0),labels = c("F","M"))
str(sex)
Factor w/ 2 levels "F","M": 2 1 2 1 1

似乎因子水平现在重置为1和2.我相信1和2是这里对因子水平的引用.但是,我丢失了原始值,即0和1.

Try2:

sex<-c(0,1,0,1,1)
sex<-factor(sex,levels = c(0,1),labels = c("F","M"))
str(sex)
Factor w/ 2 levels "F","M": 1 2 1 2 2

同上.我的0和1现在是1和2.非常令人惊讶.为什么会这样呢？

Try3 现在,我想知道1和2是否有任何不良影响回归.所以,这就是我做的:

这是我的数据:

> head(data.frame(sassign$total_,sassign$gender))
  sassign.total_ sassign.gender
1            357              M
2            138              M
3            172              F
4            272              F
5            149              F
6            113              F

myfit<-lm(sassign$total_ ~ sassign$gender)

myfit$coefficients
    (Intercept) sassign$genderM 
      200.63522        23.00606  

事实证明,手段是正确的.在运行回归时,R确实使用0和1值作为假人.

我确实在SO上检查了其他线程,但是他们主要讨论R代码如何在不告诉我原因的情况下对变量进 Stata和SPSS通常要求基变量为"0".所以,我想到了这个问题.

我很感激任何想法.

我想在R中进行序数回归,所以我想使用包中的polr函数MASS.首先,我创建一个这样的模型:

model <- polr(labels ~ var1 + var2, Hess = TRUE)  

现在我想使用该模型来预测新病例.我以为那只是:

pred <- predict(model, data = c(newVar1, newVar2))  

然而,似乎预测是以某种方式预测训练集,而不是新数据.当我的训练集是2000个例子时,我的新数据是700个例子.我仍然得到2000个预测标签.

所以我的问题是:如何使用polr新数据进行预测？

作为编写我自己的弹性网络解算器的热身,我正在尝试使用坐标下降来实现足够快的普通最小二乘法.

我相信我已经正确地实现了坐标下降算法,但是当我使用"快速"版本(见下文)时,该算法非常不稳定,输出的回归系数通常在特征数量为6时浮动64位浮点数.与样品数量相比,尺寸适中.

线性回归和OLS

如果b = A*x,其中A是矩阵,未知回归系数的xa向量,y是输出,我想找到最小化的x

|| b - Ax || ^ 2

如果A [j]是A的第j列,而A [-j]是没有列j的A,并且A的列被归一化,所以对于所有j,|| A [j] || ^ 2 = 1,坐标然后是逐步更新

坐标下降:

x[j]  <--  A[j]^T * (b - A[-j] * x[-j])

我正在关注这些注释(第9-10页),但推导是简单的微积分.

它指出,不是一直重新计算A [j] ^ T(b - A [-j]*x [-j]),更快的方法是

快速坐标下降:

x[j]  <--  A[j]^T*r + x[j]

其中总残差r = b - Ax是在环路坐标之外计算的.这些更新规则的等效性来自注意到Ax = A [j]*x [j] + A [-j]*x [-j]并重新排列术语.

我的问题是,虽然第二种方法确实更快,但只要特征数量与样本数量相比不小,它在数值上就会非常不稳定.我想知道是否有人可能会对为什么会出现这种情况有所了解.我应该注意到,第一种更稳定的方法仍然开始不同意更多标准方法,因为特征数量接近样本数量.

朱莉娅代码

以下是两个更新规则的一些Julia代码:

function OLS_builtin(A,b)
    x = A\b
    return(x)
end

function OLS_coord_descent(A,b)    
    N,P = size(A) …

我想测试我的线性回归模型的所有假设是否成立.我手动完成这个,看起来很好.但是,我想用函数gvlma仔细检查.我得到的输出是:

 gvlma(x = m_lag) 

                Value p-value                   Decision
 Global Stat        82.475 0.00000 Assumptions NOT satisfied!
 Skewness           72.378 0.00000 Assumptions NOT satisfied!
 Kurtosis            1.040 0.30778    Assumptions acceptable.
 Link Function       6.029 0.01407 Assumptions NOT satisfied!
 Heteroscedasticity  3.027 0.08187    Assumptions acceptable.

我的问题是:

1. 我如何解释Global Stat

2. 由于违反了这个假设,我现在可以做些什么呢？(与未被接受的其他2个假设相同)

最近，我开始学习sklearn，numpy和pandas，并为多元线性回归函数。我想知道，是否可以进行多元多项式回归？

这是我的多元多项式回归代码，它显示此错误：

in check_consistent_length " samples: %r" % [int(l) for l in lengths])
ValueError: Found input variables with inconsistent numbers of samples: [8, 3]

你知道是什么问题吗？

import numpy as np
import pandas as pd
import xlrd
from sklearn import linear_model
from sklearn.model_selection import train_test_split

def polynomial_prediction_of_future_strenght(input_data, cement, blast_fur_slug,fly_ash,
                                              water, superpl, coarse_aggr, fine_aggr, days):

    variables = prediction_accuracy(input_data)[4]
    results = prediction_accuracy(input_data)[5]

    var_train, var_test, res_train, res_test = train_test_split(variables, results, test_size = 0.3, random_state = 4)

    Poly_Regression = PolynomialFeatures(degree=2)
    poly_var_train = Poly_Regression.fit_transform(var_train)
    poly_var_test = Poly_Regression.fit_transform(var_test) …