In this activity, we will be analysing data from an experiment in *Drosophila serrata* conducted by Jacqueline. She set up a paternal half-sibling breeding design to estimate genetic variance in male cuticular hydrocarbon (CHC) pheromone profile, and conducted a binomal fitness assay using these same individuals to get a measure of their mating success. *D. serrata* has 8 well-studied CHCs that are known to be a target of sexual selection via female choice, and to have an important ecological role in preventing desiccation. We have measured these traits of interest in all offspring from the breeding design, and we also assayed the mating success of these individuals. Our goal is to better understand the genetic basis of these traits and whether they are under selection. In this activity we will analyse a subset of 4 CHC traits in order to estimate the genetic variance and covariance between them (ie. the **G** matrix), and to estimate direction and non-linear selection acting on these traits.

The original data were published in *Sztepanacz, J.L. and Rundle, H.D., 2012. Reduced genetic variance among high fitness individuals: inferring stabilizing selection on male sexual displays in Drosophila serrata. Evolution, 66(10), pp.3101-3110* and the full data set is available on Dryad.

If you have not done so already, open this .Rmd file in RStudio to take advantage of the RMarkdown formatting.

First, set up your R environment by setting your working directory and loading packages.

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(cache = TRUE, cache.lazy = FALSE)
knitr::opts_knit$set(root.dir = '/Users/jsztepanacz/Library/CloudStorage/OneDrive-UniversityofToronto/EQGW workshop mountain lake 2025/Lectures and Activities/Activity idea') #change this path to match your setup

library(MCMCglmm)
library(dplyr)
library(ggplot2)
library(stringr)
library(evolvability)
library(pedtricks)
library(celltrackR)

```

First, load the data. We have one file of data, contained in `Data.csv`, and one file with the pedigree, contained in `fullped.csv`. Load these data into R.

```{r load_data}
data <- read.csv("Data.csv", header = T)
ped <- read.csv("fullped.csv", header = T)
```

Before using the data, inspect the structure of the data in order to ensure that the classes of all variables are correct. Use this code block to convert variables to factors as needed. The data contains the following variables:

-   Sire: ID of the father
-   Dam: ID of the mother
-   animal: ID of offspring
-   Column: gas chromatography column the sample was run on
-   Block: temporal block of the experiment
-   Chosen/Rejected: mating success of the individual; 1=chosen, 0=rejected
-   Vial: same as Dam,
-   l2: log-contrast value of CHC2
-   l3: log-contrast value of CHC3
-   l4: log-contrast value of CHC4
-   l5: log-contrast value of CHC5

```{r factors}
head(data)
head(ped)

data$Sire<-as.factor(data$Sire)
data$Dam<-as.factor(data$Dam)
data$animal<-as.factor(data$animal)
data$Column<-as.factor(data$Column)
data$Block<-as.factor(data$Block)
data$Chosen.Rejected<-as.numeric(data$Chosen.Rejected)
data$Vial<-as.factor(data$Vial)

for (x in 1:3) ped[, x] <- as.factor(ped[, x])


#this function manipulates the pedigree to conform to the structure required for MCMCglmm
ped<-fix_ped(ped=ped)

```

# Explore data

Before jumping in to analyses, take a moment to explore the data. Produce graphs of your choice to visualize patterns in the distributions of the variables and the relationships between them.

```{r visualize_data_initially}
plot(data[,8:11])
hist(data$l2)
hist(data$l3)
hist(data$l4)
hist(data$l5)

#summarise the data
summary(data[,8:11])

# standardise the data to mean 0 variance 1

data_std<-scale(data[,8:11])
colnames(data_std)<-c("l2s", "l3s", "l4s", "l5s")
data<-cbind(data, data_std)

cor(data[,12:15])

```

**Are there any outliers in the data? How are the data distributed? Do they meet the assumptions of the model you would like to fit? Is there anything interesting to note in the summary statistics? Is there anything else we should be looking for when exploring the data**

**Do you predict that there is a genetic correlation between any of the traits? If so, why?**

If there are any outliers in the data remove them here, and then standardise the data to mean of 0 and variance of 1

```{r remove_outliers}
pairs(data[,8:11]) 

#plot l3 vs l5 with animals labelled to figure out which to remove
ggplot(data=data, aes(x = l3, y=l5))+geom_text(aes(label=animal))

