Assignment #4

Questions

5-1. 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 = α + β(x_i − z_i) \tag{3} \\ μ_i = α + β_xx_i + β_zz_i \tag{4} \\ \end{align}\]

# In the above list models 2 and 4 can be considered as multiple linear regression since they both have 2 independent variables.

5-2. 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.

# Here  considering f as amount of funding and s as size of laboratory

##\(\mu\)_i = \(\beta\)_ff_i + \(\beta\)_ss_i

5-3. 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). Plot the correlation before analysis, designate and run the ‘quap’ model, show the ‘precis’ table, plot the results and explain.

library(rethinking)
library(dplyr)
library(tibble)
library(tidyr)
n = 100
obesity <- rnorm(n, mean = 0, sd = 1)
soda <- rnorm(n, mean = obesity, sd = 2)
diabetic <- rnorm(n, mean = obesity, sd = 1)
df <- data.frame(obesity, soda, diabetic)
plot(df)

m.1 <- map(
  alist(
    obesity ~ dnorm(mu, sigma),
    mu <- a + bo * soda,
    a ~ dnorm(0, 5),
    bo ~ dnorm(0, 5),
    sigma ~ dunif(0, 5)
  ),
  data = df
)
precis(m.1)

##             mean         sd        5.5%     94.5%
## a     0.07908623 0.09071023 -0.06588623 0.2240587
## bo    0.20346056 0.04041661  0.13886702 0.2680541
## sigma 0.89022197 0.06294626  0.78962168 0.9908223

plot(precis(m.1))

m.2 <- map(
  alist(
    obesity ~ dnorm(mu, sigma),
    mu <- a + bo * soda + bi * diabetic,
    a ~ dnorm(0, 5),
    bo ~ dnorm(0, 5),
    bi ~ dnorm(0, 5),
    sigma ~ dunif(0, 5)
  ),
  data = df
)
precis(m.2)

##             mean         sd        5.5%     94.5%
## a     0.08101354 0.06444440 -0.02198106 0.1840081
## bo    0.17132878 0.02889414  0.12515037 0.2175072
## bi    0.41852972 0.04224249  0.35101806 0.4860414
## sigma 0.63239683 0.04471752  0.56092960 0.7038641

plot(precis(m.2))

5-4. 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. Plot the correlation before analysis, designate and run the ‘quap’ model, show the ‘precis’ table, plot the results and explain.

library(dplyr)
library(tibble)
library(tidyr)
library(rethinking)
N <- 100
rho <- 0.6

sports <- rnorm(n = N, mean = 0, sd = 1)
injury <- rnorm(n = N, mean = rho * sports, sd = sqrt(1 - rho^2))
health <- rnorm(n = N, mean = sports - injury, sd = 1)
d <- data.frame(sports, injury,health )

plot(d)

m5.1 <- quap(
  alist(
    sports ~ dnorm( mu , sigma ) ,
    mu <- a + bN*injury ,
    a ~ dnorm( 0 , 0.2 ) ,
    bN ~ dnorm( 0 , 0.5 ) ,
    sigma ~ dexp( 1 )
  ) ,
  data=d 
)
precis(m5.1)

##              mean         sd       5.5%     94.5%
## a     -0.02343503 0.07846400 -0.1488357 0.1019656
## bN     0.64944184 0.08201257  0.5183699 0.7805138
## sigma  0.85034332 0.05979084  0.7547860 0.9459006

plot(precis(m5.1))

m5.2 <- quap(
  alist(
    sports ~ dnorm( mu , sigma ) ,
    mu <- a + bM*injury,
    a ~ dnorm( 0 , 0.2 ) ,
    bM ~ dnorm( 0 , 0.5 ) ,
    sigma ~ dexp( 1 )
  ) ,
  data=d 
)
precis(m5.2)

##              mean         sd       5.5%     94.5%
## a     -0.02343975 0.07845386 -0.1488242 0.1019447
## bM     0.64952906 0.08200002  0.5184772 0.7805809
## sigma  0.85021347 0.05976801  0.7546926 0.9457343

plot(precis(m5.2))

m5.3 <- quap(
  alist(
    sports ~ dnorm( mu , sigma ) ,
    mu <- a + bN*injury + bM*health,
    a ~ dnorm( 0 , 0.2 ) ,
    bN ~ dnorm( 0 , 0.5 ) ,
    bM ~ dnorm( 0 , 0.5 ) ,
    sigma ~ dexp( 1 )
  ) ,
  data=d 
)
precis(m5.3)

##              mean         sd       5.5%      94.5%
## a     -0.01037494 0.06439402 -0.1132890 0.09253914
## bN     0.78558618 0.06791126  0.6770509 0.89412149
## bM     0.37048407 0.04919221  0.2918654 0.44910272
## sigma  0.67766251 0.04770376  0.6014227 0.75390233

plot(precis(m5.3))

compare(m5.1,m5.2,m5.3)

##          WAIC       SE    dWAIC      dSE    pWAIC       weight
## m5.3 215.2370 18.46153  0.00000       NA 4.299609 1.000000e+00
## m5.1 257.7318 14.18892 42.49480 13.71703 2.627388 5.920677e-10
## m5.2 258.0320 14.24775 42.79502 13.76082 2.757691 5.095420e-10

5-5. 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. Make sure to include the ‘quap’ model, the ‘precis’ table and your explanation of the results.

5-6. In the divorce example, suppose the DAG is: M → A → D. What are the implied conditional independencies of the graph? Are the data consistent with it? (Hint: use the dagitty package)

library(dagitty)

mad_dag <- dagitty("dag{ M -> A -> D}")
impliedConditionalIndependencies(mad_dag)

## D _||_ M | A

equivalentDAGs(mad_dag)

## [[1]]
## dag {
## A
## D
## M
## A -> D
## M -> A
## }
## 
## [[2]]
## dag {
## A
## D
## M
## A -> D
## A -> M
## }
## 
## [[3]]
## dag {
## A
## D
## M
## A -> M
## D -> A
## }

Assignment #4

Gayathri Mutyala

2021-12-06

Chapter 5 - Many Variables and Spurious Waffles

Questions