Assignment #4

Questions

5E1. Which of the linear models below are multiple linear regressions? \[\begin{align} {μ_i = α + βx_i} \tag{1}\\ μ_i = β_xx_i + β_zz_i \tag{2} \\ μ_i = β_xx_i + β_zz_i \tag{3} \\ μ_i = α + β(x_i − z_i) \tag{4} \\ μ_i = α + β_xx_i + β_zz_i \tag{5} \\ \end{align}\]

#2, 3, 4, 5

5E2. Write down a multiple regression to evaluate the claim: Animal diversity is linearly related to latitude, but only after controlling for plant diversity. You just need to write down the model definition.

#μi = α + βAAi + βPPi

#In the above regression, A is the animal diversity and P is the plant diversity.

5E3. Write down a multiple regression to evaluate the claim: Neither amount of funding nor size of laboratory is by itself a good predictor of time to PhD degree; but together these variables are both positively associated with time to degree. Write down the model definition and indicate which side of zero each slope parameter should be on.

#μi=α+βFFi+βSSi

#In the above regression, F is the amount of funding and S is the size of laboratory. Both slope parameters are on the right side of zero, positive.

5E4. Suppose you have a single categorical predictor with 4 levels (unique values), labeled A, B, C and D. Let Ai be an indicator variable that is 1 where case i is in category A. Also suppose Bi, Ci, and Di for the other categories. Now which of the following linear models are inferentially equivalent ways to include the categorical variable in a regression? Models are inferentially equivalent when it’s possible to compute one posterior distribution from the posterior distribution of another model. \[\begin{align} μ_i = α + β_AA_i + β_BB_i + β_DD_i \tag{1} \\ μ_i = α + β_AA_i + β_BB_i + β_CC_i + β_DD_i \tag{2} \\ μ_i = α + β_BB_i + β_CC_i + β_DD_i \tag{3} \\ μ_i = α_AA_i + α_BB_i + α_CC_i + α_DD_i \tag{4} \\ μ_i = α_A(1 − B_i − C_i − D_i) + α_BB_i + α_CC_i + α_DD_i \tag{5} \\ \end{align}\]

#1, 3, 4, 5

5M1. Invent your own example of a spurious correlation. An outcome variable should be correlated with both predictor variables. But when both predictors are entered in the same model, the correlation between the outcome and one of the predictors should mostly vanish (or at least be greatly reduced).

#An example I came up with is the spurious correlation between the sale of houses on the market and pas average price when the pandemic situation is taken into account. Assume sale is the amount of house sold (in 10k), avg_price the pas average price (in 10k), c19 the case of covid_19 cases confirmed ( in 10k).

N = 1000
c19 <- rnorm(n = N, mean = 10, sd = 3)
sale <- rnorm(n = N, mean = c19, sd = 2)
avg_price <- rnorm(n = N, mean = c19, sd = 1)
df <- data.frame(c19, sale, avg_price)
pairs(df)

# sale ~ c19
m <- map(
  alist(
    sale ~ dnorm(mu, sigma),
    mu <- a + bo * c19,
    a ~ dnorm(0, 5),
    bo ~ dnorm(0, 5),
    sigma ~ dunif(0, 5)
  ),
  data = df
)
precis(m)

##            mean         sd       5.5%     94.5%
## a     0.2974094 0.23161261 -0.0727523 0.6675711
## bo    0.9679405 0.02166465  0.9333162 1.0025648
## sigma 2.1108898 0.04719843  2.0354576 2.1863220

#sale ~ c19 + past average price. From the result, it correlation between c19 and sale seems spurious.
m2 <- map(
  alist(
    sale ~ dnorm(mu, sigma),
    mu <- a + bo * c19 + bi * avg_price,
    a ~ dnorm(0, 5),
    bo ~ dnorm(0, 5),
    bi ~ dnorm(0, 5),
    sigma ~ dunif(0, 5)
  ),
  data = df
)
precis(m2)

##             mean         sd        5.5%     94.5%
## a     0.30829275 0.23177776 -0.06213288 0.6787184
## bo    0.90008773 0.07141421  0.78595403 1.0142214
## bi    0.06694238 0.06713723 -0.04035588 0.1742406
## sigma 2.10997451 0.04718233  2.03456803 2.1853810

5M2. Invent your own example of a masked relationship. An outcome variable should be correlated with both predictor variables, but in opposite directions. And the two predictor variables should be correlated with one another.

#An example I invented is exploring the amount of weight loss and the consumed and burnt calories.

N = 100
r = 0.5
consumed = rnorm(n = N, mean = 0, sd = 1)
burnt = rnorm(n = N, mean = r*consumed, sd = sqrt(1 - r^2))
weight = rnorm(n = N, mean = consumed - burnt, sd = 1)
df2 <- data.frame(consumed, burnt, weight)
pairs(df2)

# weight ~ consumed
m3 <- map(
  alist(
    weight ~ dnorm(mu, sigma),
    mu <- a + bo * consumed,
    a ~ dnorm(0, 5),
    bo ~ dnorm(0, 5),
    sigma ~ dunif(0, 5)
  ),
  data = df2
)
precis(m3)

##            mean         sd        5.5%     94.5%
## a     0.1363377 0.12195686 -0.05857286 0.3312484
## bo    0.5693075 0.11891014  0.37926610 0.7593488
## sigma 1.2129691 0.08576959  1.07589274 1.3500455

#weight ~ burnt
m4 <- map(
  alist(
    weight ~ dnorm(mu, sigma),
    mu <- a + bi * burnt,
    a ~ dnorm(0, 5),
    bi ~ dnorm(0, 5),
    sigma ~ dunif(0, 5)
  ),
  data = df2
)
precis(m4)