#remove outliers
remove<-c("2143", "1093", "2153", "2154")
data <- data %>% filter(!animal %in% remove) 

# standardise the data to mean 0 variance 1

data_std<-scale(data[,8:11])
colnames(data_std)<-c("l2s", "l3s", "l4s", "l5s")
data<-cbind(data, data_std)

cor(data[,12:15])

```

#Explore pedigree; learn how to plot a pedigree, 

```{r draw_pedigree_matrix}
#here we are visualizing the pedigree for our interest
#this code draws the pedigree structure; unfortunately for the half-sib full-sib design, it doesn't look very interesting
draw_ped(ped)
#this code draws the additive genetic relatedness matrix-- again not very interesting
draw_pedA(ped)
```
So you can see what a different kind of pedigree might look like, let's draw an example

```{r draw_example_pedigree}
#here is an example of a pedigree that does look interesting
data(gryphons)
test.pedigree <- fix_ped(gryphons[,1:3])
## draw the gryphon pedigree by pedigree depth
draw_ped(test.pedigree)
draw_pedA(test.pedigree)
```

# Estimate G

Next, build a model to investigate the genetic variance and covariance between the log-contrast CHC traits. Even though we only have 4 traits, we still need to estimate quite a few parameters, in fact,

*n\**(*n*+1)/2 parameters, where *n* is the number of traits. So, 10 in this case.

Use this code block to generate the inverse A matrix that you can provide to `MCMCglmm()`.

```{r inv_A_mat}
Ainv <- inverseA(ped)$Ainv
```

Then, use this code block to specify your prior and run your model using `MCMCglmm()`. In this model, include `Block` and `Column` as fixed effects, and vial as a random effect. Choose a small number of iterations to run in the interest of time. Check out this tutorial https://devillemereuil.legtux.org/wp-content/uploads/2012/12/tuto_en.pdf for more information on how ti set up a a multivariate analysis.

```{r first_model}

#specify a prior for a 4x4 covariance matrix
prior <- list(
  G = list(G1 = list(V = diag(4)/3, nu = 1.002),
           G2 = list(V = diag(4)/3, nu = 1.002)),
  R = list(V = diag(4)/3, nu = 1.002)
)

#build a model
model <- MCMCglmm(cbind(l2s, l3s, l4s, l5s) ~ trait-1+Block+Column,
  random = ~ us(trait):animal + us(trait):Vial,
  rcov = ~ us(trait):units,
  family = c("gaussian", "gaussian", "gaussian", "gaussian"),
  ginverse = list(animal = Ainv),
  data = data, prior = prior, verbose = T,
  #nitt=, burnin=, thin=,
)

#summarize the model
summary(model)

```
After running your model, assess the MCMC run to see whether it ran for a sufficient length of time. What you are looking for is a trace plot that looks like a fuzzy caterpillar, and a density plot that looks like a normal distribution

```{r assess_model_b}
#plot the trace and posterior density for a variable
#plot the trace and posterior density for the `traittarsus:traittarsus.animal` variable
plot(model$Sol)
plot(model$VCV[,1:10])
plot(model$VCV[,11:20])
plot(model$VCV[,21:32])
#calculate autocorrelation between estimates
autocorr.diag(model$VCV)

```

The following bit of code makes a nice table with the estimates and effective sample sizes beside them.The effective sample size indicates the number of independent samples from you MCMC chain, which you want to be as large as possible. You are looking for effective sample sizes close to the actual sample sizes from your MCMC chain, smaller numbers indicates autocorrelation in the estimates.

```{r model_summary}
    clean.MCMC <- function(x) {
    sols <- summary(x)$solutions  ## pull out relevant info from model summary
    Gcovs <- summary(x)$Gcovariances
    Rcovs <- summary(x)$Rcovariances
    fixed <- data.frame(row.names(sols), sols, row.names = NULL)  ## convert to dataframes with the row.names as the first col
    random <- data.frame(row.names(Gcovs), Gcovs, row.names = NULL)
    residual <- data.frame(row.names(Rcovs), Rcovs, row.names = NULL)
    names(fixed)[names(fixed) == "row.names.sols."] <- "variable"  ## change the columns names to variable, so they all match
    names(random)[names(random) == "row.names.Gcovs."] <- "variable"
    names(residual)[names(residual) == "row.names.Rcovs."] <- "variable"
    fixed$effect <- "fixed"  ## add ID column for type of effect (fixed, random, residual)
    random$effect <- "random"
    residual$effect <- "residual"
    modelTerms <- as.data.frame(bind_rows(fixed, random, residual))  # merge it all together
    }
    getName.MCMC <- function(x) {deparse(substitute(x))}
    
