Chapter 8 - Conditional Manatees

This chapter introduced interactions, which allow for the association between a predictor and an outcome to depend upon the value of another predictor. While you can’t see them in a DAG, interactions can be important for making accurate inferences. Interactions can be difficult to interpret, and so the chapter also introduced triptych plots that help in visualizing the effect of an interaction. No new coding skills were introduced, but the statistical models considered were among the most complicated so far in the book.

Place each answer inside the code chunk (grey box). The code chunks should contain a text response or a code that completes/answers the question or activity requested. Problems are labeled Easy (E), Medium (M), and Hard(H).

Finally, upon completion, name your final output .html file as: YourName_ANLY505-Year-Semester.html and publish the assignment to your R Pubs account and submit the link to Canvas. Each question is worth 5 points.

Questions

8E1. For each of the causal relationships below, name a hypothetical third variable that would lead to an interaction effect:

  1. Bread dough rises because of yeast.
  2. Education leads to higher income.
  3. Gasoline makes a car go.
#1. bread dough rises because of yeast, we can add some sugar to promote it. there is an interaction between yeast and sugar.
#2. education leads to higher income. it also depends on the major/fields. 
#3. gasoline makes a car go.Type of gasoline used/overall car condition can also affect.

8E2. Which of the following explanations invokes an interaction?

  1. Caramelizing onions requires cooking over low heat and making sure the onions do not dry out.
  2. A car will go faster when it has more cylinders or when it has a better fuel injector.
  3. Most people acquire their political beliefs from their parents, unless they get them instead from their friends.
  4. Intelligent animal species tend to be either highly social or have manipulative appendages (hands, tentacles, etc.).
#only 1

#1. cooking with low heat leads to caramelizing, conditional on the onions not drying out.

#2.  a car with many cylinders can go fast, whether or not it also has a good fuel injector.

#3. It just implies that one may be influenced by either parents or friends.

#4. the "or" shows they are not interaction.

8E3. For each of the explanations in 8E2, write a linear model that expresses the stated relationship.

\[\begin{align} “extent - caramelized”, μi = α + β_{H}H_{i} + β_{M}M_{i} + β_{HM}H_{i}M_{i}\\ \end{align}\]

\[\begin{align} “maximum-speed”, μi = α + β_{C}C_{i}+ β_{F}F_{i} \end{align}\]

\[\begin{align} “extent-conservative”, μi = α + β_{P}P_{i} + β_{F}F_{i} \end{align}\]

\[\begin{align} “intelligence” μi = α + β_{S}S_{i} + β_{M}M_{i} \end{align}\]

8M1. Recall the tulips example from the chapter. Suppose another set of treatments adjusted the temperature in the greenhouse over two levels: cold and hot. The data in the chapter were collected at the cold temperature. You find none of the plants grown under the hot temperature developed any blooms at all, regardless of the water and shade levels. Can you explain this result in terms of interactions between water, shade, and temperature?

# generally tulips are winter flower. high temperatures can prevent the all blooms. it's not a linear interaction

#there is a three way interaction here: the influence of water and shade depend upon each other, and both and their interaction depend upon temperature.

8M2. Can you invent a regression equation that would make the bloom size zero, whenever the temperature is hot?

we define Li is the ordinary linear model from the chapter. Then let Hi be a 0/1 indicator of whether or not the temperature was hot. Then:

\[\begin{align} μi = L_{i}(1-H_{i}) \end{align}\]

When Hi = 1, the entire model above is zero

8M3. In parts of North America, ravens depend upon wolves for their food. This is because ravens are carnivorous but cannot usually kill or open carcasses of prey. Wolves however can and do kill and tear open animals, and they tolerate ravens co-feeding at their kills. This species relationship is generally described as a “species interaction.” Can you invent a hypothetical set of data on raven population size in which this relationship would manifest as a statistical interaction? Do you think the biological interaction could be linear? Why or why not?

#the relationship between ravens and wolves is wolves do not need ravens, but ravens very much do benefit from wolves 

