Lm над всеми возможными попарно комбинациями столбцов двух матриц

На данный момент я работаю над проблемой в R и застрял. Я искал в разных справочных списках помощь, но ничего не мог найти - но извини, если я что-то пропустил. Ниже приведен фиктивный пример моей проблемы. Я буду продолжать работать над этим, но любая помощь будет очень признательна.

Спасибо заранее за ваше время.

У меня есть матрица переменных ответа:

p<-matrix(c(rnorm(120,1),
 rnorm(120,1),
 rnorm(120,1)),
 120,3)

и две матрицы ковариатов:

g<-matrix(c(rep(1:3, each=40),
 rep(3:1, each=40),
 rep(1:3, 40)),
 120,3)

m<-matrix(c(rep(1:2, 60),
 rep(2:1, 60),
 rep(1:2, each=60)),
 120,3)

Для всех комбинаций столбцов ковариационных матриц g и m я хочу запустить эти две модели:

test <- function(uniq_m, uniq_g, p = p) { 
 full <- lm(p ~ factor(uniq_m) * factor(uniq_g))
 null <- lm(p ~ factor(uniq_m) + factor(uniq_g))
 return(list('f'=full, 'n'=null))
}

Поэтому я хочу протестировать взаимодействие между столбцом 1 m и столбцом 1 g, затем столбцом 2 m и столбцом 1 g, затем столбцом 2 m и столбцом 2 g... и так далее через все возможные попарно взаимодействия. Каждый раз переменная ответа одинакова и представляет собой матрицу, содержащую несколько столбцов.

Пока что я могу сделать это для одной комбинации столбцов:

test_1 <- test(m[ ,1], g[ ,1], p)

И я могу также запустить модель по всем столбцам m и одной coloumn из g:

test_2 <- apply(m, 2, function(uniq_m) {
 test(uniq_m, g[ ,1], p = p)
})

Затем я могу получить статистику F для каждой переменной ответа каждой модели:

sapply(summary(test_2[[1]]$f), function(x) x$fstatistic)
sapply(summary(test_2[[1]]$n), function(x) x$fstatistic)

И я могу сравнить модели для каждой переменной ответа с помощью F-теста:

d1<-colSums(matrix(residuals(test_2[[1]]$n),nrow(g),ncol(p))^2)
d2<-colSums(matrix(residuals(test_2[[2]]$f),nrow(g),ncol(p))^2)
F<-((d1-d2) / (d2/114))

Мой вопрос: как запустить lm-модели по всем комбинациям столбцов из m и g-матрицы и получить F-статистику?

Хотя это примерный пример, реальный анализ будет иметь матрицу ответов, которая составляет 700 х 8000, а матрицы ковариации будут 700 х 4000 и 700 х 100, поэтому мне нужно что-то как можно быстрее.

1 ответ

Надеюсь, это поможет, это какой-то код, который мой друг поделился со мной. Возможно, это не совсем то, что вам нужно, но может привести вас в правильном направлении (хотя, учитывая, что это на 9 месяцев позже, чем вы просили, это может быть бесполезно для вас конкретно):

#### this first function models the correlation and fixes the text size based on the strength of the correlation
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
 usr <- par("usr"); on.exit(par(usr))
 par(usr = c(0, 1, 0, 1))
 r <- abs(cor(x, y))
 txt <- format(c(r, 0.123456789), digits = digits)[1]
 txt <- paste0(prefix, txt)
 if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
 text(0.5, 0.5, txt, cex = cex.cor * r)
}

##### this function places a histogram of your data on the diagonal
panel.hist <- function(x, ...)
{
 usr <- par("usr"); on.exit(par(usr))
 par(usr = c(usr[1:2], 0, 1.5) )
 h <- hist(x, plot = FALSE)
 breaks <- h$breaks; nB <- length(breaks)
 y <- h$counts; y <- y/max(y)
 rect(breaks[-nB], 0, breaks[-1], y, col="cyan", ...)
}

### read in Fishers famous iris dataset for our example
data(iris)
head(iris)

library(corrgram)
##corrgram also gives you some nice panel options to use in pairs, but you dont necesarily need them
##e.g. panel.ellipse, panel.pie, panel.conf
library(asbio)
##asbio offers more panel options, such as a linear regression (panel.lm) etc

### run pairs() on your data
### set upper panel to panel.cor (the function we just wrote), and diagonal to panel.hist
### do what you like for the lower, add a smoother line isnt very informative
pairs(~ Sepal.Length + Sepal.Width + Petal.Length, data=iris, lower.panel=panel.lm, upper.panel=panel.cor, diag.panel = panel.hist, main="pair plots of variables")

Весь кредит Джеймсу Китингу.

licensed under cc by-sa 3.0 with attribution.