Estimating the standard errors of log-transformed response variables in Proc Mixed

[TD21]

Hey Folks:

I would like to know how to get the values of standard errors expressed in the original unit of log-transformed data that are analyzed using PROC MIXED. While the estimates of the means can be calculated using the anti-log of the results, I am not sure if the same procedure can be applied to calculate the standard errors since they are seemingly unitless (I could be wrong) based on the properties of logarithm. That is, Log x- Log y = Log (x/y). Thanks.

[PGStats]

On the original scale, standard errors are better expressed as a fraction of the means. Thus if SE=0.3 on the natural log scale, the std err is exp(0.3)-1  = 35% on the original scale.

PG

[SteveDenham]

You might consider doing the analysis in PROC GLIMMIX, with link=log, and using the ILINK option in the LSMEANS statement.  This will give the least squares mean and standard error on the original scale.  The standard error is calculated using the delta method.

Steve Denham

[TD21]

Hi Steve:

What are the equivalent PROC GLIMMIX codes for the following PROC Mixed Codes in order to get the estimates of the mean and standard error from the log (base 10) scale to the original scale? Thanks.

1.)

PROC SORT DATA = WORK.LOGKsat_MID_TRK_Pre_Post;

By Block Point Sample Depth LOC Harvest;

RUN;

ODS GRAPHICS ON;

PROC MIXED DATA = WORK.LOGKsat_MID_TRK_PrePost plots(only)=(ResidualPanel(marginal))PLOTS (MAXPOINTS = 2000000);

CLASS Block POINT SAMPLE DEPTH LOC Harvest;

MODEL LG_Ksat=Harvest|LOC|DEPTH/DDFM = Satterthwaite RESIDUAL;

RANDOM Block Block*Harvest BLOCK*LOC*Harvest;

RANDOM Depth/SUBJECT = Point(Block*LOC*Harvest) TYPE = AR(1);

LSMEANS Harvest Depth Harvest*Depth/PDIFF ADJUST = TUKEY SLICE=(Harvest Depth);

LSMEANS Harvest*LOC*Depth/PDIFF ADJUST = TUKEY SLICE=(Harvest LOC Depth);

RUN;

ODS GRAPHICS OFF;

2.)

PROC SORT DATA = Work.logKsat_Post2014;

BY Block Treatment LOC Point Depth Sample;

RUN;

ODS GRAPHICS ON;

PROC MIXED DATA =Work.logKsat_Post2014 plots(only)=(ResidualPanel(marginal)) PLOTS(MAXPOINTS= 2000000);

CLASS Block Treatment LOC Point  Depth Sample;

MODEL LG_Ksat = Treatment|Depth LOC(Treatment) Depth*LOC(Treatment) /RESIDUAL DDFM = Satterthwaite;

RANDOM Block Block*Treatment  Block*Treatment*LOC  ;

RANDOM Depth/SUBJECT = Point(Block*Treatment*LOC) TYPE = AR(1) ;

LSMEANS Treatment Depth Treatment*Depth/PDIFF ADJUST = TUKEY slice=(treatment depth);

LSMEANS  LOC*Depth(Treatment)/PDIFF ADJUST = TUKEY slice=(treatment depth LOC);

RUN;

ODS GRAPHICS OFF;

Thank you.

TD21

[SteveDenham]

For 1)

PROC GLIMMIX DATA = WORK.LOGKsat_MID_TRK_PrePost;* Not sure about plots in GLIMMIX, so this is commented out plots(only)=(ResidualPanel(marginal))PLOTS (MAXPOINTS = 2000000);

CLASS Block POINT SAMPLE DEPTH LOC Harvest;

MODEL Ksat=Harvest|LOC|DEPTH/DDFM = KR2 LINK=LOG;

OUTPUT out=residout resid=r;

RANDOM Block Block*Harvest BLOCK*LOC*Harvest;

RANDOM Depth/SUBJECT = Point(Block*LOC*Harvest) TYPE = AR(1);

LSMEANS Harvest Depth Harvest*Depth/ILINK ADJUST = TUKEY SLICEDIFF=(Harvest Depth) SLICEDIFFTYPE=ALL;

