CODEWITHSUNDEEP
rloopsmatrixindexing
coding_heart
coding_heart

Reputation: 1295

R Compare all columns in a matrix against each in loop

I have a matrix of variables and I am trying to run a loop that compares the difference between all of the variables in a regression such that a matrix is produced and filled in with the difference. Below is some simulation code to get at the problem. I would like to produce a matrix that compares x_1, x_2, and x_3 to produce a 3x3 matrix that is symmetric about the diagonal, which should all be zeros. 

y <- sample(seq(1:4), 100, replace = TRUE)
x_1 <- sample(seq(1:2), 100, replace = TRUE)
x_2 <- sample(seq(1:4), 100, replace = TRUE)
x_3 <- sample(seq(1:4), 100, replace = TRUE)

frame <- cbind(x_1, x_2, x_3)
dif <- matrix(NA, ncol = 3, nrow = 3)

for(i in 1:3){
    model_1 <- lm(y ~ frame[,i])
    model_2 <- lm(y ~ frame[,i])
    dif[i]<- (model_2$coef[2] - model_1$coef[2])
}

I'm confused on how to index the loop and refer to the matrix of x's to produce a 3x3 table with the results - any help would be much appreciated.

Upvotes: 3

Views: 2411

Answers (5)

Tyler Rinker
Tyler Rinker

Reputation: 109844

I may not understand what you're after but here's a crack using a version of outer that takes multiple vectors.

library(qdap)
FUN <- function(x1, x2, y)lm(y ~ x1)$coef[2] - lm(y ~ x2)$coef[2]
v_outer(list(x_1, x_2, x_3), FUN, y = y)

##       X1     X2     X3
## X1 0.000 -0.311 -0.079
## X2 0.311  0.000  0.232
## X3 0.079 -0.232  0.000

Upvotes: 0

agstudy
agstudy

Reputation: 121568

If I understand you try to compare models by comparing their coefficients. One idea is to use meifly package.

First I generate your data :

set.seed(1)
frame <- matrix(sample(1:4,3*100,rep=TRUE),ncol=3)
y <- sample(seq(1:4), 100, replace = TRUE)

Then I use fitbest which uses the leaps package to very rapidly find the n best models for a given number of variables.

library(meifly)
library(reshape2)
library(ggplot2)
## we look only on models with one variable
res <- fitbest(y~.,as.data.frame(frame),nvmax=1)
## get coefficients
res.coef <- coef(res)
## remove zero models 
res.coef[res.coef == 0] <- NA
res.coef <- na.omit(res.coef)

Now for each model, we have a summary of coefficients. For each variable we have the following information :

  1. Raw coefficient
  2. t value and absolute t-value
  3. Standardised coefficient.

And res.coef looks like this :

res.coef
     model observ         raw          t      abst         std
m1v1     1     V1 -0.12884211 -1.2438295 1.2438295 -0.14110975
m2v2     2     V3  0.09258638  0.8922776 0.8922776  0.10161095
m3v3     3     V2  0.01534989  0.1420060 0.1420060  0.01623527

One way to compare the models, is to plot a Scatterplot of all vs variable

colnames(res.coef)[colnames(res.coef)== "variable"] <- "observ"
dat <- melt(res.coef)
ggplot(dat) +
   geom_point(aes(observ,value,color=variable),size=5) +
   theme_bw()

enter image description here

Upvotes: 0

Doctor Dan
Doctor Dan

Reputation: 771

Try this:

model <- list()
for(i in 1:3) { 
    model[[i]] <- lm(y~frame[,i])
}

dif<-sapply( 1:3, function(i) { sapply(1:3, function(j) { model[[i]]$coef[2] - model[[j]]$coef[2] } ) } )

The matrix will be antisymmetric, i.e., dif[i,j] = -dif[j,i]

Upvotes: 0

IRTFM
IRTFM

Reputation: 263331

 vcoef <- numeric(3)
 for(i in 1:3) { 
     vcoef[i] <- coef( lm(y~frame[,i]))[2]
               }

outer(vcoef, vcoef, "-")
#----------
          [,1]        [,2]        [,3]
[1,] 0.0000000 -0.15208933 -0.17302592
[2,] 0.1520893  0.00000000 -0.02093659
[3,] 0.1730259  0.02093659  0.00000000

If you didn't want the redundant information you could get all the pairwise differences with combn:

> combcos  <- combn(vcoef,2)
> combcos[1, ] -combcos[2, ]
[1] -0.15208933 -0.17302592 -0.02093659

Upvotes: 2

Simon O&#39;Hanlon
Simon O&#39;Hanlon

Reputation: 59970

I prefer the eval and parse route and like @Tyler I like base:::outer...

#  Make your data into a data.frame
df <- data.frame( y , x_1 , x_2 , x_3 )

#  The variables we want to test
x <- c("x_1","x_2","x_3") 

#  Make the text for each model to parse and evalaute
mods <- paste0( "lm( y ~ " , x , " , data = df )" )

#  Evaluate the lm for each variable
coefs <- unlist( lapply( mods , function(x) eval(parse(text=x))$coef[2] ) )
#       x_1         x_2         x_3 
# -0.52140856  0.04662379  0.08694344 

#  Combine the results with outer
outer( coefs , coefs ,  "-")
#         x_1         x_2         x_3
# x_1 0.0000000 -0.56803236 -0.60835201
# x_2 0.5680324  0.00000000 -0.04031965
# x_3 0.6083520  0.04031965  0.00000000

Upvotes: 0

Related Questions