@@ -249,7 +249,6 @@ def ABC(contr_time=None, contr_volume=None, treat_time=None, treat_volume=None):
249249 return {"metric" : "abc" , "value" : abc ,'time' :np .max (treat_time )}#, "control": con, "treatment": tre}
250250
251251
252- ###LMM CODE
253252def lmm (time , volume , treatment , drug_name ):
254253 """
255254 Compute the linear mixed model (lmm) statistics for a PDX batch.
@@ -261,41 +260,21 @@ def lmm(time, volume, treatment, drug_name):
261260 dict: A dictionary containing the fit object and the specific coefficient value.
262261 """
263262
264- # --- DEBUG: inputs to LMM ---
265- print ("\n [DBG][LMM] incoming arrays" )
266- print (f"[DBG][LMM] len(time)={ len (time )} , len(volume)={ len (volume )} , len(treatment)={ len (treatment )} , drug_name={ drug_name } " )
267-
268-
269263 data = pd .DataFrame ({'model_id' :['model' ]* len (time ),\
270264 'volume' :volume ,\
271265 'time' :time ,\
272266 'exp_type' :treatment })
273267
274-
275- print ("[DBG][LMM] raw data shape:" , data .shape )
276- print ("[DBG][LMM] raw dtypes:\n " , data .dtypes )
277- print ("[DBG][LMM] NA counts:\n " , data .isna ().sum ())
278268
279269 data ['log_volume' ] = np .log (data ['volume' ])
280270 n_nonfinite_log = (~ np .isfinite (data ['log_volume' ])).sum ()
281- print (f"[DBG][LMM] non-finite log_volume count: { n_nonfinite_log } " )
282271
283272 # categories
284273 data ['exp_type' ] = data ['exp_type' ].astype ('category' )
285274 data ['exp_type' ] = pd .Categorical (data ['exp_type' ],
286275 categories = ['control' , drug_name ],
287276 ordered = True )
288- print ("[DBG][LMM] exp_type categories:" , list (data ['exp_type' ].cat .categories ))
289- print ("[DBG][LMM] exp_type value_counts:\n " , data ['exp_type' ].value_counts (dropna = False ))
290277
291- # time variation overall and by treatment
292- print (f"[DBG][LMM] time min/max: { data ['time' ].min ()} / { data ['time' ].max ()} " )
293- print ("[DBG][LMM] time describe by exp_type:\n " ,
294- data .groupby ('exp_type' , dropna = False )['time' ].describe ())
295-
296- # groups
297- print ("[DBG][LMM] unique model_id count:" , data ['model_id' ].nunique ())
298- print ("[DBG][LMM] model_id value_counts:\n " , data ['model_id' ].value_counts ())
299278
300279 ##create data frame from these 4 vectors
301280 required_columns = ["model_id" , "volume" , "time" , "exp_type" ]
@@ -315,15 +294,10 @@ def lmm(time, volume, treatment, drug_name):
315294 #print(data['exp_type'].cat.categories)
316295 # Fit the model
317296 model = mixedlm (formula , data , groups = data ['model_id' ])
318- print ("[DBG][LMM] exp_type counts:\n " , data ['exp_type' ].value_counts (dropna = False ))
319- print ("[DBG][LMM] time min/max:" , data ['time' ].min (), data ['time' ].max ())
320- print ("[DBG][LMM] volume<=0:" , (data ['volume' ]<= 0 ).sum (), " nonfinite log_volume:" , (~ np .isfinite (np .log (data ['volume' ]))).sum ())
321- print ("[DBG][LMM] groups:" , data ['model_id' ].nunique ())
322297
323- print ( "[DBG][LMM] DATA:" , data )
298+
324299 fit = model .fit ()
325300
326- print ("[DBG][LMM] fit:" , fit .params )
327301
328302 # Get the coefficient for the interaction term 'time:exp_type'
329303 #interaction_term = 'time:exp_type'
@@ -332,7 +306,6 @@ def lmm(time, volume, treatment, drug_name):
332306 #print(fit.params)
333307 i_coef_value = fit .params ['time:exp_type[T.' + drug_name + ']' ]
334308
335- print ("[DBG][LMM] i_coef_value:" , i_coef_value )
336309 #i_coef_value = fit.params['time:exp_type['+drug_name+']']
337310 # else:
338311 # coef_value = None # Handle the case when the interaction term is not present
@@ -423,24 +396,8 @@ def get_drug_stats(df, control='control'):
423396 else :
424397 singleres .append (treat_abc )
425398
426- # --- DEBUG: before LMM ---
427- print ("\n [DBG] -------- Drug block --------" )
428- print (f"[DBG] experiment={ name [0 ]} model_id={ mod } drug={ d } " )
429- print (f"[DBG] control rows: { len (ctl_data )} treated rows: { len (d_data )} " )
430-
431- # Quick sanity on time / volume
432- print (f"[DBG] control time min/max: { ctl_time .min () if len (ctl_time )> 0 else None } / { ctl_time .max () if len (ctl_time )> 0 else None } " )
433- print (f"[DBG] treat time min/max: { treat_time .min () if len (treat_time )> 0 else None } / { treat_time .max () if len (treat_time )> 0 else None } " )
434-
435- print (f"[DBG] control volume<=0: { (ctl_volume <= 0 ).sum () if len (ctl_volume )> 0 else None } " )
436- print (f"[DBG] treat volume<=0: { (treat_volume <= 0 ).sum () if len (treat_volume )> 0 else None } " )
437-
438- # Show the first few rows that will go into LMM
439399 comb = pd .concat ([ctl_data , d_data ])
440- print ("[DBG] comb dtypes:\n " , comb .dtypes )
441- print ("[DBG] comb treatment counts:\n " , comb ['treatment' ].value_counts (dropna = False ))
442- print ("[DBG] comb head:\n " , comb [['treatment' ,'time' ,'volume' ]].head ())
443- # --- END DEBUG ---
400+
444401
445402 lmm_res = lmm (comb .time , comb .volume , comb .treatment , d )
446403 lmm_res .update ({'sample' : mod , 'drug' : d , 'time' : np .max (treat_time ), 'time_unit' : 'days' })
0 commit comments