The script begins by creating a 'pump' dataset containing information on pump reliability (number of failures 'y', operating time 't', and 'group'). Subsequently, two statistical procedures are employed: 'PROC GLIMMIX' to fit a generalized linear mixed model with a Poisson distribution and a logarithmic link function, including a random effect per pump and estimating differences in intercepts and slopes between groups. The second procedure, 'PROC NLMIXED', fits a nonlinear mixed model for the same data, defining a conditional linear predictor 'eta' and a rate parameter 'lambda' for a Poisson distribution, with a normally distributed random effect and estimates for differences between group parameters and the variance of the random effect. The entire analysis is generated in HTML format using ODS.
Data Analysis
Type : CREATION_INTERNE
The 'pump' dataset's data is directly created and integrated into the script via a DATALINES statement.
1 Code Block
DATA STEP Data
Explanation : This DATA STEP block creates the 'pump' dataset. 'pump' is a unique identifier based on the observation number. 'y' represents the number of failures, 't' the operating time, and 'group' is a categorical variable. 'logtstd' is a transformed variable (logarithm of 't' minus a constant), potentially standardized or centered. The data is provided directly in the script via the DATALINES statement.
Explanation : This block uses 'PROC GLIMMIX' to fit a generalized linear mixed model to the 'pump' data. The quadrature integration method ('method=quad') is specified. The 'group' variable is declared as a classification variable. The model specifies 'y' as the response variable, 'group' and the 'logtstd*group' interaction as predictors. The 'noint' option suppresses the automatic intercept, 'dist=poisson' indicates a Poisson distribution for the response, and 'link=log' uses a logarithmic link function. 'ddfm=residual' specifies the method for calculating degrees of freedom. A random effect ('random int') is included for each 'pump'. Finally, estimates are made for the differences between the intercepts and slopes of the two groups. ODS commands open and close the HTML output.
Copied!
ods html;
title "Pump Reliability at a Pressurized Water Reactor Nuclear Power Plant";
proc glimmix data=pump method=quad;
class group;
model y = group logtstd*group / noint dist=poisson link=log
dfml=residual;
random int / subject=pump;
estimate "Difference Y-intercepts" group 1 -1;
estimate "Difference Slopes" logtstd*group 1 -1;
run;
ods html close;
1
ods html;
2
title "Pump Reliability at a Pressurized Water Reactor Nuclear Power Plant";
3
PROC GLIMMIXDATA=pump method=quad;
4
class group;
5
model y = group logtstd*group / noint dist=poisson link=log
6
dfml=residual;
7
random int / subject=pump;
8
estimate "Difference Y-intercepts" group 1 -1;
9
estimate "Difference Slopes" logtstd*group 1 -1;
10
RUN;
11
ods html close;
3 Code Block
PROC NLMIXED
Explanation : This block uses 'PROC NLMIXED' to fit a nonlinear mixed model. Initial parameters ('parms') are defined for 'logsig', 'beta1', 'beta2', 'alpha1', and 'alpha2'. The 'eta' variable is conditionally defined based on the 'group' value, incorporating fixed parameters ('alpha1', 'alpha2', 'beta1', 'beta2') and a random effect 'e'. 'lambda' is calculated as the exponential of 'eta'. The model specifies that 'y' follows a Poisson distribution with the parameter 'lambda'. The random effect 'e' is assumed to follow a normal distribution with a mean of 0 and a variance of 'exp(2*logsig)', with 'pump' as the subject for the random effects. Estimates are provided for the differences between the 'alpha' and 'beta' parameters of the groups, as well as for the variance ('Sigma**2') of the random effect.
Copied!
proc nlmixed data=pump;
parms logsig 0 beta1 1 beta2 1 alpha1 1 alpha2 1;
if (group = 1) then eta = alpha1 + beta1*logtstd + e;
else eta = alpha2 + beta2*logtstd + e;
lambda = exp(eta);
model y ~ poisson(lambda);
random e ~ normal(0,exp(2*logsig)) subject=pump;
estimate 'alpha1-alpha2' alpha1-alpha2;
estimate 'beta1-beta2' beta1-beta2;
estimate 'Sigma**2' exp(2*logsig);
run;
1
PROC NLMIXEDDATA=pump;
2
parms logsig 0 beta1 1 beta2 1 alpha1 1 alpha2 1;
3
IF (group = 1) THEN eta = alpha1 + beta1*logtstd + e;
4
ELSE eta = alpha2 + beta2*logtstd + e;
5
lambda = exp(eta);
6
model y ~ poisson(lambda);
7
random e ~ normal(0,exp(2*logsig)) subject=pump;
8
estimate 'alpha1-alpha2' alpha1-alpha2;
9
estimate 'beta1-beta2' beta1-beta2;
10
estimate 'Sigma**2' exp(2*logsig);
11
RUN;
This material is provided "as is" by We Are Cas. There are no warranties, expressed or implied, as to merchantability or fitness for a particular purpose regarding the materials or code contained herein. We Are Cas is not responsible for errors in this material as it now exists or will exist, nor does We Are Cas provide technical support for it.
Copyright Info : Gaver and O’Muircheartaigh (1987), Draper (1996)
SAS and all other SAS Institute Inc. product or service names are registered trademarks or trademarks of SAS Institute Inc. in the USA and other countries. ® indicates USA registration. WeAreCAS is an independent community site and is not affiliated with SAS Institute Inc.
This site uses technical and analytical cookies to improve your experience.
Read more.