oneModel <- clean.MCMC(model)  # get all the info from summary(modelName)
oneModel$modelName <- getName.MCMC(model)  # add the model's name in a new column
oneModel 
```


**Do you think that your model ran for a sufficient number of MCMC iterations? Describe what features you used to reach your conclusion.**



If you model did not run for a sufficient amount of time, import the output from an analysis that I ran previously and continue your analysis using this output.

```{r load_mcmc_output}
#load(file = "./activity/MCMCglmm_model_z_Longrun.rda")
```

Now, let's look at the estimates for the variance components in the model.

```{r model_variance, results = "hold"}
writeLines("Variance components:")
#calculate the modes of the posterior distribution for the variance components
posterior.mode(model$VCV)
```

Estimate the heritability of each trait

```{r heritability_estimate, results = "hold"}
#estimate heritability of l2-l5
heritability.l2 <- model$VCV[, "traitl2s:traitl2s.animal"] / (model$VCV[, "traitl2s:traitl2s.animal"] +model$VCV[, "traitl2s:traitl2s.Vial"]+ model$VCV[, "traitl2s:traitl2s.units"])
heritability.l3 <- model$VCV[, "traitl3s:traitl3s.animal"] / (model$VCV[, "traitl3s:traitl3s.animal"] +model$VCV[, "traitl3s:traitl3s.Vial"]+ model$VCV[, "traitl3s:traitl3s.units"])
heritability.l4 <- model$VCV[, "traitl4s:traitl4s.animal"] / (model$VCV[, "traitl4s:traitl4s.animal"] + model$VCV[, "traitl4s:traitl4s.Vial"]+  model$VCV[, "traitl4s:traitl4s.units"])
heritability.l5 <- model$VCV[, "traitl5s:traitl5s.animal"] / (model$VCV[, "traitl5s:traitl5s.animal"] +model$VCV[, "traitl5s:traitl5s.Vial"]+ model$VCV[, "traitl5s:traitl5s.units"])


#calculate the mode of the posterior distributions for heritability and the HPD intervals
h2.l2<-c(posterior.mode(heritability.l2), HPDinterval(heritability.l2))

h2.l3<-c(posterior.mode(heritability.l3), HPDinterval(heritability.l3))

h2.l4<-c(posterior.mode(heritability.l4), HPDinterval(heritability.l4))

h2.l5<-c(posterior.mode(heritability.l5), HPDinterval(heritability.l5))

#arrange in a data frame and produce a table
heritabilities<-rbind(h2.l2, h2.l3, h2.l4, h2.l5)
colnames(heritabilities)<-c("h2", "lowerHPD", "upperHPD")
knitr::kable(heritabilities[,1:3], "simple", digits=2)


```

**What is your estimate of the heritability for each trait? Include the 95% HPD.** 


Let's arrange the estimates from our model into a G-matrix

```{r genetic_covariance_matrix}

#put the posterior mode of the VCV elements in a matrix

pointG<-matrix(posterior.mode(model$VCV[,1:16]),4,4)
pointG

#turn into correlation matrix

pointGcor<-cov2cor(pointG)
pointGcor

```

**What is the genetic covariance between each pair of traits? Include the 95% HPD**.

**What are the genetic correlations between all traits?**

**What would you conclude about the genetic architecture of these traits?**


## Quantify potential genetic constraints

Now that we have estimated **G**, we would like to determine how genetic variance is distributed across multivariate space.

Start by doing an eigenanalysis of **G** and produce a scree plot of the eigenvalues.

```{r eigenanalysis of G}

eigen(pointG)

#plot the eigenvalues
plot(eigen(pointG)$values)

#plot the eigenvalues with uncertainty

#save an object of the MCMC samples of G
Gmcmc<-model$VCV[,1:16]

#eigenanalysis of each posterior sample
eigs<-array(dim=c(1000,4))
for (i in 1:1000){
 eigs[i,]<-eigen(matrix(Gmcmc[i,],4,4))$values
}

eig1l<-quantile(eigs[,1], 0.025)
eig1u<-quantile(eigs[,1], 0.975)