LSMEANS Harvest*LOC*Depth/ILINK ADJUST = TUKEY SLICEDIFF=(Harvest LOC Depth) SLICEDIFFTYPE=ALL;

RUN;

For 2)

PROC GLIMMIX DATA =Work.logKsat_Post2014;* see above plots(only)=(ResidualPanel(marginal)) PLOTS(MAXPOINTS= 2000000);

CLASS Block Treatment LOC Point  Depth Sample;

MODEL Ksat = Treatment|Depth LOC(Treatment) Depth*LOC(Treatment) /LINK=LOG DDFM = KR2;

OUTPUT out=residout resid=r;

RANDOM Block Block*Treatment  Block*Treatment*LOC  ;

RANDOM Depth/SUBJECT = Point(Block*Treatment*LOC) TYPE = AR(1) ;

LSMEANS Treatment Depth Treatment*Depth/ILINK ADJUST = TUKEY SLICEDIFF=(Treatment Depth) SLICEDIFFTYPE=ALL;

LSMEANS  LOC*Depth(Treatment)/ILINK ADJUST = TUKEY SLICEDIFF=(LOC TREATMENT DEPTH) SLICEDIFFTYPE=ALL;

RUN;

For this one, I see Sample as a CLASS variable, but not referenced in the MODEL or RANDOM statements.  My main concern is that perhaps it is needed to define the subject so that there is only one value for depth per point*block*treatment*loc combination.

Steve Denham

[TD21]

Thank you Steve.

[TD21]

Hi Steve:

I tried both 1 and 2, but I was not able to have a successful run due to the following: "Warning: obtaining minimum variance quadratic unbiased estimates as starting values for covariance parameters failed."

What do I need to change to avoid this error?

Thanks,

TD21

[SteveDenham]

We will probably have several iterations on fixing this error.  My first thought would be as follows.  Since you are fitting this as having a gaussian distribution with additive errors on the log scale, the marginal model should work.  Try changing the RANDOM statement that models the repeated nature in each of the two models to include the RESIDUAL option:

1)

RANDOM Depth/SUBJECT = Point(Block*LOC*Harvest) TYPE = AR(1) RESIDUAL;

2)

RANDOM Depth/SUBJECT = Point(Block*Treatment*LOC) TYPE = AR(1) RESIDUAL;

Hopefully, this will do the trick.  If not, I am afraid you might have insufficient data to support the complexity of the models you have chosen, so I will steal from one of the best - - and say try fitting something simpler, and see if you can "sneak up" on the solution.

Oh, and is the variable DEPTH equally spaced?  If not, AR(1) may not be a good candidate for the covariance structure.

Steve Denham

[cen]

Hi Steve

I tried this and compare the outcome for three approaches: 1) using link=log for the data transformation 2) compared to non-transformed data 3) compare to the result with taking log transformation before fitting the model.

Here is my finding: 1 and 2 approaches got exactly the same output.  1 and 3 doesn't agree (not what I have expected...).

**Orignial Variable with Link=log**;
ods graphics on;
proc glimmix data=sub method=laplace namelen=50 plots=STUDENTPANEL;
nloptions maxiter=100;  
    format &exp. &expf..;
    class &exp. (ref="&expref.") ;
    weight discwt;
    model los = &exp.  /solution link=log ;
    lsmeans &exp. /ilink;
title "";
run;
ods graphics off;

**Not transformed**;
ods graphics on;
proc glimmix data=sub method=laplace namelen=50 plots=STUDENTPANEL;
nloptions maxiter=100;  
    format &exp. &expf..;
    class &exp. (ref="&expref.") ;
    weight discwt;
    model los = &exp.  /solution;
    lsmeans &exp.;
title "";
run;
ods graphics off;

**Take log transformation before fit the model***;
proc glimmix data=sub method=laplace namelen=50;
nloptions maxiter=100;  
    format &exp. &expf.. ;
    class &exp. (ref="&expref.") ;
    weight discwt;
    model log_los = &exp.  /solution;
    lsmeans &exp.;
title "";
run;

Thanks

Steven