Phenotypic selection analysis.

[Sofia Ochoa]

Hello everyone,

I'm trying to use aster models to make phenotypic selection analysis. I want to see if the amount of defensive traits (HCN, DT and NEF) in different ontogenetic stages (P, J and R) reflects in fitness. The fitness measures that I'm using are Survival to reproduction (Rep), number of flowers (Flores), and number of seeds (sem).

I'm using the following graphical model:

1 --> survival to reproduction --> number of flowers --> number of seeds.

pred<-c(0,1,2)

survival to reproduction is Bernoulli, number of flowers and number of seeds are truncated Poisson = 0.

fam<-c(1,3,3)

The model I think works for me is:

out1<-aster(resp~varb+Edo:(HCN+NEF+DT),pred, fam, varb, id, root, data = redata)

The other option is to analyze each defensive attribute independently:

out2<-aster(resp~varb+HCN:Edo,pred, fam, varb, id, root, data = redata)

out3<-aster(resp~varb+NEF:Edo,pred, fam, varb, id, root, data = redata)

out4<-aster(resp~varb+DT:Edo,pred, fam, varb, id, root, data = redata)

In both cases I get the same error when asking for the summary of the models:

> out2<-aster(resp~varb+Edo:HCN,pred, fam, varb, id, root, data = redata2)

> summary(out2, show.graph = TRUE)

apparent null eigenvectors of information matrix

directions of recession or constancy of log likelihood

           [,1]       [,2]

[1,] -0.6465051 -0.0188301

[2,] -0.0266545  0.9996425

[3,]  0.7624439  0.0189800

[4,]  0.0000004 -0.0000001

[5,]  0.0000026 -0.0000004

[6,] -0.0001000  0.0000106

Error in summary.aster(out2, show.graph = TRUE) : 

  cannot compute standard errors

> out1<-aster(resp~varb+Edo:(HCN+NEF+DT),pred, fam, varb, id, root, data = redata)

> summary(out1, show.graph = TRUE)

apparent null eigenvectors of information matrix

directions of recession or constancy of log likelihood

            [,1]       [,2] [,3]

 [1,]  0.6461344  0.0180836    0

 [2,]  0.0253716 -0.9996755    0

 [3,] -0.7627954 -0.0179347    0

 [4,]  0.0000006 -0.0000001    0

 [5,] -0.0000208  0.0000035    0

 [6,]  0.0000952 -0.0000099    0

 [7,]  0.0000000  0.0000000    1

 [8,] -0.0030719  0.0004909    0

 [9,] -0.0000110  0.0000011    0

[10,]  0.0005539 -0.0001013    0

[11,]  0.0001449 -0.0000266    0

Error in summary.aster(out1, show.graph = TRUE) : 

  cannot compute standard errors

I don't know if my approach in using these models is correct, any thoughts would be greatly appreciated

Thank you

Sofia Ochoa

UNAM, Mexico

[geyer]

Firstly, I think you want fam to be c(1, 3, 2) rather than c(1, 3, 3).  See the post When to use zero-truncated Poisson or zero-truncated negative binomial at the top level.

Secondly, I do not understand why you do not want a formula resp ~ varb + seeds : something following slides for aster models course, deck 1, slides 57 ff.  If you follow that logic, then I suppose you want the formula resp ~ varb + seeds : Edo : (HCN + DT + NEF), but that really does have a direction of recession because there are no individuals with nonzero flower count and Edo == J and NEF != 0.  (Actually this seems to be a direction of constantcy so one parameter cannot be estimated).

But I don't understand why you aren't wanting any main effects for any of these parameters.

I await answers to this question before I try to do more.

By the way the model resp ~ seeds : (Edo : (HCN + DT) + NEF) does seem to not have any directions of recession, but I don't want to recommend that when I don't understand why no main effects.

[Sofia Ochoa]

Hi, 

Thanks for such rapid response.

I agree that the fam vector must be c(1,3,2).

In regards to the formula, what I want to test is if the change of defense among ontogenetic stages (Edo) affects the fitness. So, reading the slides I think that the model should be 

out2<-aster(resp~varb+seeds:Edo:(HCN+DT+NEF),pred, fam, varb, id, root, data = redata2)

However, when I run the model I get an error

out2<-aster(resp~varb+seeds:Edo:(HCN+DT+NEF),pred, fam, varb, id, root, data = redata2)

> summary(out2, show.graph = TRUE)

apparent null eigenvectors of information matrix

directions of recession or constancy of log likelihood

            [,1]       [,2] [,3]

 [1,]  0.6456795  0.0172959    0

 [2,]  0.0241342 -0.9997055    0

 [3,] -0.7631929 -0.0169906    0

 [4,]  0.0000016 -0.0000002    0

 [5,] -0.0000376  0.0000055    0

 [6,]  0.0001189 -0.0000105    0

 [7,]  0.0000099  0.0000016    0

 [8,]  0.0011696 -0.0001537    0

 [9,]  0.0002943 -0.0000391    0

[10,]  0.0000000  0.0000000   -1

[11,] -0.0071182  0.0010363    0

Error in summary.aster(out2, show.graph = TRUE) : 

  cannot compute standard errors

I try the model you suggest, and it runs right, but I don't understand how to interpret the summary.

> out4<-aster(resp~seeds:(Edo:(HCN+DT)+NEF),pred, fam, varb, id, root, data = redata2)

> summary(out4, show.graph = TRUE)

Call:

aster.formula(formula = resp ~ seeds:(Edo:(HCN + DT) + NEF), 

    pred = pred, fam = fam, varvar = varb, idvar = id, root = root, 

    data = redata2)