##             mean         sd        5.5%      94.5%
## a      0.1664104 0.12861395 -0.03913957  0.3719603
## bi    -0.3938134 0.12479534 -0.59326041 -0.1943663
## sigma  1.2825104 0.09068702  1.13757504  1.4274458

#weight ~ consumed + burnt
m5 <- map(
  alist(
    weight ~ dnorm(mu, sigma),
    mu <- a + bo * consumed + bi * burnt,
    a ~ dnorm(0, 5),
    bo ~ dnorm(0, 5),
    bi ~ dnorm(0, 5),
    sigma ~ dunif(0, 5)
  ),
  data = df2
)
precis(m5)

##              mean         sd       5.5%      94.5%
## a     -0.00164302 0.08820577 -0.1426129  0.1393268
## bo     1.09888856 0.10061443  0.9380873  1.2596899
## bi    -0.97929153 0.09987154 -1.1389055 -0.8196775
## sigma  0.86591857 0.06122926  0.7680624  0.9637747

5M3. It is sometimes observed that the best predictor of fire risk is the presence of firefighters— States and localities with many firefighters also have more fires. Presumably firefighters do not cause fires. Nevertheless, this is not a spurious correlation. Instead fires cause firefighters. Consider the same reversal of causal inference in the context of the divorce and marriage data. How might a high divorce rate cause a higher marriage rate? Can you think of a way to evaluate this relationship, using multiple regression?

#The more divorce occurs, the more single people are available in the marriage market. With more people in the market, there is a higher chance that they find their Mr./Ms. Right to marry to. Maybe we can use the rate of people who have married more than once into the multiple regression and see if its influence takes over that of the divorce rate.

5M4. In the divorce data, States with high numbers of members of the Church of Jesus Christ of Latter-day Saints (LDS) have much lower divorce rates than the regression models expected. Find a list of LDS population by State and use those numbers as a predictor variable, predicting divorce rate using marriage rate, median age at marriage, and percent LDS population (possibly standardized). You may want to consider transformations of the raw percent LDS variable.

data(WaffleDivorce)
wd <- WaffleDivorce

wd$LDS <- c(0.0077, 0.0453, 0.0610, 0.0104, 0.0194, 0.0270, 0.0044, 0.0057, 0.0041, 0.0075, 0.0082, 0.0520, 0.2623, 0.0045, 0.0067, 0.0090, 0.0130, 0.0079, 0.0064, 0.0082, 0.0072, 0.0040, 0.0045, 0.0059, 0.0073, 0.0116, 0.0480, 0.0130, 0.0065, 0.0037, 0.0333, 0.0041, 0.0084, 0.0149, 0.0053, 0.0122, 0.0372, 0.0040, 0.0039, 0.0081, 0.0122, 0.0076, 0.0125, 0.6739, 0.0074, 0.0113, 0.0390, 0.0093, 0.0046, 0.1161)
wd$logLDS = log(wd$LDS)
wd$logLDS.s = (wd$logLDS - mean(wd$logLDS)) / sd(wd$logLDS)

simplehist(wd$LDS)

simplehist(wd$logLDS)

simplehist(wd$logLDS.s)

#m_std <- (wd$Marriage - mean(wd$Marriage)) / sd(wd$Marriage)
#m_med <- (wd$MedianAgeMarriage - mean(wd$MedianAgeMarriage)) / sd(wd$MedianAgeMarriage)
#per_lds <- (lds - mean(lds)) / sd(lds)

m6 <-map(
  alist(
    Divorce ~ dnorm(mu, sigma),
    mu<- a + bm * Marriage + ba * MedianAgeMarriage + bl * logLDS.s,
    a ~ dnorm(10, 20),
    bm ~ dnorm(0, 10),
    ba ~ dnorm(0, 10),
    bl ~ dnorm(0, 10),
    sigma ~ dunif(0, 5)
  ),
  data = wd
)
precis(m6)

##              mean         sd        5.5%      94.5%
## a     35.43685973 6.77531065 24.60860473 46.2651147
## bm     0.05341282 0.08261844 -0.07862741  0.1854530
## ba    -1.02961628 0.22469488 -1.38872210 -0.6705105
## bl    -0.60807737 0.29057047 -1.07246511 -0.1436896
## sigma  1.37872564 0.13838911  1.15755312  1.5998982

#The coefficient of Marriage is positive; ba and bl are negative and their intervals do not include 0. That means states with lower median age at marriage or lower percentages of LDS have higher divorce rate.

5M5. One way to reason through multiple causation hypotheses is to imagine detailed mechanisms through which predictor variables may influence outcomes. For example, it is sometimes argued that the price of gasoline (predictor variable) is positively associated with lower obesity rates (outcome variable). However, there are at least two important mechanisms by which the price of gas could reduce obesity. First, it could lead to less driving and therefore more exercise. Second, it could lead to less driving, which leads to less eating out, which leads to less consumption of huge restaurant meals. Can you outline one or more multiple regressions that address these two mechanisms? Assume you can have any predictor data you need.

#μ_i = α + β_GGasoline_i + β_EExercise_i + β_RRestaurant_i
#Gasoline: price of gasoline
#Exercise: amount of exercise (per day, per week, bi-weekly......)
#Restaurant: frequency of going eating out (per day, per week, bi-weekly......)

Assignment #4

Shidong Li

2020-12-05

Chapter 5 - Many Variables and Spurious Waffles

Questions