eig2l<-quantile(eigs[,2], 0.025)
eig2u<-quantile(eigs[,2], 0.975)

eig3l<-quantile(eigs[,3], 0.025)
eig3u<-quantile(eigs[,3], 0.975)

eig4l<-quantile(eigs[,4], 0.025)
eig4u<-quantile(eigs[,4], 0.975)

#plot

eigen_plot<-data.frame(seq(1:4),eigen(pointG)$values, rbind(eig1l, eig2l, eig3l, eig4l), rbind(eig1u, eig2u, eig3u, eig4u))
colnames(eigen_plot)<-c("eigenvalue", "value", "lci", "uci")

p<-ggplot(data = eigen_plot, aes(x = eigenvalue, y = value)) +
  geom_point(cex=5) + 
  geom_errorbar(aes(x = eigenvalue, ymin = lci, ymax = uci))+
   scale_y_continuous(expand = c(0, 0), limits = c(-1, 2)) 

p


```

**Are there any multivariate trait combinations that have low genetic variance? How many?**

**What is the effective dimensionality of the matrix?**

## Fit a null-model

We would like to know if our estimates of genetic variance and covariance are different from sampling error (ie. are they statistically significant). MCMCglmm constrains estimates of **G** to be positive-definite, which means that estimates of variances, including eigenvalues, must be greater than or equal to 0. Therefore, the lower bound of an HPD cannot be less than 0, which means we can't compare it to 0 to determine statistical support.

We need to fit a null-model to compare our estimates to. There are different ways one could go about doing this, but the most straightforward is to randomise the phenotypes across the pedigree to keep the same data structure, but remove any genetic signal, and then re-fit the model. There are technical questions about the best way to do this randomisation in the presence of fixed effects, but let's ignore those questions for now.

```{r randomise_phenotypes}
#randomise the phenotypes to break any genetic links in the data

shuffled_data<-data[sample(1:nrow(data)), 12:15]
colnames(shuffled_data)<-c("l2sr", "l3sr", "l4sr", "l5sr")
data<-cbind(data, shuffled_data)

#fit the same model that you used to estimate G

#build a model
model.r <- MCMCglmm(cbind(l2sr, l3sr, l4sr, l5sr) ~ trait - 1+Block+Column,
  random = ~ us(trait):animal + us(trait):Vial,
  rcov = ~ us(trait):units,
  family = c("gaussian", "gaussian", "gaussian", "gaussian"),
  pedigree= ped,
  data = data, prior = prior, verbose = T
)

#summarize the model
summary(model.r)

#this will take a while to run, so load the output from a model that I have run for you

```

Produce a scree plot of the eigenvalues from your estimated G and your null G with 95% HPD on the same figure.

```{r eigenanalysis_of_G}

#put the posterior mode of the VCV elements in a matrix

pointG.r<-matrix(posterior.mode(model.r$VCV[,1:16]),4,4)
pointG.r

#turn into correlation matrix

pointG.rcor<-cov2cor(pointG.r)
pointG.rcor

#save an object of the MCMC samples of G
Gmcmc.r<-model.r$VCV[,1:16]

#eigenanalysis of each posterior sample
eigs.r<-array(dim=c(1000,4))
for (i in 1:1000){
 eigs.r[i,]<-eigen(matrix(Gmcmc.r[i,],4,4))$values
}

eig1l.r<-quantile(eigs.r[,1], 0.025)
eig1u.r<-quantile(eigs.r[,1], 0.975)

eig2l.r<-quantile(eigs.r[,2], 0.025)
eig2u.r<-quantile(eigs.r[,2], 0.975)

eig3l.r<-quantile(eigs.r[,3], 0.025)
eig3u.r<-quantile(eigs.r[,3], 0.975)

eig4l.r<-quantile(eigs.r[,4], 0.025)
eig4u.r<-quantile(eigs.r[,4], 0.975)

#plot the eigenvalues


eigen_plot.r<-data.frame(seq(1:4),eigen(pointG.r)$values, rbind(eig1l.r, eig2l.r, eig3l.r, eig4l.r), rbind(eig1u.r, eig2u.r, eig3u.r, eig4u.r))
colnames(eigen_plot.r)<-c("eigenvalue", "value", "lci", "uci")

