[ADMB Users] adpool unit size error
Robert H Brown
rhbrown1 at netzero.com
Tue Jul 1 06:00:58 PDT 2014
Dear ADMB Users,
I am trying to fit a fairly simple 2-piece linear growth curve to some
longitudinal data with different curves in different subsets of the data
as defined by the design matrix. Hope to extend to random effects mixed
model. The tpl file compiles without error but when I try to run the
executable I get the message "error in adpool object you must set the
unit size". This message is written by the adpool.cpp. There is no
object name, just a blank. Not sure what this means or what to do.
Any suggestions would be appreciated. tpl file is attached.
Thank you.
Bob Brown
-------------- next part --------------
DATA_SECTION
init_int n // Number of patients
init_matrix obs(1,n,1,18) // Subject number, trtpn, fmgrp and Average weekly pain scores
init_matrix design(1,n,1,9) // Design matrix
vector op(1,n)
vector ds(1,n)
PARAMETER_SECTION
init_number acta1
init_number actb1
init_number acta2
init_number actb2
init_number acta3
init_number actb3
init_number acta4
init_number actb4
init_number actc
init_number plaa1
init_number plab1
init_number plaa2
init_number plab2
init_number plaa3
init_number plab3
init_number plaa4
init_number plab4
init_number plac
init_number s
init_number d1
init_number d2
objective_function_value g // Joint log-likelihood
TOP_OF_MAIN_SECTION
// arrmblsize = 40000000;
// gradient_structure::set_MAX_NVAR_OFFSET(450000);
GLOBALS_SECTION
#include <fvar.hpp>
#include <df1b2fun.h>
PRELIMINARY_CALCS_SECTION
op = column(obs,17);
ds = column(obs,18);
g = 0;
PROCEDURE_SECTION
dvariable s2, va2, vb2, mu, x, z, q, beta, a, b, c ; // Local variables
s2 = s*s;
beta = 0.90;
for(int i=1;i<=n;i++)
{
value(a) = value(design(i,2)*acta1 + design(i,3)*acta2 + design(i,4)*acta3 + design(i,5)*acta4 +
design(i,6)*plaa1 + design(i,7)*plaa2 + design(i,8)*plaa3 + design(i,9)*plaa4);
value(b) = value(design(i,2)*actb1 + design(i,3)*actb2 + design(i,4)*actb3 + design(i,5)*actb4 +
design(i,6)*plab1 + design(i,7)*plab2 + design(i,8)*plab3 + design(i,9)*plab4);
value(c) = value((design(i,2) + design(i,3) + design(i,4) + design(i,5))*actc
+ (design(i,6) + design(i,7) + design(i,8) + design(i,9))*plac);
for(int j=1;j<=13;j++)
{ x = j-1;
z = x - c;
mu = d1*op(i) + d2*ds(i) + a + b*x;
if (value(z) <= value(-beta)) {
q = 0;
}
if ((value(z) > value(-beta)) & (value(z) < value(beta))) {
q = (z+beta)*(z+beta)/(4*beta);
}
if (value(z) >= value(beta)) {
q = z;
}
mu = mu - b*q;
if (value(obs(i,3+j)) > -99.0) {
g -= (0.5*log(s2) + 0.5*(obs(i,3+j) - mu)*(obs(i,3+j) - mu)/s2);
}
}
}
g *= -1;
More information about the Users
mailing list