[ADMB Users] next try: "Incompatible shapes in df1b2vector functionoperator +"

[Bolker,Benjamin Michael]

  [Continuing with my saga of trying to get a
nonlinear size-density functional response model with
random effects running ...]

  It does look like initial values for the random effects
vector are indeed required.  Is this actually documented somewhere
that I missed ... ??

  I got simple2.tpl to run OK, so I am not completely hopeless (maybe).

 Re Dave Fournier's previous message: I agree that there are potential issues with the scale of
the random effects/ bounds on the probabilities, but I think that in reality for this problem
the random-effects variance is quite small, so I'm hoping to get away
(using good starting values and bounding the random effects variance)
with using this form for now.  Or I will use posfun() and penalties to make things
behave. If it seems likely that this issue is causing my current problems I will
go ahead and try to address it now -- otherwise I would like to try to get
the model running with a highly constrained RE variance, then relax
the constraints a bit and see if I run into trouble.

  I have instrumented my TPL file a little bit so that it tells
me what calculations it is doing.  It makes a little bit of headway
(prints out finite values for the function value, claims to have
gotten to at least iteration 50 -- although if I'm reading the output
correctly some of the parameter values are negative despite being
bounded below at zero??? never mind ...)

   It merrily prints out my debug statements for the calculations
a few hundred times, then gets to 

inner maxg = 0.103391
  Inner second time = 0  Inner f = 11.6341
Newton raphson 1

  and then stops with 

Incompatible shapes in df1b2vector functionoperator +

  ??? hard to tell where this is coming from ... ???

  The attached data file is very small (only the first 10 records),
but I get similar results with the full data file.

  I could try simulating well-behaved data, but (a) these data are
pretty well-behaved anyway and (b) I have a feeling that it won't
make a difference (but will give it a shot if someone thinks it will).

   Any ideas?

  
  sincerely
    Ben Bolker



PIN file:
------------------
# "mccoypred5.pin" produced by pin_write() from ADMButils; Thu Aug  6 01:39:03 2009
# c 
 0.506685 

# d 
 4.08488 

# h 
 0.008644705 

# g 
 3.019749 

# sigma_c 
 0.1 

# u 
 0 0 0 0 0 0 


DAT file:
----------------  
# "mccoypred5.dat" produced by dat_write() from ADMButils; Thu Aug  6 01:39:03 2009
# nobs 
 10 

# nblock 
 6 

# killed 
 3 0 1 3 1 0 1 2 2 2 

# size 
 8.830333 9.155 7.991667 8.000833 8.140333 8.711167 16.5265 14.1508 18.8445 14.87525 

# initial 
 6 6 6 6 6 6 6 6 6 6 

# Z 
1 0 0 0 0 0
0 1 0 0 0 0
0 0 1 0 0 0
0 0 0 1 0 0
0 0 0 0 1 0
0 0 0 0 0 1
1 0 0 0 0 0
0 1 0 0 0 0
0 0 1 0 0 0
0 0 0 1 0 0

TPL file:

DATA_SECTION

  init_int nobs               // # of observations
  init_int nblock             // # of blocks
  init_vector killed(1,nobs)  // # killed per trial
  init_vector size(1,nobs)    // size of individuals
  init_vector initial(1,nobs) // starting density (# individuals)
  init_matrix Z(1,nobs,1,nblock) // random-effects model matrix

PARAMETER_SECTION

  init_bounded_number c(0.0,1.0,1)
  init_bounded_number d(0.0,10,1)
  init_bounded_number h(0.0,1.0,1) 
  init_bounded_number g(0.0,10,1)
  init_bounded_number sigma_c(0.00001,0.1,1)
  random_effects_vector u(1,nblock,2)
  vector prob(1,nobs)
  vector cvec(1,nobs)
  vector a1vec(1,nobs)
  objective_function_value f

PROCEDURE_SECTION

  cout << "cvec calc\n";
  cvec = sigma_c*(Z*u)+c;
  a1vec = elem_prod(pow(size/d,g),exp(-size/d));
  cout << "prob calc\n";
  prob = 1/(1/(elem_prod(cvec,a1vec))+h*initial);

  cout << "loglik calc\n";
  f -= sum( log_comb(initial,killed)+
            elem_prod(killed,log(prob))+
            elem_prod(initial-killed,log(1-prob)));


GLOBALS_SECTION

  #include <admodel.h>
  #include <df1b2fun.h>
  #include <adrndeff.h>

  df1b2vector pow(const df1b2vector& v,const df1b2variable & x)
  {
    int mmin=v.indexmin();
    int mmax=v.indexmax();
    df1b2vector tmp(mmin,mmax);
    for (int i=mmin;i<=mmax;i++)
    {
      tmp(i)=pow(v(i),x);
    }
    return tmp;
  }