Блог

Итак, я хочу протестировать все возможные модели линейной регрессии, доступные с 1-5 независимыми переменными и одной из 18 зависимых переменных. У меня есть код, работающий для генерации всех моделей линейной регрессии с 1-й зависимой переменной и 5 независимыми, но я не уверен, как запустить этот код для каждой из 18 зависимых переменных, которые я хочу проверить.

GClist <- data.frame(GC1, GC2, GC3, GC4, GC5, GC6, GC7, GC8, GC9, GC10, GC11, GC12, GC13, GC14, GC15, GC16, GC17, GC18)

До сих пор я составлял список из 18 DV, я также пытался выполнить цикл, используя цикл foreach, а затем пытался проанализировать список, содержащий 18 имен DV.

Я пытался использовать:

 for(value in GClist)
{
the code below
}
// but then I did not manage to make it work and include the "value" in the code

// I also tried to use foreach, but I using df[j] and having df containing all my 18 dependent variables did not seem to work.


foreach (j=1:18) %do% {the code below }
  

В любом случае, код, который работает, это:

 df1 <- data.frame(GC1, D1, D2, D3, D4, D5)
library(foreach)

#create the linear models list 

xcomb <- foreach(i=1:5, .combine=c) %do% {combn(names(df1)[-1], i, simplify=FALSE) }
formlist <- lapply(xcomb, function(l) formula(paste(names(df1), paste(l, collapse=" "), sep="~")))```

 # get the p value for each model
models.p <- sapply(formlist, function (i)  {
  f <- summary(lm(i))$fstatistic
  p <- pf(f[1],f[2],f[3],lower.tail=F)
  attributes(p) <- NULL
  return(p)
})
# R squared for each model
models.r.sq <- sapply(formlist, function(i) summary(lm(i))$r.squared)
# adjusted R squared for each model
models.adj.r.sq <- sapply(formlist, function(i) summary(lm(i))$adj.r.squared)
# MSEp squared for each model
models.MSEp <- sapply(formlist, function(i) anova(lm(i))['Mean Sq']['Residuals',])

# Full model MSE for each linear model
MSE <- anova(lm(formlist[[length(formlist)]]))['Mean Sq']['Residuals',]

# Mallow's Cp - skipped for now 

df.model.eval <- data.frame(model=as.character(formlist), pF=models.p,
                            r.sq=models.r.sq, adj.r.sq=models.adj.r.sq, MSEp=models.MSEp)
  

Как я могу запустить этот код для каждой из 18 зависимых переменных, чтобы я мог собрать всю информацию в df.model.eval? Прямо сейчас у меня есть все, что мне нужно знать о моделях, которые используют GC1 в качестве зависимой переменной. Моя цель — просмотреть все модели (от GC1 до GC18) и выделить все, которые имеют статистическую значимость, и те, которые этого не делают.

Мне удалось это сделать, создав фрейм данных с возможными DVS;

GC <- data.frame(GC1, GC2, GC3, GC4, GC5, GC6, GC7, GC8, GC9, GC10, GC11, GC12, GC13, GC14, GC15, GC16, GC17, GC18)

Затем, используя код Джорджа Донтаса, найденный в https://stats.stackexchange.com/a/6862 , Я готовлю свой первый фрейм данных для генерации моделей линейной регрессии.

 library(foreach)

finallist <- list()
xcomb <- foreach(i=1:5, .combine=c) %do% {combn(names(df1)[-1], i, simplify=FALSE) }
  

Я заменяю первый элемент моего фрейма данных здесь:

 for(j in c(names(GC)))
{
  names(df1)[1] <- j #switch DVs
  #George's code to create a list with all the possible combinations of the DV and the IVs
  formlist <- lapply(xcomb, function(l) formula(paste(names(df1), paste(l, collapse=" "), sep="~"))) 
  #append the list of linear regression models for each DV to the main list
  finallist <- append(finalist, formlist)
}
  

Теперь у меня есть список всех моделей линейной регрессии, которые я затем могу использовать для оценки и интерпретации соответствующих статистических показателей, с
как.character(финалист)

Мой окончательный код:

 df1 <- data.frame(GC1, D1, D2, D3, D4, D5)
library(foreach)
library(rlist)
finallist <- list()
xcomb <- foreach(i=1:5, .combine=c) %do% {combn(names(df1)[-1], i, simplify=FALSE) }

#iterate through DV names
for(j in c(names(GC)))
{
  names(df1)[1] <- j #switch DVs (I know 1st step replaces same DV name)
  #create a list with all the possible combinations of the DV and the IVs
  formlist <- lapply(xcomb, function(l) formula(paste(names(df1), paste(l, collapse=" "), sep="~"))) 
  #append the list of linear regression models for each DV to the main list
  finallist <- append(finalist, formlist)
}

# get the p value for each model
models.p <- sapply(formlist, function (i)  {
  f <- summary(lm(i))$fstatistic
  p <- pf(f[1],f[2],f[3],lower.tail=F)
  attributes(p) <- NULL
  return(p)
})
# R squared for each model
models.r.sq <- sapply(finallist, function(i) summary(lm(i))$r.squared)
# adjusted R squared for each model
models.adj.r.sq <- sapply(finallist, function(i) summary(lm(i))$adj.r.squared)
# MSEp squared for each model
models.MSEp <- sapply(finallist, function(i) anova(lm(i))['Mean Sq']['Residuals',])

# Full model MSE for each linear model
MSE <- anova(lm(finallist[[length(finallist)]]))['Mean Sq']['Residuals',]

# Mallow's Cp - skipped for now 

df.model.eval <- data.frame(model=as.character(finallist), pF=models.p,
                            r.sq=models.r.sq, adj.r.sq=models.adj.r.sq, MSEp=models.MSEp)