2

假设我有一个由以下代码生成的系数列表,并且我希望它采用下面评论中显示的形式。有什么直接的方法吗?

library(leaps)

set.seed(1)

N = 1000
P = 20

x=matrix(rnorm(N*P),N,P)
eps=rnorm(N)

beta = sample(c(0,1), P, replace=T)
y = x %*% beta + eps

regfit.full= regsubsets(y~., data=data.frame(x=x[train,], y=y[train]), nvmax=20)
coefi = coef(regfit.full, id=3)

#Output:
#> coefi
#(Intercept)         x.5        x.10        x.20
# 0.03730904  1.39039580  1.68618982  1.15607983

# How do I generate from coef a vector of the form 
# transformed = [0 0 0 0 1.39039580 0 0 0 0 1.68618982 0...0 1.15607983]
4

3 回答 3

1

尝试类似:

# generate colnames with paste, since x is matrix without column names
# and compare with names in coefficient vector
v <- paste("x", seq_len(ncol(x)), sep=".") %in% names(coefi)
coef_full <- numeric(ncol(x))
coef_full[v] <- coefi[-1] # remove intercept
coef_full
于 2013-11-06T09:57:15.563 回答
1

我不确定train您的特定示例的向量是什么,但这是我自己的一种可能性train

set.seed(2)
train <- sample(length(y), length(y)*0.5)

regfit.full= regsubsets(y~., data=data.frame(x=x[train,], y=y[train]), nvmax=20)
coefi = coef(regfit.full, id=3)

#> coefi
#(Intercept)         x.6         x.7        x.16 
# 0.09013856  1.10080511  0.97903517  1.37892870 


vars <- paste("x", seq(ncol(x)), sep=".")
res <- rep(0, length(vars))
put <- which(vars %in% names(coefi))
take <- which(names(coefi) %in% vars)
res[put] <- coefi[take]
res
#> res
# [1] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 1.1008051 0.9790352
# [8] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
#[15] 0.0000000 1.3789287 0.0000000 0.0000000 0.0000000 0.0000000
于 2013-11-06T10:06:24.517 回答
0

这是另一种方法

> pos <- na.omit(as.numeric(gsub("[[:punct:]]*[[:alpha:]]*", "" ,names(coefi))))
> transformed  <- rep(0, pos[length(pos)])
> transformed[pos] <- coefi[-1]
> 
> coefi       # original
(Intercept)         x.2         x.7        x.12 
 0.08130501  1.16965068  1.24027061  1.22437825 
> transformed # desired
 [1] 0.000000 1.169651 0.000000 0.000000 0.000000 0.000000 1.240271 0.000000 0.000000 0.000000 0.000000 1.224378
于 2013-11-06T10:34:01.427 回答