#Region Wolves Ravens
#1        13     44
#2        16     47
#3         8     29
#4        32     98
#5        16     61
#6        72     212

#this “interaction” is not a statical interaction effect. ravens depend on wolves means that we can partially predict raven density with wolf density. 

#A statistical interaction requires that some other third variable regulate the dependency of ravens on wolves. 
#for example: in some places, some small prey for ravens

8M4. Repeat the tulips analysis, but this time use priors that constrain the effect of water to be positive and the effect of shade to be negative. Use prior predictive simulation. What do these prior assumptions mean for the interaction prior, if anything?

we can us normal and exponential distributions. they are all positive, but can let the shade effect to be negative

prior information shows more water and more effect of light/ more water and less impact of shade. they are negative. we can use normal and exponential and let them to be minus.

these are the code with normal, mean of normal is exp(μ+σ2/2) we can load the data

library(rethinking) 
## Loading required package: rstan
## Warning: package 'rstan' was built under R version 4.0.3
## Loading required package: StanHeaders
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 4.0.3
## rstan (Version 2.21.2, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## Do not specify '-march=native' in 'LOCAL_CPPFLAGS' or a Makevars file
## Loading required package: parallel
## Loading required package: dagitty
## Warning: package 'dagitty' was built under R version 4.0.3
## rethinking (Version 2.01)
## 
## Attaching package: 'rethinking'
## The following object is masked from 'package:stats':
## 
##     rstudent
data(tulips)
d <- tulips
d$blooms_std <- d$blooms / max(d$blooms)
d$water_cent <- d$water - mean(d$water)
d$shade_cent <- d$shade - mean(d$shade)
m8M4a <- quap(
alist(
blooms_std ~ dnorm( mu , sigma ) ,
mu <- a + bw*water_cent - bs*shade_cent - bws*water_cent*shade_cent ,
a ~ dnorm( 0.5 , 0.25 ) ,
bw ~ dlnorm( 0 , 0.25 ) ,
bs ~ dlnorm( 0 , 0.25 ) ,
bws ~ dlnorm( 0 , 0.25 ) ,
sigma ~ dexp( 1 )
) , data=d )

we can do prior simulation to check if these priors make sense.

p <- extract.prior( m8M4a )
par(mfrow=c(1,3),cex=1.1) # 3 plots in 1 row
for ( s in -1:1 ) {
  idx <- which( d$shade_cent==s )
plot( d$water_cent[idx] , d$blooms_std[idx] , xlim=c(-1,1) , ylim=c(0,1) ,
xlab="water" , ylab="blooms" , pch=16 , col=rangi2 )
mtext( concat( "shade = " , s ) )
mu <- link( m8M4a , post=p , data=data.frame( shade_cent=s , water_cent=-1:1 ) )
for ( i in 1:20 ) lines( -1:1 , mu[i,] , col=col.alpha("black",0.3) )
}

it’s not good. we hope pior associations within possibility. we make the prior means closer to zero.

m8M4b <- quap(
alist(
blooms_std ~ dnorm( mu , sigma ) ,
mu <- a + bw*water_cent - bs*shade_cent - bws*water_cent*shade_cent ,
a ~ dnorm( 0.5 , 0.25 ) ,
bw ~ dlnorm( -2 , 0.25 ) ,
bs ~ dlnorm( -2 , 0.25 ) ,
bws ~ dlnorm( -2 , 0.25 ) ,
sigma ~ dexp( 1 )
) , data=d )

plot the code

p <- extract.prior( m8M4b )
par(mfrow=c(1,3),cex=1.1) # 3 plots in 1 row
for ( s in -1:1 ) {
  idx <- which( d$shade_cent==s )
plot( d$water_cent[idx] , d$blooms_std[idx] , xlim=c(-1,1) , ylim=c(0,1) ,
xlab="water" , ylab="blooms" , pch=16 , col=rangi2 )
mtext( concat( "shade = " , s ) )
mu <- link( m8M4b , post=p , data=data.frame( shade_cent=s , water_cent=-1:1 ) )
for ( i in 1:20 ) lines( -1:1 , mu[i,] , col=col.alpha("black",0.3) )
}

