JacobSimonsen Tracker

Re: Marginal likelihood function of a fragile survival model and the Laplace transform

JacobSimonsen — Fri, 19 Sep 2025 13:35:02 GMT

Hi Tom,

I prefer the textbook. It is most logical for me.

and btw, it is not gemini that came up with the solution gemini shows. Gemini has it from somewhere. You would also not say cite google for a result, even that you found it using google.

Best, Jacob

Re: Train and test split for the proc phreg command

JacobSimonsen — Wed, 04 Jun 2025 11:53:36 GMT

Train and test data is relevant for prediction. Then it is probably not a Cox model you want to fit, as the Cox model does not predict anything. It only estimate the hazard rate ratios.

Re: comparing disease rate of a population with a disease rate for a subset of that population

JacobSimonsen — Tue, 06 Aug 2024 06:39:31 GMT

I agree with @StatDave .

You can not test if the rates are similar when there is no information about the certainty of each rate in your table. That can be either in the form of counts (numerator and denominator), standard errors or confidence limits. Do you have any of these?

Re: Counting Process and Interactions

JacobSimonsen — Mon, 22 Apr 2024 18:03:55 GMT

I have not really understood what that confuse you. It should make no difference whether you use counting process style ( (entry exit)*event(0) ), or delayed entry style (using the entry as an option).

Also, it is not really an interaction term you have made. The x1*stop term just add an effect that is b multiplied on x1 multiplied on time. Then it test b = 0.

Re: Counting Process and Interactions

JacobSimonsen — Mon, 22 Apr 2024 10:50:17 GMT

Hi,

I think instead you should use "stop" directly in the model line : "Model (start stop) * event(0) = x1 x1*stop x2 x3/rl alpha=0.05;" and you will then estimate x1 as being a constant + something linear with time. And you can then test if there is a significant time dependent term.

- though this work only under assumption that if there is time dependency, then its linear.

Re: How does the covariance matrix gets the confidence intervals of regression coefficients?

JacobSimonsen — Fri, 05 Apr 2024 11:09:27 GMT

I assume that with "covariance matrix" you mean the covariance between the parameter estimates.

Then you find the standard errors by taking squareroot to the diagonal elements in the covariance matrix.

Then you find CL by calculating the estimates +/- 1.96 x std errors.

Re: Interaction term in the PHREG procedure

JacobSimonsen — Tue, 05 Dec 2023 08:41:20 GMT

Hi Maria,

It is not necessary to have the "param=glm" as an option in the class statement. You will get the same estimates whether or not you use that option.

I though think it is easier to see what the reference group is then the option is used, as that will give a row with 0s.

Btw, I made a typo. It was the intention not to include the main effect of rx in the first program.

Jacob

Re: Interaction term in the PHREG procedure

JacobSimonsen — Thu, 30 Nov 2023 14:58:15 GMT

You can indeed add an interaction term. This is also possible even when one of the terms is used as a stratification variable as it is in your case.

You should just specify that enum variable is a class variable. Otherwise, phreg will just make regression on the product.

proc phreg data=survival.bladder2 covs(aggregate);
class enum rx/param=glm;
model futime*status(0)=rx rx*enum size number / ties=efron;
strata enum;
id id;
run;

And you can test the rx-effect is the same across the enum's by keeping in the main effect of rx. The test will come as in the table with type 3 tests.

proc phreg data=survival.bladder2 covs(aggregate);
class enum rx/param=glm;
model futime*status(0)=rx rx*enum size number / ties=efron;
strata enum;
id id;
run;

Good luck

Re: Proc Phreg for time varying exposure variable- different output

JacobSimonsen — Wed, 25 Oct 2023 11:57:44 GMT

Both can be correct, but most likely you should use (start, stop).

It has to do with how you define the baseline hazard. If the baseline hazard should be a function of time since 0 then use (start, stop).

If the hazard functions should be reset, such it is a functio of time since last stop time then you can use your model 2.

In far most cases, model 1 is the right way to go.

For your second question about how the hazard ratio is calculated, the simple answer is that it is the maximum point of the Cox's partial likelihood function. The hazard ratio has no closed analytical form. Except in simple cases with very few events. In your case the estimated hazard ratio is the maximum point of the function (x/(1+2*x))*(1/(1+x)) which is sqrt(1/2). I dont get exactly same number as you mention btw, but I have tried run the SAS code and the solution is indeed sqrt(0.5)=0.707.

