title: How to code a multiparameter log-likelihood function in R
body: <p>I would like to estimate power of the following problem. I am interested in comparing two groups that both follow Weibull distribution. So, group A has two parameters (shape par = a1,scale par = b1) and two parameters has group B (a2, b2). By simulating random variables from distribution of interest (for example assuming different scale and shape parameters, i.e. a1=1.5*a2, and b1=b2*0.5; or either differences between groups are just in either shape or scale parameters), apply log-likelihood ratio test to test if a1=a2 and b1=b2 (or e.g. a1=a1, when we know that b1=b2), and estimate power of the test.</p>
  
  <p>The questions would be what are log-likelihoods for the full models, and how to code it in R when 
  a)having exact data, 
  and b) for interval-censored data ?</p>
  
  <p>That is, for reduced model (when a1=a2,b1=b2) log-likelihoods for exact and interval-censored data are:</p>
  
  <pre><code>LL.reduced.exact &lt;- function(par,data){sum(log(dweibull(data,shape=par[1],scale=par[2])))};
  LL.reduced.interval.censored&lt;-function(par, data.lower, data.upper) {sum(log((1-pweibull(data.lower, par[1], par[2])) – (1-pweibull(data.upper, par[1],par[2]))))}
  </code></pre>
  
  <p>What is it for the full model, when a1!=a2, b1!=b2, taking into account two different observational schemes, i.e. when 4 parameters have to be estimated (or, in case when interested in looking at diferences in shape parameters, 3 parameters have to be estimated)? </p>
  
  <p>Is it possible to estimate it buy building two log-likelihoods for separate groups and add it together (i.e.<strong>LL.full&lt;-LL.group1+LL.group2</strong>)? </p>
  
  <p>Regarding log-likelihood for interval-censored data, censoring is non-informative and all observations are interval-censored. Any better ideas how to perform this task will be appreciated.</p>
  
  <p>Please, find the R Code for exact data below to illustrate the problem. Thank you very much in advance.</p>
  
  <pre><code>R Code:    
  # n (sample size) = 500
  # sim (number of simulations) = 1000
  # alpha  = .05
  # Parameters of Weibull distributions: 
     #group 1: a1=1, b1=20
     #group 2: a2=1*1.5 b2=b1
  
  n=500
  sim=1000
  alpha=.05
  a1=1
  b1=20
  a2=a1*1.5
  b2=b1
  #OR: a1=1, b1=20, a2=a1*1.5, b2=b1*0.5 
  
  # the main question is how to build this log-likelihood model, when a1!=a2, and b1=b2
  # (or a1!=a2, and b1!=b2)
  LL.full&lt;-????? 
  LL.reduced &lt;- function(par,data){sum(log(dweibull(data,shape=par[1],scale=par[2])))}
  
  LR.test&lt;-function(red,full,df) {
  lrt&lt;-(-2)*(red-full)
  pvalue&lt;-1-pchisq(lrt,df)
  return(data.frame(lrt,pvalue))
  }
  
  rejections&lt;-NULL
  
  for (i in 1:sim) {
  
  RV1&lt;-rweibull (n, a1, b1)
  RV2&lt;-rweibull (n, a2, b2)
  RV.Total&lt;-c(RV1, RV2)
  
  par.start&lt;-c(1, 15)
  
  mle.full&lt;- ????????????  
  mle.reduced&lt;-optim(par.start, LL, data=RV.Total, control=list(fnscale=-1))
  
  LL.full&lt;-????? 
  LL.reduced&lt;-mle.reduced$value
  
  LRT&lt;-LR.test(LL.reduced, LL.full, 1)
  
  rejections1&lt;-ifelse(LRT$pvalue&lt;alpha,1,0)
  rejections&lt;-c(rejections, rejections1)
  }
  
  table(rejections)
  sum(table(rejections)[[2]])/sim   # estimated power
  </code></pre>