Added parametric bootstrap for phylolm

[yuqing19118]

Also, I added the description of model for measurement error in the manual. I disabled lambda plus measurement error, because there are identifiability problems. If you want to see for yourself, Below is an example. You would need to re-enable lambda + measurement error to run this example (comment out lines 12 and 13 in phylolm.R). The example shows that the likelihood profile is flat where it's maximum.

first: simulate some data

set.seed(1238)
tre = rcoal(60)
taxa = sort(tre$tip.label)
b0=0; b1=1;
x <- rTrait(n=1, phy=tre,model="EB",
            parameters=list(ancestral.state=0,sigma2=1))
response <- b0 + b1*x + 
  rTrait(n=1,phy=tre,model="BM",parameters=list(ancestral.state=0,sigma2=1))
dat = data.frame(trait=response[taxa],pred=x[taxa])
z <- response+rnorm(60,0,.4)
dat = data.frame(trait=z[taxa],pred=x[taxa])

next: use phylolm with a fixed lambda, to get the profile likelihood at many fixed lambda values

lam = seq(0,1,by=.001)
loglk = rep(NA, length(lam))
s2e = rep(NA,length(lam))
cat("\n")
for ( i in 1:length(lam)){
  fit=phylolm(trait~pred,data=dat,phy=tre,model="lambda",measurement_error=TRUE,
                 lower.bound = lam[i], upper.bound = lam[i], starting.value = lam[i])
  loglk[i] = fit$logLik
  s2e[i] = fit$sigma2_error
  cat("done with replicate",i,"\r")
}
plot(loglk~lam, type="l") # full profile: lambda from 0 to 1
start=350 # next: zooming in the ML zone: flat.
plot(loglk[start:1001]~lam[start:1001], type="l")

the plot above showed the flat likelihood. The plot below shows how lambda hat varies with sigma_err hat:

plot(lam[start:1001], s2e[start:1001], type="l") # 

[lamho86]

Thanks. I would like to add you to the contributor list of phylolm. If you agree, please let me know your full name.

Best, Lam

[yuqing19118]

Hello Lam,

Thanks so much for your prompt response. My name is Qing Yu (Sabrina)

Thanks

Sent from my iPhone

On Jun 13, 2016, at 2:44 PM, Lam Ho notifications@github.com<mailto:notifications@github.com> wrote:

Thanks. I would like to add you to the contributor list of phylolm. If you agree, please let me know your full name.

Best, Lam

You are receiving this because you authored the thread. Reply to this email directly, view it on GitHubhttps://github.com/lamho86/phylolm/pull/8#issuecomment-225686915, or mute the threadhttps://github.com/notifications/unsubscribe/ASf_3Qq4T-0UKB09a_9NTnzrYflkvHqeks5qLbMvgaJpZM4I0eIb.