Re: Cox regression vs Poisson regression for analysis

JacobSimonsen — Tue, 05 May 2020 10:28:47 GMT

If you assume piecewise constant hazard rates, then the likelihood function (as a function of parameteres) has same form as if the number of events had been poisson distributed. Therefore, you can use poisson regression on time to event data where you have counts on left side in the model statement. If you chop the timeaxis into finer and finer pieces, then the model will be equivalent to a cox-regression, and in that case the difference is only that the parameter of the time-effect is non-parametric in the cox-regression while it will be estimated together with other parametes in the Poisson regression model.

Be aware, there is no assumption that the counts actually are poisson distributed - and they are obviously not since the count are limited by the number of subjects in the trial, whereas if it was poisson distributed then there would not be an upper limit. Still, it is called poisson regression because of the similarity of the likelihood functions. Since data are not poisson distributed, it will not give meaning to apply model check that rely on the poison distribution ("variance = mean" does not give sense here).

Re: Visualizing an interaction between 2 continuous predictors w/confidence bands (Proc Mixed model)

JacobSimonsen — Wed, 25 Mar 2020 10:26:55 GMT

Hi again,

With some simpler simulated I was able to plot the effect with the proc plm using the effectplot with slicefit option.

*simulate a dataset with predictors c and x, and outcome y;
data simulation;
  do i=1 to 1000;
    c=int(rand('uniform',0,5));
    x=int(rand('uniform',0,5));
	y=rand('normal',(1/2-c)*x+x*x+2,1);
    output; 
  end;
run; 

*estimate the model and save the model to "mymodel";
proc mixed data=simulation;
  model y= c*x c c*x*x;
  store out=mymodel;
run;

*plot the prediction with confidence limits;
ods graphics;
ods hmtl;
proc plm restore=mymodel;
  effectplot slicefit(x=x sliceby= c=2,3) /clm;
run;

Since you have many more covariates in your model, all the other covariates will be set at some level. I hope you find it usefull.

Re: Visualizing an interaction between 2 continuous predictors w/confidence bands (Proc Mixed model)

JacobSimonsen — Tue, 24 Mar 2020 07:54:42 GMT

I think you are right by using "store" and proc plm. Can you provide your program?

Re: Comparing stratified logistic regression models to an unstratified model

JacobSimonsen — Mon, 23 Mar 2020 13:56:04 GMT

It depends on what you mean by "stratified model". If you mean that you run your model with or without a effect modifier on some covariate of interest, then you can indeed test if the effect modifier is significant (in terms of p-value). It is just to test if there is an interaction effect.

But, If stratified model means that you run a conditional logistic regression model by using the strata-statement, then I don't see any way to test if the variable used in the strata-statement change the result significant (in terms of p-value). You will have to judge if you find the strata variable has changed the result so much that your find the change "significant" (but then done use the word significant as it indicate you use a p-value). Also, It will most often not make sense to remove the variable in the strata statement, as it is given by the design of the model that the strata-variable has to be there (as in forexample a 1:k sampled case control study) .

Re: trying to check proportional hazards for cox model

JacobSimonsen — Tue, 17 Dec 2019 14:21:27 GMT

Hello,

you are doing some mistakes.

1) it doesnt give meaing to multiply the character variable sex with time. That should be replaced with something that can intepretate and translate characters into numeric. Try this
timeXsex=time*(sex='Male');

2) This is the most important! It doesnt give meaing to make the interaction terms between covariates and the time variable in the dataset. These interaction terms should be made inside phreg. The reason for this is that interaction should based on the time at each riskset, and not the time for event at each observation. it will be something like this (simplified)


proc phreg data=simulation;
  class sex;
  timeXsex=time*(sex='Male');
  model time=exposure sex timeXsex;
run;

3) Unless you have very low statistical power, I dont see any reason for why you have all covariates in the model line, and none in the strata statement. Have the covariates in the strata statement gives a more robust model, and you will then then only need to check the proportional hazard assumption on the covariate of interest which should be in the model line.

