0
$\begingroup$

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:

  1. Am I addressing problems with normality correctly (see below)? I will ask my questions surrounding this excluding $N_e$ from my models.
  2. 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 enter image description here

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)

enter image description here

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)

enter image description here

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. enter image description here

Ho_compa = lm(HoLogit ~ Year + Pop + Locus + Year:Pop + Year:Locus+ Pop:Locus+ Year:Pop:Locus, data = r12gNoOddsGenets)

enter image description here

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
|cite|improve this question
$\endgroup$
  • $\begingroup$ Thank you for the response. In this case, do you have any thoughts about whether there is a statistical approach that I can use to address the question differently in this case? GLMMs, different distributions/intercepts? Letting locus be a random effect? And, do you have any thoughts on the zeros thing that is happening with that specific variable combination during model selection? $\endgroup$ – Ella Bowles Jun 5 '18 at 15:02
  • $\begingroup$ In, e.g., the comparison of Ho_mod8 and Ho_mod9, the levels of Year:Ne are completely dependent on the predictors in Ho_mod9. That's why adding that term makes absolutely no difference. This is also referred to as a multicollinearity situation. If you look at summary(Ho_mod8), you will find that a whole bunch of coefficients are NA, meaning that those variables were omitted in fitting the model. $\endgroup$ – rvl Jun 5 '18 at 19:04
1
$\begingroup$

Given the lumpiness of the residuals, no response transformation in a model fitted using lm is going to overcome that. Your attempted logit transformation appeared to lessen the problem, but only because you discarded all the 0s and 1s, thus undoubtedly biasing the estimates.

From here, it seems that the logit may indeed be an effective modeling choice, but not the way you are going about it. And you mention $N_e$ several times as if we know what it is, and don't explain it.

I'm guessing that your response values are proportions of the form $x/N_e$, where $x$ is some count of a binary outcome, and that $N_e$ is the number of outcomes total for each observation. If that is the case, and under a model that says these outcomes are independent, I think a logistic regression model is more appropriate. It would be fitted using the glm function; something like:

mod <- glm(ObsHet ~ Year + Pop + Locus + Year:Pop + Year:Locus + 
    Pop:Locus + Year:Pop:Locus, 
    weights = Ne, family = binomial,
    data = r12gNoOddsGenets)

This does incorporate the logit transformation, in the sense that the logit is the default link function for a binomial model.

Once you settle on a model, you should use the emmeans package for post-hoc comparisons, rather than lsmeans, which is being deprecated.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ Thank you for this response. My apologies for not explaining Ne. Ne is it's own variable, and is a point estimate of effective population size. My Hobs (heterozygosity) estimates are completely separate from that. Think of these two questions of mine separately, without Ne as a part of the discussion about the best model type. Heterozygosity is a proportion bounded by 0 and 1, but is not a binary outcome. I think this means that the binomial is not an appropriate distribution. I think the gamma distribution is best, but is really still in development. Thoughts on a different distribution? $\endgroup$ – Ella Bowles Jun 5 '18 at 16:08
  • $\begingroup$ Gamma would not ne appropriate. Beta would be, because it is a distribution on (0,1). However, you have a large number of hard 0s and 1s, so you need something that can fit a 0- and 1-inflated beta model. Maybe somebody can give you advice. $\endgroup$ – rvl Jun 5 '18 at 16:21
  • $\begingroup$ Maybe gamlss can handle it. But emmeans support is minimal and flaky. $\endgroup$ – rvl Jun 5 '18 at 16:29
  • $\begingroup$ Oops, yes, Beta is what I meant. Thank you for the thoughts on this. $\endgroup$ – Ella Bowles Jun 5 '18 at 17:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.