Calculate half life from linear regression

[mariko5797]

I want to calculate half-life by fitting a linear regression model to log transformed PRNT titers, then obtain p-values using Mann-Whitney U Test. I'm not confident in my methodology, and would appreciate some verification/feedback.

1) Most of what I've seen in calculating half-life use natural log-transformations. Can log10 transformations be used instead?

2) If so, then the model equation would be log10(AVAL)=log10(AVAL0)+k*t, where AVAL is the PRNT results, AVAL0 is the PRNT results at baseline, k is the reaction rate constant, and t is the time (days). 

3) Then the slope would be k? So half-life can be calculated using t_1/2=log10(0.5)/k. 

4) Is the following code an accurate way to calculate half-life and p-value?

/*H0: Humoral immune response in Group 1 as assessed by PRNT half-life will be similar to Group 2*/

data have;
 label 	USUBJID	= "Unique Subject Identifier"
 		ACTARMN	= "Actual Treatment Arm (N)"
		ATPTN	= "Analysis Time Point (N)"
		AVAL	= "Analysis Value (PRNT Result)"
		;
 input USUBJID ACTARMN ATPTN AVAL @@;
 cards;
 1	1	1	10 	1	1	15	49 	1	1	29	79 
 1	1	43	271	1	1	57	216	1	1	90	142
 1	1	181	86 	1	1	365	20 
 2	1	1	12 	2	1	15	54 	2	1	29	89 
 2	1	43	267	2	1	57	209	2	1	90	131
 2	1	181	80 	2	1	365	26 
 3	1	1	21 	3	1	15	60 	3	1	29	101
 3	1	43	282	3	1	57	220	3	1	90	151
 3	1	181	85 	3	1	365	29 
 4	2	1	17 	4	2	15	55 	4	2	29	90 
 4	2	43	291	4	2	57	232	4	2	90	160
 4	2	181	100	4	2	365	32 
 5	2	1	17 	5	2	15	59  5	2	29	99 
 5	2	43	286	5	2	57	203	5	2	90	151
 5	2	181	91 	5	2	365	18 
 6	2	1	20 	6	2	15	56 	6	2	29	92 
 6	2	43	288	6	2	57	226	6	2	90	151
 6	2	181	90 	6	2	365	18 
 ;
data have;
 set have;
 	log_AVAL = log10(AVAL);
run;

/*linear regression model*/
proc reg data= have plots= none;
 by USUBJID ACTARMN;
 model log_AVAL=ATPTN;
 ods output ParameterEstimates= reg_mod;
run;

/*mean half-life*/
data reg_mod;
 set reg_mod;
	if Variable = 'ATPTN'	then T_HALF = log10(0.5)/Estimate;
proc means data= reg_mod(where= (Variable='ATPTN'));
 class ACTARMN;
 var T_HALF;
 output out= t_half mean= Mean_T_HALF min= Min_T_HALF max= Max_T_HALF;
run;

/*Mann-Whitney U Test*/
proc npar1way data= reg_mod wilcoxon;
 class ACTARMN;
 var T_HALF;
 ods output WilcoxonTest= pvalue;
run;

[sbxkoenk]

Hello,

I am not very much at ease in the Pharmaco-Kinetics (PK) world , but have you seen below papers?

• PharmaSUG 2023 - Paper SA-284
  SAS-based Method for PK Noncompartmental Analysis and Validation
  Hui Mao, BioPier Inc.
  Lixin Gao, BioPier Inc.
  https://www.pharmasug.org/proceedings/2023/SA/PharmaSUG-2023-SA-284.pdf
• PhUSE 2012 -- Paper SP08
  Estimating terminal half life by non-compartmental methods with some data below the limit of quantification
  Jochen Müller-Cohrs, CSL Behring, Marburg, Germany
  https://www.lexjansen.com/phuse/2012/sp/SP08.pdf

Calling @kessler (M.)

Koen

[mariko5797]

The following code was pulled from a previous study. Note, the analysis value has been transformed to LOG2 instead of LOG10. The issue with this method is it does not give me individual half-life results, i.e., for each subject. I need to find the minimum, maximum, and p-value (Mann Whitney U) of the treatment groups. Is there a way to tweak the program to obtain that?

proc genmod data= DSNIN; 
 by ACTARMN;
 class USUBJID;
 model log2_AVAL=ATPTN	/ dist= normal link= identity corrb type3;
 ods output GEEEmpPEst= Reg_Mod;
 repeated subject= USUBJID / corr= ar(1) corrw;
run;

proc sql;
 create table Results as
 	select 'PRNT' as ASSAY, ACTARMN, -1*(1/Estimate) as T_HALF, -1*(1/LowerCL) as T_HALF_LB, -1*(1/UpperCL) as T_HALF_UB
		from Reg_Mod
		where PARM= 'ATPTN';
quit;

[SteveDenham]

Just some miscellaneous comments - it doesn't make any difference what base log transform is used, the estimate for t1/2 should be the same. This is because the various log values differ only by a multiplicative constant, which factors out. 

Next, log transforming the data makes some assumptions regarding the variance (that the variance is proportional to the mean, for instance). Consequently, this could be done on untransformed data using PROC NLIN or PROC NLMIXED. By including an indicator variable and some clever programming steps, you could get an estimate of the difference in t1/2 between treatments using an ESTIMATE statement in NLMIXED. If you are interested in this approach, let us know.

Finally, I had no trouble running your code. I assume that the method is historically what has been done in your field, so I would stay with it.