Re: Model is too large to be fit in a reasonable amount of time - Glimmix

JacobSimonsen — Tue, 17 Dec 2019 13:03:55 GMT

I agree with @Ksharp . It will be impossible to do with glimmix, because this procedure will estimate a parameter for each family.

It should be possible to estimate with generelized estimation equations (genmod+repeated) or GEE, because this method works in a non-parametric way.

I will though not use type=un, but rather type=cs, as cs assume same correlation between any two individuals within a familiy.

Re: Negative Binomial and Incidence

JacobSimonsen — Sat, 09 Nov 2019 20:18:59 GMT

My suggestion is to use the poisson regression model - no matter whether you see overdisperson or not.

The reason for that that overdispersion is something that measure variance relative to mean for count-data. But here the data was time-to-event, which has been summarized by number of events and person-time. It can be analyzed with poisson regression assuming piecewice constant incidence rates, because the likelihood function for will be similar to what you would have if you assumed the count-data was poisson distributed. As you dont need the assumption of distribution it you also dont need to check for overdispersion.

It quite easy to simulate time to event data with some incidence rate, and estimate the incidencerate with confidence interval correctly, even that you may have a big over dispersion. The dispersion will depend on how much you have aggreated (for example 5 agegroup instead of 10), which can be completely arbitrary and irrelevant for the rate estimate and it std error. So, just use proc genmod with the dist=poisson, and remember the offset=logtime.

Re: Multipass option in PHREG

JacobSimonsen — Thu, 07 Nov 2019 13:55:49 GMT

Yes, the Cox-aggregate macro can also be used instead of the fast-option. It willdo the aggregation on riskset. The weigt statement in phreg should then be used to specify how many person there are at risk. Examples of how to do this is below the macro.

It will make a much smaller utility file so you shouldnt get problems with disk-space.

I dont see any reason for experimenting with using multipass or not, because, if the utility file can be created once, then it almost sure can be recreated again for each iteration. It may save a bit time to use it, or not use it, but it will not do the major difference. The major difference is obtained by aggregating.

Re: Multipass option in PHREG

JacobSimonsen — Tue, 05 Nov 2019 07:49:30 GMT

The problem is that PROC PHREG will not per default make aggregation on each riskset when time-dependent variables are present. Its also not always that it will be possible. Instead, it for each event go evaluate every individual whether it is at risk at that eventtime or not. That makes the running time quadratic relative to the cohort size. If multipass option is used, then it do this exercise again for each iteration.

If your timedependent covariate is piecewise constant, then you can chop your data into intervals. Such that you will have both an entry and exittime variable. Then specify the "fast" option to tell PROC PHREG that it should use the aggregetion technique when it calculate the likelihood function. The calculation time will then only be linear relative to the cohort size, and you will probably be surprised how fast go.

Good luck.

Re: Wald Chi Square statistics - Logistic Regression

JacobSimonsen — Tue, 09 Jul 2019 12:02:12 GMT

You can use wald statistics, and likelihood ratio test that have asymptotically chi-squared distributions in linear regression. But, when data is normal distributed, then it is possible to use the exact distributions (not relying on asymptotic results). Therefore, you use t-statistics and F-test in linear regression as it is more exact. Actually, if you use proc genmod instead of proc glm/proc mixed for normal distributed data then you will get the wald and chi-square statistics.

In logistic regression it is not possible (or in best case very difficult) to find test statistics with a known exact distribution, therefore you use chi-square and wald statistics because then you at least know their asymptotic distribution. And actually, n doesnt need to be very large before the chi-square statistic are practically indistinguishable from a χ²distribution.

Re: proc mcmc: An exception has been encountered - Read Access Violation - Write Access Violation

JacobSimonsen — Mon, 06 May 2019 11:30:54 GMT

The error comes because of your way of reading in data. You dont give the procedure your data in the usual way (data=..). Instead, you attemp to read in data with the read_array function. According to the documentation, the read_array function can only be used inside a function or call routine, which in means that you only can use read_array inside proc fcmp.

This is from the documentation:

"The READ_ARRAY function attempts to dynamically resize the array to match the dimensions of the input data set. This means that the array must be dynamic. That is, the array must be declared either in a function or CALL routine or declared with the /NOSYMBOLS option. "