p + geom_point(data = eigen_plot.r, aes(x = eigenvalue, y = value), color = "purple", cex=5) + 
  geom_errorbar(data = eigen_plot.r, aes(x = eigenvalue, ymin = lci, ymax = uci))


```

**Do any of your previous interpretations change now that you have compared your estimates to the null model?**


# Estimate selection

Next we will estimate how selection acts on these 4 CHC traits. The component of fitness that we will use is mating success. The standard Lande-Arnold approach to estimate selection regresses trait values on relative fitness.

First calculate the relative fitness of individuals in the data

```{r relative_fitness}
#calculate relative fitness and add it to the data
str(data)

data$relfit<-data$Chosen.Rejected/mean(data$Chosen.Rejected)

```

## Directional selection

Fit the regression model using the lm function to estimate the selection differential for each trait. Don't forget to include block and column as fixed effects

```{r selection_differential}
#regress each trait on relative fitness separately and save the estimate and standard error

sel.dif<-rbind(
summary(lm(relfit~Block+Column+l2s, data=data))$coefficients[4,1:2],
summary(lm(relfit~Block+Column+l3s, data=data))$coefficients[4,1:2],
summary(lm(relfit~Block+Column+l4s, data=data))$coefficients[4,1:2],
summary(lm(relfit~Block+Column+l5s, data=data))$coefficients[4,1:2])

```

We fit a linear regression to estimate selection, even though the response variable was not normally distributed. That is totally fine, and these are the correct estimates of the gradients! However, to determine which of these coefficients are statistically significant, we should fit a binomial model with a logit link function to get the appropriate p-values.

```{r statistical_significance}
#fit a model using the glm function
sel.dif.p<-rbind(
summary(glm(Chosen.Rejected ~Block+Column+l2s, family=binomial(link=logit), data=data))$coefficients[4,4],
summary(glm(Chosen.Rejected ~Block+Column+l3s, family=binomial(link=logit), data=data))$coefficients[4,4],
summary(glm(Chosen.Rejected ~Block+Column+l4s, family=binomial(link=logit), data=data))$coefficients[4,4],
summary(glm(Chosen.Rejected ~Block+Column+l5s, family=binomial(link=logit), data=data))$coefficients[4,4])

sel.dif.table<-data.frame(cbind(sel.dif, sel.dif.p))
colnames(sel.dif.table)<-c("estimate", "se", "p-value")
knitr::kable(sel.dif.table, "simple", digits=2)

```

**Which trait is under the strongest directional selection? Weakest?**

Now let's use multiple regression to disentangle the effects of direct versus indirect selection. Use the lm function to fit a multiple regression of all 4 traits on relative fitness. Don't forget to include block and column as fixed effects.

```{r selection_gradient}
#fit multiple regression model
sel.grad<-lm(relfit~Block+Column+l2s+l3s+l4s+l5s, data=data)
beta<-sel.grad$coefficients[4:7]
print(beta)

#fit a model using glm to determine statistical significance for directional selection overall
a3<-glm(Chosen.Rejected ~Block+Column+l2s+l3s+l4s+l5s, family=binomial(link=logit), data=data)
a4<-glm(Chosen.Rejected ~Block+Column, family=binomial(link=logit), data=data)
anova(a3,a4)


#make a table of the selection gradients and p-values for individual coefficients
sel.grad.table<-cbind(beta, summary(a3)$coefficients[4:7,4])
colnames(sel.grad.table)<-c("selection gradient", "p-value")
knitr::kable(sel.grad.table, "simple", digits=2)

```

**Now which trait is under the strongest directional selection? Weakest?**

**What is the total strength of directional selection on CHCs? hint: calculate the length of the vector**


## Quadratic selection

We might also be interested in whether there is stabilising or disruptive selection on each CHC trait. We will fit the same type of model as before, but now estimate the quadratic terms (ie. the square of each trait). Don't forget to include the linear terms as covariates in the model, or to double the quadratic coefficients! (see Stinchcombe J. R., A. F. Agrawal, P. A. Hohenlohe, S. J. Arnold, and M. W. Blows. 2008. Estimating nonlinear selection gradients using quadratic regression coefficients: double or nothing? Evolution 62:2435-2440).

```{r quadratic_selection}

a5<-lm(relfit~Block+Column + l2s+l3s+l4s+l5s+I(l2s^2)+I(l3s^2)+I(l4s^2)+I(l5s^2), data=data)
summary(lm(relfit~Block+Column + l2s+l3s+l4s+l5s+I(l2s^2)+I(l3s^2)+I(l4s^2)+I(l5s^2), data=data))

