6

I'm trying to conduct a multivariate multiple regression analysis. Fortunately, I've found an excellent page demonstrating how to do so in Stata:

http://www.ats.ucla.edu/stat/stata/dae/mvreg.htm

The issue is that I'm using R, while I've figured out how to do the basics of running a multivariate multiple regression model in R, I'm still not sure how to see if coefficients are different for each dependent variable (as shown in the link). Does anyone know how to compute this analysis in R? Seems like seeing if the same independent variable exerts different effects on each dependent variable is an incredibly useful tool, and I'd love to be able to do it!

UPDATE: Here is a reproducible example of what I've done with my own data thus far:

# data
data(mtcars)

# fitting a multivariate multiple regression where mpg and cyl are predicted by wt and hp
car.mod <- lm(cbind(mpg,cyl) ~ wt + hp,data=mtcars)

# see if there is a multivariate effect of wt and hp
summary(manova(car.mod),test="W")

# Get coefficients for each dependent variable
summary(car.mod)

What I want to know in this example is how I can test the equivalence of "wt" on both "mpg" and "cyl". Apparently, this is possible in Stata using the test command.

Nick Cox
  • 35,529
  • 6
  • 31
  • 47
user2917781
  • 273
  • 2
  • 10
  • 1
    Check out the [r tag info page](http://stackoverflow.com/tags/r/info) for lots of R-specific resources, some of which cover this topic. What do you have so far, though? You need to add [a reproducible example](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example#5963610). – alistaire Dec 03 '16 at 21:42
  • 1
    Ok, this is done! – user2917781 Dec 03 '16 at 21:51

2 Answers2

2

AFAIK, there is no package that does this, so I would do the manual solution. The manual solution is

z = (b1 - b2) / (b1^2 + b2^2)^(1/2).

Here is the (sloppy) code. There may be a more elegant solution to extract the coefficients and standard errors.

# extract coefficients and SEs
coef_mpg <- coef(summary(car.mod))[[1]]
coef_cyl <- coef(summary(car.mod))[[2]]

# test equality of wt coefficients
z_wt <- (coef_mpg[2, 1] - coef_cyl[2, 1]) / (coef_mpg[2, 2]^2 + coef_cyl[2, 2]^2)^(1/2)
p_wt <- 2*pnorm(-abs(z_wt)) 
p_wt

But I would feel better about a bootstrap solution, which makes fewer assumptions.

require(boot)
b_b <- function(data=mtcars, indices) {
  d <- data[indices, ]
  model_mpg <- lm(mpg ~ wt + hp, data=d)
  model_cyl <- lm(cyl ~ wt + hp, data=d)
  return(coef(model_mpg)["wt"] - coef(model_cyl)["wt"])
}
results <- boot(data=mtcars, statistic=b_b, R=1000)
results
Richard Herron
  • 9,760
  • 12
  • 69
  • 116
0

So here is an example assuming you have two outcomes and multiple predictors. The formula was suggested by Clifford Clogg et al. (1995) and cited by Ray Paternoster et al. (1998).

Below is some data and script in the R language to illustrate the process by @Richard Herron.

df1 = data.frame(
  estimate = c(15.2418519, 2.2215987, 0.3889724, 0.5289710),
  `std.error` = c(1.0958919, 0.2487793, 0.1973446, 0.1639074),
  row.names = c('(Intercept)', 'psychoticism', 'extraversion', 'neuroticism')
  ); df1

df2 = data.frame(
  estimate = c(17.2373874, 0.8350460, -0.3714803, 1.0382513),
  `std.error` = c(1.0987151, 0.2494201, 0.1978530, 0.1643297),
  row.names = c('(Intercept)', 'psychoticism', 'extraversion', 'neuroticism')
); df2

So in the data examples above, we are assuming that you have a multivariate multiple regression model. df1 and df2 relates to two outcome variables and their relevant predictors (which are the same).

The script below iterates over the relevant data, and outputs the comparisons of the coefficients in a dataframe. I have commented the script, so it's easy to follow along and see what is happening.

compare_coefs <- function(.data1, .data2){
  
  # import map_dbl() and pluck() functions from the purrr library
  import::here(map_dbl, pluck, .from = purrr)

  # extract the coef and std error from df1
  b1 = map_dbl(.data1[-1, 1], pluck)
  se1 = map_dbl(.data1[-1, 2], pluck)

  # extract the coef and std error from df2
  b2 = map_dbl(.data2[-1, 1], pluck)
  se2 = map_dbl(.data2[-1, 2], pluck)
  
  # calculate the difference per Clogg et al (1995) formula 
  # as cited by Paternoster et al. (1998)
  b = b1 - b2
  s1 = se1^2
  s2 = se2^2
  sc = s1 + s2
  v = b / sqrt(sc)
  
  data.frame(diff=b, zdiff=v, `p-value`=format(2*pnorm(-abs(v)), scientific=FALSE))
  
  }

# call the function
compare_coefs(df1, df2)
GSA
  • 751
  • 8
  • 12