[adegenet-forum] Inbreeding calculation

[Jombart, Thibaut]

Hello, 

This is a scaling problem (I think I mentioned it before, now is probably a good time to clarify it). Densities are very small and when using them to get probabilities for 'sample' all probabilities are zero. You need to use functions directly. 

Note that I don't get your error, but informative warnings. I suspect your version of adegenet is outdated.

## what you've done ##

U1<- inbreeding(U, N=100)
There were 50 or more warnings (use warnings() to see the first 50)
> warnings()
Warning messages:
1: In FUN(X[[81L]], ...) :
  Likelihood uniformly zero likely reflecting precision issue
returning uniformly distributed sample
2: In FUN(X[[81L]], ...) :
  Likelihood uniformly zero likely reflecting precision issue
returning uniformly distributed sample
....


## the alternative ##
res <- inbreeding(U, res.type="function")

## look at the first function ##
plot(res[[1]]) # see, the y scale is 1e-245 (!!); still the graph tells you there's no inbreeding here

## look at the first 80 functions ##
par(mfrow=c(8,10), mar=rep(.1,4))
for(i in 1:80) plot(res[[i]], yaxt="n", col="red")


You'll see that a density can be estimated for most individuals. Flat bars or missing curves indicate NA problems.

If you really, really need to derive a sample for one given distribution, you have to scale the y axis:
####
x <- seq(0,1,le=1000)
mySamp <- sample(x, 100, prob=res[[1]](x)*1e250, replace=TRUE)
hist(mySamp, xlim=c(0,1))
####

But still, note that:
##
hist(mySamp, xlim=c(0,1))
##
is just an approximation of:
##
plot(res[[1]])
##

Cheers
Thibaut

--
######################################
Dr Thibaut JOMBART
MRC Centre for Outbreak Analysis and Modelling
Department of Infectious Disease Epidemiology
Imperial College - School of Public Health
St Mary’s Campus
Norfolk Place
London W2 1PG
United Kingdom
Tel. : 0044 (0)20 7594 3658
t.jombart at imperial.ac.uk
http://sites.google.com/site/thibautjombart/
http://adegenet.r-forge.r-project.org/
________________________________________
From: adegenet-forum-bounces at lists.r-forge.r-project.org [adegenet-forum-bounces at lists.r-forge.r-project.org] on behalf of Tohamy Yousef [tohamyy at yahoo.com]
Sent: 13 May 2014 17:18
To: adegenet-forum at lists.r-forge.r-project.org
Subject: [adegenet-forum] Inbreeding calculation

Dear Dr. Jombart
Firsst of all, I am sorry for disrupting you again. But, really I need your help. I wrote you before and asked about inbreeding calculation. You suggested that I use less number of markers and other advices. SO, I tried  to follow all your comments. For example I decided to use only 2247 SNP markers without missing data (before I had 43231 with missing data). But I am trying a lot and still have the same problem.
 I am attaching to you the new data (new_inbreeding.stru).


Thanks in advance.

######### my script

> library(adegenet)
> toto<-read.structure(file="new_inbreeding.stru",n.ind=154,n.loc=2247,onerowperind=FALSE,col.lab=1,row.marknames=1,NA.char="-9",missing="mean")

Which column contains the population factor ('0' if absent)? 2

 Which other optional columns should be read (press 'return' when done)? 1:
> toto

   #####################
   ### Genind object ###
   #####################
- genotypes of individuals -

S4 class:  genind
@call: read.structure(file = "inbreeding.stru",
    n.ind = 154, n.loc = 2247, onerowperind = FALSE, col.lab = 1,
    row.marknames = 1, NA.char = "-9", missing = "mean")

@tab:  154 x 3610 matrix of genotypes

@ind.names: vector of  154 individual names
@loc.names: vector of  2247 locus names
@loc.nall: number of alleles per locus
@loc.fac: locus factor for the  3610 columns of @tab
@all.names: list of  2247 components yielding allele names for each locus
@ploidy:  2
@type:  codom

Optionnal contents:
@pop:  factor giving the population of each individual
@pop.names:  factor giving the population of each individual

@other: - empty -
> toto$pop.names
  P1   P2
"US" "IP"
>U<- seppop(toto)$US
> U

   #####################
   ### Genind object ###
   #####################
- genotypes of individuals -

S4 class:  genind
@call: .local(x = x, i = i, j = j, treatOther = ..1, quiet = ..2, drop = drop)

@tab:  81 x 3610 matrix of genotypes

@ind.names: vector of  81 individual names
@loc.names: vector of  2247 locus names
@loc.nall: number of alleles per locus
@loc.fac: locus factor for the  3610 columns of @tab
@all.names: list of  2247 components yielding allele names for each locus
@ploidy:  2
@type:  codom

Optionnal contents:
@pop:  factor giving the population of each individual
@pop.names:  factor giving the population of each individual

@other: a list containing: elements without names

> U1<- inbreeding(U, N=100)
Error in sample.int(length(x), size, replace, prob) :
  NA in probability vector