quad_coeffs<-a5$coefficients[8:11]*2
quad_coeffs

## compare the fit of a model with and without quadratic terms

a6<-glm(Chosen.Rejected~Block+Column + l2s+l3s+l4s+l5s+I(l2s^2)+I(l3s^2)+I(l4s^2)+I(l5s^2), family=binomial(link=logit), data=data)

quad.sel.table<-cbind(quad_coeffs, summary(a6)$coefficients[8:11,4])
colnames(quad.sel.table)<-c("quad.grad", "p-value")
knitr::kable(quad.sel.table, "simple", digits=2)


```

**Are any traits under stabilising selection?**

**Are there any traits under disruptive selection?**

## Quadratic and correlational selection

The previous analysis looks at individual traits in isolation, but there might be selection on combinations of traits (ie. correlational selection). To estimate quadratic and correlational selection we need to expand our model to include both quadratic and correlational terms. With 4 traits, that is already quite a few parameters! 4 linear + 10 quadratic and correlational terms = 14. The matrix of quadratic and correlational selection coefficients that we are going to estimate is called 'gamma'.

```{r correlational_selection}
a6<-lm(relfit~Block+Column + l2s+l3s+l4s+l5s+I(l2s^2)+I(l3s^2)+I(l4s^2)+I(l5s^2)+I(l2s*l3s)+I(l2s*l4s)+I(l2s*l5s)+I(l3s*l4s)+I(l3s*l5s)+I(l4s*l5s), data=data)
summary(a6)

# arrange the coefficients into the gamma matrix
gamma<-diag(a6$coefficients[8:11]*2)
gamma[2:4, 1]<-a6$coefficients[12:14]
gamma[3:4, 2]<-a6$coefficients[15:16]
gamma[4, 3]<-a6$coefficients[16]
gamma<-gamma+t(gamma)- diag(diag(gamma))

```

**What can you say about patterns of correlational selection on CHCs?**

## Canonical rotation of gamma

In the same way that it can be difficult to interpret **G** just by looking at it, it can also be difficult to interpret the multivariate pattern of selection just by looking at gamma. We can conduct an eigenanalysis of gamma the same way that we did with **G** to get a better understanding of multivariate selection. This is called the canonical rotation (see Blows for original paper, and x for a criticism of the approach).

```{r canonical_rotation}

eigen(gamma)

#plot the eigenvalues of gamma

plot(eigen(gamma)$values)
abline(h=0, lty=2)

```

**How many multivariate trait combinations are under disruptive selection?**

**How many multivariate trait combinations are under stabilising selection?**

**Overall, is there more disruptive or stabilising selection?**

# Combine G and selection

Note: In a real analysis you will need to account for the uncertainty of your estimates of **G** and selection. For this exercise, however, we will use the point estimates only. Use the posterior mode for **G**, and use the estimated selection vector from your linear model for beta.

## What is the predicted response to selection?

First we will calculate the predicted response to directional selection on the 4 CHC traits

```{r}
deltazbar<-pointG %*% beta

```

## Alignment of directional selection and genetic variance

To calculate the correlations between vectors they need to be standardised to unit length first. Check whether they are, and if not, standardise to unit length before answering the follwing questions.

```{r standardise_unit_length}

l2_norm<-sqrt(sum(beta^2))
l2_norm

beta_std<-beta/l2_norm
beta_std<-as.matrix(beta_std, 4,1)



gmax<-eigen(pointG)$vectors[,1]
#check whether it's unit length
sqrt(sum(gmax^2))
#check beta_std for unit length
sqrt(sum(beta_std^2))

```

**Is directional selection aligned with the major axis of genetic variance? hint: calculate the vector correlation between gmax and the directional selection gradient**

```{r}

#calculate the angle     
angle<-acos(cor(beta_std, gmax))*180/pi

if(angle<90){angle}else{180-angle}
```

**How much genetic variance is there in the direction of linear selection? hint: project the selection gradient through G**

```{r}
#project beta through G
#project beta through G

t(beta_std) %*% pointG %*% beta
```

**Is that a lot or not much? How did you come to your conclusion?**


# Congratulations! You have now completed a multivariate genetic and selection analysis
