The biological question that I am trying to address is: does observed heterozygosity change between years, and is it influenced by population and locus, and possibly $N_e$, but $N_e$ seems to be mucking up my model selection. My plan is to do this using linear models (lm
function in R), using hierarchical model selection and lsmeans.
I have heterozygosity at a given locus, for each of three populations for each of two years. So, the heterozygosity value at a locus is the average for a population for a given year. The $N_e$ values that I have are just point estimates. I have these for each population for each year (so six values in total)
The questions that I have for this forum are:
- Am I addressing problems with normality correctly (see below)? I will ask my questions surrounding this excluding $N_e$ from my models.
- Do you have any suggestions for why $N_e$ is causing errors (see below)? Suggestions for how to troubleshoot would be much appreciated.
A snippet of the dataframe that I read in:
Locus PopID Pop Year ObsHet HoLogit
13 6 Chalifour 2003 0.2857 -0.39797
13 2 Chalifour 2015 0.2439 -0.49137
13 7 IcPe 2003 0.425 -0.13128
13 5 IcPe 2015 0.3418 -0.28459
13 9 Takwa 2003 0.4516 -0.08434
13 3 Takwa 2015 0.2059 -0.58622
Here’s the structure of my data, as read in by R after making year a factor and the new HoLogitR column:
> str(r12gNoOddsGenets)
'data.frame': 51140 obs. of 7 variables:
$ Locus : int 13 13 13 13 13 13 51 51 51 51 ...
$ PopID : int 6 2 7 5 9 3 2 7 5 9 ...
$ Pop : Factor w/ 3 levels "Chalifour","IcPe",..: 1 1 2 2 3 3 1 2 2 3 ...
$ Year : Factor w/ 2 levels "2003","2015": 1 2 1 2 1 2 2 1 2 1 ...
$ ObsHet : num 0.286 0.244 0.425 0.342 0.452 ...
$ HoLogit : num -0.398 -0.4914 -0.1313 -0.2846 -0.0843 ...
$ HoLogitR: num -0.864 -1.063 -0.287 -0.62 -0.184 ...
ObsHet
is untransformed, HoLogit
is Ho transformed in excel and imported (transformation is =(H0/1-H0), and I had to change all zeros and 1s to NAs. I’ll discuss this below, as I am sure it is not correct to do this, but anyway, HoLogitR
is the logit transformation that I did in the car
package of R. I’m not sure that this is correct either, and I will discuss that below.
First,
Plots of untransformed ObsHet
These look very clearly right-skewed, not normally distributed
And the full model, with residual plots, using untransformed data.
Ho_comp = lm(ObsHet ~ Year + Pop + Locus + Year:Pop + Year:Locus + Pop:Locus + Year:Pop:Locus, data = r12gNoOddsGenets)
I am really not an expert with this stuff, but my reading of these plots is that homoscedasticity is basically okay, but normality is not.
Based on positive skewness evident in the histograms above, a square root transformation could be appropriate. With a square root transformation
Ho_comp = lm(sqrt(ObsHet) ~ Year + Pop + Locus + Year:Pop + Year:Locus + Pop:Locus + Year:Pop:Locus, data = r12gNoOddsGenets)
The QQ plot still looks pretty skewed. So, I tried the logit transformation.
Here are my first issues.
• As a cursory check, I did a logit transformation in excel (Ho/1-Ho
). But there is then no way to deal with 0s and 1s in the data. So, I changed them to NAs. This value corresponds to the HoLogit
column in the snippet of my dataframe below.
• Since getting rid of the zeros and 1s is an issue, I also tried to use the logit
function in the car
package in R. This corresponds to the HoLogitR
column. You can see that the values are very different. Presumably this has something to do with the “percent” and the “adjust” parameter, which I have set to the default 0.025. Heterozygosity is obviously bounded by 0 and 1. Here’s my code for the R transformation.
> r12gNoOddsGenets["HoLogitR"] <- logit(r12gNoOddsGenets$ObsHet, percents=max(r12gNoOddsGenets$ObsHet, na.rm = TRUE) > 1, 0.025)
> r12gNoOddsGenets[1:25,]
Locus PopID Pop Year ObsHet HoLogit HoLogitR
1 13 6 Chalifour 2003 0.2857 -0.39797041 -0.8644282
2 13 2 Chalifour 2015 0.2439 -0.49136744 -1.0631654
Thoughts for how I can deal address this issue, the zero and 1 thing, and how to deal with the percents and adjustment parameters?
Here are the plots using the HoLogit column, so the one that I did in excel, removing the zeros and ones.
Ho_compa = lm(HoLogit ~ Year + Pop + Locus + Year:Pop + Year:Locus+ Pop:Locus+ Year:Pop:Locus, data = r12gNoOddsGenets)
I’m really no expert at this stuff, but I think the residuals look best with these plots. QQ still looks skewed, but thoughts on this?
Now, addressing the $N_e$ issue. Regardless of the fit of the data, I’ll just show what happens with the untransformed data when I include $N_e$ in the models.
> str(r12gNoOddsGenets)
'data.frame': 51140 obs. of 8 variables:
$ Locus : int 13 13 13 13 13 13 51 51 51 51 ...
$ PopID : int 6 2 7 5 9 3 2 7 5 9 ...
$ Pop : Factor w/ 3 levels "Chalifour","IcPe",..: 1 1 2 2 3 3 1 2 2 3 ...
$ Year : Factor w/ 2 levels "2003","2015": 1 2 1 2 1 2 2 1 2 1 ...
$ Ne : int 2366 2760 2403 2463 1741 3146 2760 2403 2463 1741 ...
$ ObsHet : num 0.286 0.244 0.425 0.342 0.452 ..
> #Hobs hierarchical model selection
> Ho_comp = lm(ObsHet ~ Year + Pop + Locus + Ne + Year:Pop + Year:Locus + Year:Ne + Pop:Locus + Pop:Ne + Locus:Ne + Year:Pop:Locus:Ne, data = r12gNoOddsGenets)
> Ho_mod1 = lm(ObsHet ~ Year + Pop + Locus + Ne + Year:Pop + Year:Locus + Year:Ne + Pop:Locus + Pop:Ne + Locus:Ne, data = r12gNoOddsGenets)
> anova(Ho_comp, Ho_mod1)
Analysis of Variance Table
Model 1: ObsHet ~ Year + Pop + Locus + Ne + Year:Pop + Year:Locus + Year:Ne +
Pop:Locus + Pop:Ne + Locus:Ne + Year:Pop:Locus:Ne
Model 2: ObsHet ~ Year + Pop + Locus + Ne + Year:Pop + Year:Locus + Year:Ne +
Pop:Locus + Pop:Ne + Locus:Ne
Res.Df RSS Df Sum of Sq F Pr(>F)
1 51128 1530.3
2 51129 1530.3 -1 -0.0018427 0.0616 0.804
> Ho_mod2 = lm(ObsHet ~ Year + Pop + Locus + Ne + Year:Pop + Year:Locus + Year:Ne + Pop:Locus + Pop:Ne + Locus:Ne, data = r12gNoOddsGenets)
> Ho_mod3 = lm(ObsHet ~ Year + Pop + Locus + Ne + Year:Pop + Year:Locus + Year:Ne + Pop:Locus + Pop:Ne, data = r12gNoOddsGenets)
> anova(Ho_mod2, Ho_mod3)
Analysis of Variance Table
Model 1: ObsHet ~ Year + Pop + Locus + Ne + Year:Pop + Year:Locus + Year:Ne +
Pop:Locus + Pop:Ne + Locus:Ne
Model 2: ObsHet ~ Year + Pop + Locus + Ne + Year:Pop + Year:Locus + Year:Ne +
Pop:Locus + Pop:Ne
Res.Df RSS Df Sum of Sq F Pr(>F)
1 51129 1530.3
2 51130 1530.3 -1 -0.0022629 0.0756 0.7833
> Ho_mod4 = lm(ObsHet ~ Year + Pop + Locus + Ne + Year:Pop + Year:Locus + Year:Ne + Pop:Locus + Pop:Ne, data = r12gNoOddsGenets)
> Ho_mod5 = lm(ObsHet ~ Year + Pop + Locus + Ne + Year:Pop + Year:Locus + Year:Ne + Pop:Locus, data = r12gNoOddsGenets)
> anova(Ho_mod4, Ho_mod5)
Analysis of Variance Table
Model 1: ObsHet ~ Year + Pop + Locus + Ne + Year:Pop + Year:Locus + Year:Ne +
Pop:Locus + Pop:Ne
Model 2: ObsHet ~ Year + Pop + Locus + Ne + Year:Pop + Year:Locus + Year:Ne +
Pop:Locus
Res.Df RSS Df Sum of Sq F Pr(>F)
1 51130 1530.3
2 51130 1530.3 0 0
I get zeros for the F and p values whenever there is a Ne locus combination as the last variables shown in the list for the models. See below for a second time this happens. Any thoughts about how to trouble shoot this?
> Ho_mod6 = lm(ObsHet ~ Year + Pop + Locus + Ne + Year:Pop + Year:Locus + Year:Ne + Pop:Locus, data = r12gNoOddsGenets)
> Ho_mod7 = lm(ObsHet ~ Year + Pop + Locus + Ne + Year:Pop + Year:Locus + Year:Ne, data = r12gNoOddsGenets)
> anova(Ho_mod6, Ho_mod7)
Analysis of Variance Table
Model 1: ObsHet ~ Year + Pop + Locus + Ne + Year:Pop + Year:Locus + Year:Ne +
Pop:Locus
Model 2: ObsHet ~ Year + Pop + Locus + Ne + Year:Pop + Year:Locus + Year:Ne
Res.Df RSS Df Sum of Sq F Pr(>F)
1 51130 1530.3
2 51132 1530.5 -2 -0.15553 2.5982 0.07442 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> Ho_mod8 = lm(ObsHet ~ Year + Pop + Locus + Ne + Year:Pop + Year:Locus + Year:Ne, data = r12gNoOddsGenets)
> Ho_mod9 = lm(ObsHet ~ Year + Pop + Locus + Ne + Year:Pop + Year:Locus, data = r12gNoOddsGenets)
> anova(Ho_mod8, Ho_mod9)
Analysis of Variance Table
Model 1: ObsHet ~ Year + Pop + Locus + Ne + Year:Pop + Year:Locus + Year:Ne
Model 2: ObsHet ~ Year + Pop + Locus + Ne + Year:Pop + Year:Locus
Res.Df RSS Df Sum of Sq F Pr(>F)
1 51132 1530.5
2 51132 1530.5 0 0
> Ho_mod10 = lm(ObsHet ~ Year + Pop + Locus + Ne + Year:Pop + Year:Locus, data = r12gNoOddsGenets)
> Ho_mod11 = lm(ObsHet ~ Year + Pop + Locus + Ne + Year:Pop, data = r12gNoOddsGenets)
> anova(Ho_mod10, Ho_mod11)
Analysis of Variance Table
Model 1: ObsHet ~ Year + Pop + Locus + Ne + Year:Pop + Year:Locus
Model 2: ObsHet ~ Year + Pop + Locus + Ne + Year:Pop
Res.Df RSS Df Sum of Sq F Pr(>F)
1 51132 1530.5
2 51133 1530.5 -1 -0.010752 0.3592 0.549
> Ho_mod12 = lm(ObsHet ~ Year + Pop + Locus + Ne + Year:Pop, data = r12gNoOddsGenets)
> Ho_mod13 = lm(ObsHet ~ Year + Pop + Locus + Ne, data = r12gNoOddsGenets)
> anova(Ho_mod12, Ho_mod13)
Analysis of Variance Table
Model 1: ObsHet ~ Year + Pop + Locus + Ne + Year:Pop
Model 2: ObsHet ~ Year + Pop + Locus + Ne
Res.Df RSS Df Sum of Sq F Pr(>F)
1 51133 1530.5
2 51134 1530.5 -1 -0.029893 0.9987 0.3176
Ho_mod8
andHo_mod9
, the levels ofYear:Ne
are completely dependent on the predictors inHo_mod9
. That's why adding that term makes absolutely no difference. This is also referred to as a multicollinearity situation. If you look atsummary(Ho_mod8)
, you will find that a whole bunch of coefficients areNA
, meaning that those variables were omitted in fitting the model. $\endgroup$ – Russ Lenth Jun 5 '18 at 19:04