Graphical Model:

 variable predecessor family                           

 Rep      root        bernoulli                        

 Flores   Rep         truncated.poisson(truncation = 0)

 sem      Flores      poisson                          

                 Estimate Std. Error z value Pr(>|z|)    

(Intercept)     1.168e+00  3.745e-03 312.043  < 2e-16 ***

seeds:NEF      -1.563e-01  1.967e-02  -7.943 1.97e-15 ***

seeds:EdoJ:HCN  1.900e-04  9.708e-06  19.576  < 2e-16 ***

seeds:EdoP:HCN -1.841e-03  1.243e-04 -14.803  < 2e-16 ***

seeds:EdoR:HCN  1.156e-03  2.866e-05  40.320  < 2e-16 ***

seeds:EdoJ:DT  -4.316e-03  3.626e-04 -11.903  < 2e-16 ***

seeds:EdoP:DT   5.534e-02  1.354e-03  40.881  < 2e-16 ***

seeds:EdoR:DT   1.027e-02  2.299e-04  44.686  < 2e-16 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Thank you so much for the help.

Sofía O

[geyer]

Yes.  As I explained in my previous post, the error really is an error.  The first two directions of recession are spurious, but the third which is all zeros except for -1 for one component is a direction of constancy that the aster function does not pick up.  There are no individuals that have seeds, have EDO == "J", and have NEF > 0.  So an "interaction" between EDO and NEF makes no sense (cannot be estimated from these data).  That is why you cannot have a term EDO:NEF in your model.

As to how to interpret the summary from your out4, the same way you "interpret" any other.  First you want to compare with other simpler models to see whether an even simpler model "explains" the data just as well.  But I cannot help you with that yet because you didn't answer my question about why no "main effects".   In most applications one would expect to use a formula like

    resp ~ seeds : (Edo * (HCN + DT + NEF))

which would (among other terms) have "main effects"

    seeds:EdoJ

    seeds:EdiR

    seeds:EdoP

    seeds:HCN

    seeds:DT

which you apparently do not want for some scientific reason, but I don't know what that is.  Having these in the model would make it less likely that the "interactions" are statistically significant.

Also I should say that you and other users mostly don't understand the output of "summary" for any statistical models, not just aster models, because the R function summary is an invitation to misinterpretation.  The way to compare models is by likelihood ratio tests done by the R function anova.

> out.biggest <- aster(resp ~ varb + seeds : Edo : (HCN + DT + NEF),

+     pred, fam, varb, id, root, data = redata2)

> out.drop.nef <- aster(resp ~ varb + seeds : (Edo : (HCN + DT) + NEF),

+     pred, fam, varb, id, root, data = redata2)

> anova(out.drop.nef, out.biggest)

Analysis of Deviance Table

Model 1: resp ~ varb + seeds:(Edo:(HCN + DT) + NEF)

Model 2: resp ~ varb + seeds:Edo:(HCN + DT + NEF)

  Model Df Model Dev Df Deviance P(>|Chi|)

1       10    148148                      

2       11    148150  1    1.241    0.2653

This clearly shows that the term seeds : Edo : NEF is not statistically significant (the model without fits the data just as well as the model with it).  And that takes you down to the model you are wanting to explain.  Analogous tests for the other covariates

> out.drop.hcn <- aster(resp ~ varb + seeds : (Edo : (DT + NEF) + HCN),

+     pred, fam, varb, id, root, data = redata2)

> anova(out.drop.hcn, out.biggest)

Analysis of Deviance Table

Model 1: resp ~ varb + seeds:(Edo:(DT + NEF) + HCN)

Model 2: resp ~ varb + seeds:Edo:(HCN + DT + NEF)

  Model Df Model Dev Df Deviance P(>|Chi|)    

1        9    147823                          

2       11    148150  2   326.61 < 2.2e-16 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> out.drop.dt <- aster(resp ~ varb + seeds : (Edo : (HCN + NEF) + DT),

+     pred, fam, varb, id, root, data = redata2)

> anova(out.drop.dt, out.biggest)

Analysis of Deviance Table

Model 1: resp ~ varb + seeds:(Edo:(HCN + NEF) + DT)

Model 2: resp ~ varb + seeds:Edo:(HCN + DT + NEF)

  Model Df Model Dev Df Deviance P(>|Chi|)    

1        9    148060                          

2       11    148150  2   89.692 < 2.2e-16 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

show that the other "interaction" terms seeds : Edo : HCN and seeds : Edo : DT are highly statistically significant and cannot be dropped (models without them do not fit the data).

But I cannot stand behind those conclusions because I do not know why these are the scientifically interesting models to test!  And that is because of the unanswered question above

about "main effects".

The attached files have all the calculations.

[Sofia Ochoa]

Hi,

Sorry for the long hiatus but I was trying other models to analyze my data.

I think I finally understand my approach. We don't want to test for the main effects because we're not interested in kwon if defense traits (NEF, HCN or DT) or ontogenetic stage (P, J and R) alone have an effect on fitness, we are interested in the effect of the change of defense along ontogeny.

So, I interpret the analysis that you provided me in your last response as that the change in defense of trichome density and cyanogenetic potential have an effect over fitness, Is that accurate?

Thanks again for the help

[geyer]

OK but I still don't understand why your null model is not

    resp ~ varb + seeds : (Edo * (DT + NEF + HCN))

The point of the star (*) is to include main effects.  The fact that you are not interested in the main effects means you want them in your null model.  You would only leave them out if you assumed they were known to be exactly zero, and you don't (I wouldn't think).  And I apologize for taking a long time to answer.

[geyer]

Oops!  Typed too fast.  You would leave out the main effects only if you assumed the corresponding canonical parameters were known to be exactly zero.