SteveDenham

[mariko5797]

Thank you! I would be interested in your proposed approach, so I can better assess my options.

The method used here, isn't standard (to my knowledge). I created the program based on my own understanding of how half-life is calculated.

[SteveDenham]

So I started to work on this with NLMIXED when I noticed something critical - these data do not describe a single pool decaying process of the type  c = c0 * exp (-k * time). It appears that at time=0 the response value is at or near zero, then increases to about time = 43, followed by a decrease. There are two logical models that fit this - a time independent two-pool model with sampling from the second pool, so that there is an entry rate constant, or a time-dependent model with sampling from a single pool.  

Neither of those models fit the linearized (log transform) approach you have. One thing that might help would be to estimate the removal coefficient and estimate t1/2 from that. That would be easiest by removing time points less than 43 from the analysis and rescaling the time axis. However, that estimate of t1/2 will be biased by not having an accurate initial estimate of the maximum concentration and when it occurs.  However, in many pharmacokinetic studies, this is the standard approach, and it gives a good estimate of the removal coefficient and the resulting half-life in the pool.  Here is code to accomplish that:

proc nlmixed data=have2(where=(atptn>40)) maxiter=5000  ;
    parms
        c0 =  1.5790
        d1 =  0
        k0  = -0.5291
        d2 = 0    
        s2       =  0.2069
    ;

    pred   = c0*trt1*exp(k0*trt1*(ATPTN-43)) + (c0 + d1)*trt2*exp((k0 + d2)*trt2*(ATPTN-43));
    

    model aval ~ normal(pred,s2);
    predict pred out=pred;
estimate 'initial for trt1' c0 ;
estimate 'removal coefficient for trt1' k0;
estimate 'initial for trt2' c0 + d1;
estimate 'removal coefficient for trt2' k0 + d2;

estimate 'half-life for trt1' -log(2)/k0;
estimate 'half-life for trt 2' -log(2)/(k0 + d2);
    ods output parameterestimates=parms additionalestimates=addl;
run;

I think the use of the CMPTMODEL statement would enable fitting all of the data, but I will have to learn more about it before I post anything.

SteveDenham

[Rick_SAS]

> The method used here, isn't standard (to my knowledge). I created the program based on my own understanding of how half-life is calculated.

The method is standard. It is the same idea as for estimating the doubling time of an exponentially increasing function, but here you are estimating the "halving time" for a decreasing function. For the mathematics behind estimating doubling time, see

https://blogs.sas.com/content/iml/2020/04/01/estimate-doubling-time-exponential-growth.html

I think you will see that the math behind the two methods are similar, except for the factor of 2 vs 1/2.

[SteveDenham]

Just be sure that you use the monotonically decreasing part of these curves to estimate the half-life if you implement the method in @Rick_SAS 's blog. Including the increasing phase will lead to a best fit line that has a much, much smaller slope than the part of the curve dominated by elimination from the pool sampled. As a result, the half life will be greatly over-estimated.

SteveDenham 

PS I am working through the paper by Kurada & Chen to use compartment modeling on this data. Here is a link: https://support.sas.com/resources/papers/proceedings18/1883-2018.pdf . Then on to  the paper by Hannion and Boulanger https://www.lexjansen.com/phuse/2020/as/PRE_AS06.pdf . From these, all of the data can be used, not just the monotonically decreasing data.

[Rick_SAS]

Steve: your point is well-taken: You want to make sure you are modeling a steady-state phenomenon and do not pick up any initial transitory behavior. 

Unless the data perfectly fits an exponential distribution, the estimate depends on the location at which you apply it. In my blog, I used the five most recent records to estimate the slope to the model at the most recent data point. Other estimates of the slope will yield different results. For example, if you use only the two most recent data points, you will get an estimate that has more variance.

[mariko5797]

Thanks Steve.
What if the measurements are not monotonically decreasing? Is there another method?
I tested on the actual data I have, and while most had reasonable estimates, some were way larger than expected, especially those consistently <LLOD (imputed 1/2 LLOD).

[SteveDenham]

LLOD = can of worms.

There is a lot of statistical literature out there on how to handle values that might be left-censored at a lower limit of quantitation, but few that apply to LLOD, as that can be defined as the value below which measurement noise is greater than the signal. I would impute values below LLOD as coming from a uniform distribution on [0, LLOD], which happens to have an expected value of LLOD/2. So, your imputation is probably not the issue. What I suspect is that the plot of the data has a hockey-stick appearance, with a long flat area out to the left, where most of the values are either very small or imputed. If you fit a line out there, it is going to have a slope approaching zero, with a relatively large residual SD. When you plug that into the back calculation for t1/2 , you get a really large value.

But think about it, the concentrations out in that area are due to non-single compartment behavior, either slow release from a precursor compartment or some sort of flip-flop kinetics.  What my colleagues in the pharmacokinetics group do is to examine the time course curve, and only use the data from the descending portion prior to entering LLOD territory. They set values <LLOD to zero, and any values after the first zero are also set to zero. That more or less restricts the analysis to a single compartment with a single elimination rate constant that is time-independent. I can see their point - I would rather estimate the rate constant with 4 points descending than 104 points where the last 100 are effectively zero.

Now if your time-concentration curve isn't nicely behaved, such as having two (or more) peaks, long flat periods, or increasing concentrations well after initial dose, then noncompartmental analysis may not give you a believable half-life.

SteveDenham