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}\]

# tag{2}, {3},{4} are linear regressions

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.

#μ_i = β_fF_i + β_sS_i
#F is amount of funding and S is size of laboratory.

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.

#An example of a spurious correlation between the score on a standardized achievement test and study hours put in that vanishes when happiness is included.
N <- 100
studyhours <- rnorm(n = N, mean = 6, sd = 1)
scores <- rnorm(n = N, mean = 12 * studyhours, sd = sqrt(1 - 0.3^2))
happiness <- rnorm(n = N, mean = scores-1.5* studyhours, sd = 10)
df <- data.frame(studyhours, scores, happiness)
pairs(df)

df1 <- lm(scores~happiness, df)
precis(df1)

##                   mean         sd       5.5%      94.5%
## (Intercept) 34.8112181 3.68185134 28.9269085 40.6955276
## happiness    0.5795067 0.05699576  0.4884164  0.6705969

df2 <- lm(scores ~ happiness + studyhours, df)
precis(df2)

##                    mean          sd         5.5%       94.5%
## (Intercept) -0.53799276 0.606379772 -1.507104749  0.43111923
## happiness    0.01313228 0.009548744 -0.002128458  0.02839302
## studyhours  11.96318525 0.143379934 11.734036419 12.19233407

#Happiness predicts scores when the study hours are already known but the study hours do not predict scores when happiness is already known, which means the association  between study hours and scores is spurious.

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.

#An example of a masked relationship involving the prediction of happiness ratings from the amount of alcohol one drinks and the amount of time one spends feeling ill.
N <- 100
rho <- 0.6
alcohol <- rnorm(n = N, mean = 0, sd = 1)
illness <- rnorm(n = N, mean = rho * alcohol, sd = sqrt(1 - rho^2))
happiness <- rnorm(n = N, mean = alcohol - illness, sd = 1)
d <- data.frame(happiness, alcohol, illness)
pairs(d)

# Masking relationship, the bivariate associations should be weak (widely variable).
m <- map(
  alist(
    happiness ~ dnorm(mu, sigma),
    mu <- a + ba * alcohol,
    a ~ dnorm(0, 5),
    ba ~ dnorm(0, 5),
    sigma ~ dunif(0, 5)
  ),
  data = d
)
precis(m)

##              mean         sd       5.5%      94.5%
## a     -0.09355758 0.10622002 -0.2633177 0.07620253
## ba     0.70087143 0.10691338  0.5300032 0.87173966
## sigma  1.06197277 0.07509043  0.9419638 1.18198178

m <- map(
  alist(
    happiness ~ dnorm(mu, sigma),
    mu <- a + bi * illness,
    a ~ dnorm(0, 5),
    bi ~ dnorm(0, 5),
    sigma ~ dunif(0, 5)
  ),
  data = d
)
precis(m)

##              mean         sd       5.5%       94.5%
## a     -0.08470341 0.12418503 -0.2831751  0.11376826
## bi    -0.28019914 0.12915566 -0.4866148 -0.07378345
## sigma  1.24104260 0.08775518  1.1007929  1.38129232

#In the multivariate regression, the masking relationship should be resolved.
m <- map(
  alist(
    happiness ~ dnorm(mu, sigma),
    mu <- a + ba * alcohol + bi * illness,
    a ~ dnorm(0, 5),
    ba ~ dnorm(0, 5),
    bi ~ dnorm(0, 5),
    sigma ~ dunif(0, 5)
  ),
  data = d
)
precis(m)

##             mean         sd       5.5%        94.5%
## a     -0.1312379 0.08699578 -0.2702739  0.007798201
## ba     0.9803105 0.09597154  0.8269295  1.133691572
## bi    -0.6989993 0.09921828 -0.8575693 -0.540429327
## sigma  0.8680636 0.06138105  0.7699648  0.966162384

# slopes for alcohol and illness became much larger in the multivariate model as the alcohol increases happiness while feelings of illness decrease happiness, 
#The bi-variate relationships of alcohol and feelings of illness to happiness tends to be masked.

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.

data(WaffleDivorce)
d <- WaffleDivorce
d$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)
d$logLDS <- log(d$LDS)
d$logLDS.s <- (d$logLDS - mean(d$logLDS)) / sd(d$logLDS)
hist(d$LDS)

hist(d$logLDS)

hist(d$logLDS.s)

#we can see that LDS has a high negative correlation with the divorce rate.

m <- 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 = d
)
precis(m)

##              mean         sd        5.5%      94.5%
## a     35.43510232 6.77534695 24.60678931 46.2634153
## bm     0.05342887 0.08261873 -0.07861182  0.1854696
## ba    -1.02956075 0.22469602 -1.38866839 -0.6704531
## bl    -0.60806453 0.29057096 -1.07245305 -0.1436760
## sigma  1.37872785 0.13838975  1.15755431  1.5999014

#Results show that slopes for median age at marriage vs percentage of LDS population are negative shows that are negatively related to divorce rates.
# Also the marriage rate isn’t a good predictor for divorce rate and states with an older median age of marriage had lower divorce rates.

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)

maddog = dagitty("dag{M -> A -> D}")
impliedConditionalIndependencies(maddog)

## D _||_ M | A

equivalentDAGs(maddog)

## [[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
## }

#We see that DAG implies the divorce rate (D) is independent of marriage rate on a conditional on the median age of marriage.

Assignment #4

Nishitha Vinnakota_ANLY505

2022-02-28

Chapter 5 - Many Variables and Spurious Waffles

Questions