Published on :
Statistical CREATION_INTERNE

Pump Reliability Analysis using Mixed Models

This code is also available in: Deutsch Español Français
Awaiting validation
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.
Copied!
1DATA pump;
2 pump = _n_;
3 INPUT y t group;
4 logtstd = log(t) - 2.4564900;
5 DATALINES;
6 5 94.320 1
7 1 15.720 2
8 5 62.880 1
9 14 125.760 1
10 3 5.240 2
11 19 31.440 1
12 1 1.048 2
13 1 1.048 2
14 4 2.096 2
15 22 10.480 2
16 ;
2 Code Block
PROC GLIMMIX
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!
1ods html;
2title "Pump Reliability at a Pressurized Water Reactor Nuclear Power Plant";
3PROC GLIMMIX DATA=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;
10RUN;
11ods 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!
1PROC NLMIXED DATA=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);
11RUN;
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)