Вы, вероятно, захотите исключить что-либо о наивысшем полиноме, где p <= .7
(следует сохранить более низкие степени). Предположим, вы знаете, что делаете, вы могли бы написать функцию degAna()
, которая анализирует степени каждого полинома и применяет ее к матрице коэффициентов, полученной summary
.
REG19 <- lm(R10 ~ poly(M1, 3) + poly(M2, 3) + poly(M3, 3) + poly(M4, 3) +
poly(M5, 3) + poly(M6, 3) + poly(M7, 3) + poly(M8, 3) +
poly(M9, 3) + poly(M10, 3), WorkData)
rr <- summary(REG19)$coefficients
Функция, которая определяет наивысшую степень с p <= .7
:
degAna <- function(d) {
out <- as.matrix(rr[grep(paste0(")", d), rownames(rr)), "Pr(>|t|)"] <= .7)
dimnames(out) <- list(c(gsub("^.*\\((.*)\\,.+", "\\1", rownames(out))), d)
return(out)
}
lapply
degAna
в матрицу коэффициентов:
dM <- do.call(cbind, lapply(1:3, degAna)) # max. degree always 3 as in example
# 1 2 3
# M1 TRUE TRUE TRUE
# M2 TRUE TRUE TRUE
# M3 FALSE TRUE TRUE
# M4 TRUE TRUE TRUE
# M5 TRUE TRUE TRUE
# M6 TRUE FALSE TRUE
# M7 TRUE FALSE FALSE
# M8 TRUE TRUE TRUE
# M9 TRUE TRUE FALSE
# M10 TRUE FALSE TRUE
Теперь нам нужна последняя степень многочленов, где p <= .7
:
tM <- apply(dM, 1, function(x) max(which(x != 0)))
tM <- tM[tM > 0] # excludes polynomes where every p < .7
# M1 M2 M3 M4 M5 M6 M7 M8 M9 M10
# 3 3 3 3 3 3 1 3 2 3
(Обратите внимание, что apply
выдаст предупреждение, если многочлен полностью имеет p <= .7
, т.е. строка полностью соответствует FALSE
. Поскольку мы выбрасываем их в следующей строке, мы можем игнорировать предупреждение с помощью apply(dM, 1, function(x) suppressWarnings(max(which(x != 0))))
.)
Используя эту информацию, мы можем составить новую формулу с reformulate
,
terms.new <- paste0("poly(", names(tM), ", ", tM, ")")
FO <- reformulate(terms.new, response="R10")
# R10 ~ poly(M1, 3) + poly(M2, 3) + poly(M3, 3) + poly(M4, 3) +
# poly(M5, 3) + poly(M6, 3) + poly(M7, 1) + poly(M8, 3) + poly(M9,
# 2) + poly(M10, 3)
с которой мы, наконец, можем сделать желаемую сокращенную регрессию.
REG19.2 <- lm(FO, WorkData)
n <- length(resid(REG19.2))
REG20.2 <- step(REG19.2, direction="backward", k=log(n))
# [...]
Моделированные данные
set.seed(42)
M1 <- rnorm(1e3)
M2 <- rnorm(1e3)
M3 <- rnorm(1e3)
M4 <- rnorm(1e3)
M5 <- rnorm(1e3)
M6 <- rnorm(1e3)
M7 <- rnorm(1e3)
M8 <- rnorm(1e3)
M9 <- rnorm(1e3)
M10 <- rnorm(1e3)
R10 <- 6 + 5*M1^3 + 4.5*M2^3 + 4*M3^2 + 3.5*M4 + 3*M5 + 2.5*M6 + 2*M7 +
.5*rnorm(1e3, 1, sd=20)
WorkData <- data.frame(M1, M2, M3, M4, M5, M6, M7, M8, M9, M10, R10)
person
jay.sf
schedule
16.06.2019
dput(head(df,n))
? Выберитеn
, поскольку может быть достаточно для воспроизводимости. - person NelsonGon   schedule 16.06.2019coef(summary(lm("Petal.Length~Sepal.Length+Petal.Width",data = iris))) %>% as_tibble() %>% filter(bac_tick_here_Pr(>|t|) another_back_tick_here < 0.00000002)
? - person NelsonGon   schedule 16.06.2019