it’s better. it’s within the possible range of the outcome. there is still some uncertainty for the data.

8H1. Return to the data(tulips) example in the chapter. Now include the bed variable as a predictor in the interaction model. Don’t interact bed with the other predictors; just include it as a main effect. Note that bed is categorical. So to use it properly, you will need to either construct dummy variables or rather an index variable, as explained in Chapter 5.

load the data to check the bed variable

library(rethinking)
data(tulips)
d <- tulips
d$bed
##  [1] a a a a a a a a a b b b b b b b b b c c c c c c c c c
## Levels: a b c

there are three levels factor. we can code two or one index variable. dummy variable approach: need two dummy variables, it’s one less than dht enumber categories.

d$bed_b <- ifelse( d$bed=="b" , 1 , 0 ) 
d$bed_c <- ifelse( d$bed=="c" , 1 , 0 )

bed a is into the intercept, and each dummy variable’s coefficients will include the contrasts with bed a. the model using these dummy variables.

d$blooms_std <- d$blooms / max(d$blooms) 
d$water_cent <- d$water - mean(d$water)
d$shade_cent <- d$shade - mean(d$shade)
m1 <- quap(
alist(
blooms_std ~ dnorm( mu , sigma ),
mu <- a + bw*water_cent + bs*shade_cent +
bws*water_cent*shade_cent +
b_bed_b*bed_b + b_bed_c*bed_c ,
a ~ dnorm( 0 , 0.25 ),
c(bw,bs,bws,b_bed_b,b_bed_c) ~ dnorm(0,0.25),
sigma ~ dexp( 1 )
) , data=d )
precis(m1)
##               mean         sd        5.5%       94.5%
## a        0.2683617 0.03499307  0.21243598  0.32428733
## bw       0.2074364 0.02536895  0.16689189  0.24798085
## bs      -0.1138455 0.02536433 -0.15438264 -0.07330845
## bws     -0.1438934 0.03098856 -0.19341907 -0.09436765
## b_bed_b  0.1233159 0.04948615  0.04422749  0.20240435
## b_bed_c  0.1360289 0.04948707  0.05693902  0.21511880
## sigma    0.1081613 0.01468409  0.08469323  0.13162927

it looks that beds b and c are better than bed a. do the model using an index variable . estimate the intercept for every category, construct the index variable, use the coerce index function.

( d$bed_idx <- coerce_index( d$bed ) )
##  [1] 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3

bed a has index 1, bed b– index 2, bed c – index 3

m2 <- quap(
alist(
blooms_std ~ dnorm( mu , sigma ),
mu <- a[bed_idx] + bw*water_cent + bs*shade_cent +
bws*water_cent*shade_cent ,
a[bed_idx] ~ dnorm( 0 , 0.25 ),
c(bw,bs,bws) ~ dnorm( 0 , 0.25 ),
sigma ~ dexp( 1 )
) , data=d )
precis(m2,depth=2)
##             mean         sd        5.5%       94.5%
## a[1]   0.2630551 0.03578884  0.20585759  0.32025256
## a[2]   0.3861773 0.03582361  0.32892427  0.44343035
## a[3]   0.3988885 0.03582795  0.34162847  0.45614844
## bw     0.2074297 0.02542229  0.16679996  0.24805942
## bs    -0.1138431 0.02541761 -0.15446531 -0.07322082
## bws   -0.1438818 0.03105338 -0.19351114 -0.09425253
## sigma  0.1083907 0.01477657  0.08477494  0.13200656

beds b and c are better than bed a. and bed c a little better than bed b we can calculate the posterior distribution of the difference.

post <- extract.samples(m2)
diff_b_c <- post$a[,2] - post$a[,3]
PI( diff_b_c )
##          5%         94% 
## -0.09261304  0.06838092

we can see the expected difference. the posterior distribution of the difference has many probability on both sides of zero.

how to make all these? even though there were differences between the beds, can’t change the qualitative inference about the experiment in the bed analysis.