The plm function computes a piecewise regression model (see Huitema &
McKean, 2000).
Usage
plm(
data,
dvar,
pvar,
mvar,
AR = 0,
model = c("W", "H-M", "B&L-B", "JW"),
family = "gaussian",
trend = TRUE,
level = TRUE,
slope = TRUE,
contrast = c("first", "preceding"),
contrast_level = c(NA, "first", "preceding"),
contrast_slope = c(NA, "first", "preceding"),
formula = NULL,
update = NULL,
na.action = na.omit,
r_squared = TRUE,
var_trials = NULL,
dvar_percentage = FALSE,
...
)Arguments
- data
A single-case data frame. See
scdf()to learn about this format.- dvar
Character string with the name of the dependent variable. Defaults to the attributes in the scdf file.
- pvar
Character string with the name of the phase variable. Defaults to the attributes in the scdf file.
- mvar
Character string with the name of the measurement time variable. Defaults to the attributes in the scdf file.
- AR
Maximal lag of autoregression. Modelled based on the Autoregressive-Moving Average (ARMA) function. When AR is set, the family argument must be set to
family = "gaussian".- model
Model used for calculating the dummy parameters (see Huitema & McKean, 2000). Default is
model = "W". Possible values are:"B&L-B","H-M","W", and deprecated"JW".- family
Set the distribution family. Defaults to a gaussian distribution. See the
familyfunction for more details.- trend
A logical indicating if a trend parameters is included in the model.
- level
A logical indicating if a level parameters is included in the model.
- slope
A logical indicating if a slope parameters is included in the model.
- contrast
Sets contrast_level and contrast_slope. Either "first", "preceding" or a contrast matrix.
- contrast_level
Either "first", "preceding" or a contrast matrix. If NA contrast_level is a copy of contrast.
- contrast_slope
Either "first", "preceding" or a contrast matrix. If NA contrast_level is a copy of contrast.
- formula
Defaults to the standard piecewise regression model. The parameter phase followed by the phase name (e.g., phaseB) indicates the level effect of the corresponding phase. The parameter 'inter' followed by the phase name (e.g., interB) adresses the slope effect based on the method provide in the model argument (e.g., "B&L-B"). The formula can be changed for example to include further variables into the regression model.
- update
An easier way to change the regression formula (e.g.,
. ~ . + newvariable).- na.action
Defines how to deal with missing values.
- r_squared
Logical. If TRUE, delta r_squares will be calculated for each predictor.
- var_trials
Name of the variable containing the number of trials (only for binomial regressions). If a single integer is provided this is considered to be a the constant number of trials across all measurements.
- dvar_percentage
Only for binomial distribution. If set TRUE, the dependent variable is assumed to represent proportions
[0,1]. Otherwise dvar is assumed to represent counts.- ...
Further arguments passed to the glm function.
Value
- formula
plm formula. Uselful if you want to use the update or formula argument and you don't know the names of the parameters.
- model
Character string from function call (see
Argumentsabove).- F.test
F-test values of modelfit.
- r.squares
Explained variance R squared for each model parameter.
- ar
Autoregression lag from function call (see
Argumentsabove).- family
Distribution family from function call (see
Argumentsabove).- full.model
Full regression model list from the gls or glm function.
References
Beretvas, S., & Chung, H. (2008). An evaluation of modified R2-change effect size indices for single-subject experimental designs. Evidence-Based Communication Assessment and Intervention, 2, 120-128.
Huitema, B. E., & McKean, J. W. (2000). Design specification issues in time-series intervention models. Educational and Psychological Measurement, 60, 38-58.
See also
Other regression functions:
autocorr(),
corrected_tau(),
hplm(),
mplm(),
trend()
Examples
## Compute a piecewise regression model for a random single-case
set.seed(123)
AB <- design(
phase_design = list(A = 10, B = 20),
level = list(A = 0, B = 1), slope = list(A = 0, B = 0.05),
trend = 0.05
)
dat <- random_scdf(design = AB)
plm(dat, AR = 3)
#> Piecewise Regression Analysis
#>
#> Contrast model: W / level = first, slope = first
#>
#> Fitted a gaussian distribution.
#> Correlated residuals up to autoregressions of lag 3 are modelled
#>
#> F(3, 26) = 47.38; p = 0.000; R² = 0.845; Adjusted R² = 0.828
#>
#> B 2.5% 97.5% SE t p delta R²
#> Intercept 53.068 48.188 57.948 2.490 21.315 0.000
#> Trend mt -0.166 -1.032 0.701 0.442 -0.374 0.711 -0.001
#> Level phase B 16.266 10.061 22.470 3.166 5.138 0.000 0.069
#> Slope phase B 0.913 -0.045 1.871 0.489 1.869 0.073 0.016
#>
#> Autocorrelations of the residuals
#> lag cr
#> 1 -0.12
#> 2 -0.23
#> 3 0.38
#> Formula: values ~ 1 + mt + phaseB + interB
#>
## Another example with a more complex design
A1B1A2B2 <- design(
phase_design = list(A1 = 15, B1 = 20, A2 = 15, B2 = 20),
level = list(A1 = 0, B1 = 1, A2 = -1, B2 = 1),
slope = list(A1 = 0, B1 = 0.0, A2 = 0, B2 = 0.0),
trend = 0.0)
dat <- random_scdf(design = A1B1A2B2, seed = 123)
plm(dat, contrast = "preceding")
#> Piecewise Regression Analysis
#>
#> Contrast model: W / level = preceding, slope = preceding
#>
#> Fitted a gaussian distribution.
#> F(7, 62) = 11.05; p = 0.000; R² = 0.555; Adjusted R² = 0.505
#>
#> B 2.5% 97.5% SE t p delta R²
#> Intercept 51.354 46.820 55.889 2.314 22.196 0.000
#> Trend mt -0.085 -0.636 0.467 0.281 -0.301 0.764 0.001
#> Level phase B1 7.816 1.418 14.213 3.264 2.395 0.020 0.041
#> Level phase A2 -11.360 -17.599 -5.121 3.183 -3.569 0.001 0.091
#> Level phase B2 10.380 3.983 16.777 3.264 3.180 0.002 0.073
#> Slope phase B1 0.280 -0.377 0.937 0.335 0.835 0.407 0.005
#> Slope phase A2 -0.240 -0.897 0.417 0.335 -0.716 0.477 0.004
#> Slope phase B2 0.119 -0.538 0.776 0.335 0.354 0.724 0.001
#>
#> Autocorrelations of the residuals
#> lag cr
#> 1 -0.05
#> 2 -0.12
#> 3 0.12
#> Ljung-Box test: X²(3) = 2.27; p = 0.519
#>
#> Formula: values ~ 1 + mt + phaseB1 + phaseA2 + phaseB2 + interB1 + interA2 +
#> interB2
#>
## no slope effects were found. Therefore, you might want to the drop slope
## estimation:
plm(dat, slope = FALSE, contrast = "preceding")
#> Piecewise Regression Analysis
#>
#> Contrast model: W / level = preceding, slope = preceding
#>
#> Fitted a gaussian distribution.
#> F(4, 65) = 19.73; p = 0.000; R² = 0.548; Adjusted R² = 0.521
#>
#> B 2.5% 97.5% SE t p delta R²
#> Intercept 50.232 47.470 52.994 1.409 35.645 0.000
#> Trend mt 0.076 -0.133 0.284 0.107 0.710 0.480 0.004
#> Level phase B1 7.671 2.879 12.463 2.445 3.137 0.003 0.068
#> Level phase A2 -10.945 -15.737 -6.153 2.445 -4.477 0.000 0.139
#> Level phase B2 9.402 4.610 14.194 2.445 3.846 0.000 0.103
#>
#> Autocorrelations of the residuals
#> lag cr
#> 1 -0.03
#> 2 -0.09
#> 3 0.13
#> Ljung-Box test: X²(3) = 2.02; p = 0.568
#>
#> Formula: values ~ 1 + mt + phaseB1 + phaseA2 + phaseB2
#>
## and now drop the trend estimation as well
plm(dat, slope = FALSE, trend = FALSE, contrast = "preceding")
#> Piecewise Regression Analysis
#>
#> Contrast model: W / level = preceding, slope = preceding
#>
#> Fitted a gaussian distribution.
#> F(3, 66) = 26.34; p = 0.000; R² = 0.545; Adjusted R² = 0.524
#>
#> B 2.5% 97.5% SE t p delta R²
#> Intercept 50.762 48.427 53.097 1.191 42.611 0
#> Level phase B1 8.995 5.906 12.084 1.576 5.708 0 0.225
#> Level phase A2 -9.621 -12.710 -6.532 1.576 -6.105 0 0.257
#> Level phase B2 10.726 7.637 13.815 1.576 6.806 0 0.319
#>
#> Autocorrelations of the residuals
#> lag cr
#> 1 -0.01
#> 2 -0.08
#> 3 0.13
#> Ljung-Box test: X²(3) = 1.82; p = 0.611
#>
#> Formula: values ~ 1 + phaseB1 + phaseA2 + phaseB2
#>
## A poisson regression
example_A24 |>
plm(family = "poisson")
#> Piecewise Regression Analysis
#>
#> Contrast model: W / level = first, slope = first
#>
#> Fitted a poisson distribution.
#> X²(3) = 547.67; p = 0.000; AIC = 261
#>
#> B 2.5% 97.5% SE t p Odds Ratio 2.5%
#> Intercept 5.556 5.472 5.638 0.042 131.693 0.000 258.786 237.936
#> Trend year 0.007 -0.016 0.030 0.012 0.604 0.546 1.007 0.984
#> Level phase B -0.806 -0.938 -0.674 0.067 -11.964 0.000 0.447 0.391
#> Slope phase B -0.006 -0.031 0.019 0.013 -0.472 0.637 0.994 0.969
#> 97.5% Yule's Q 2.5% 97.5%
#> Intercept 280.900 0.99 0.99 0.99
#> Trend year 1.030 0.00 -0.01 0.01
#> Level phase B 0.510 -0.38 -0.44 -0.32
#> Slope phase B 1.019 0.00 -0.02 0.01
#>
#> Formula: injuries ~ 1 + year + phaseB + interB
#>
## A binomial regression (frequencies as dependent variable)
plm(exampleAB_score$Christiano, family = "binomial", var_trials = "trials")
#> Piecewise Regression Analysis
#>
#> Contrast model: W / level = first, slope = first
#>
#> Fitted a binomial distribution.
#> X²(3) = 240.66; p = 0.000; AIC = 120
#>
#> B 2.5% 97.5% SE t p Odds Ratio 2.5% 97.5%
#> Intercept -1.964 -2.793 -1.239 0.394 -4.991 0.000 0.140 0.061 0.290
#> Trend mt 0.023 -0.118 0.166 0.072 0.324 0.746 1.023 0.889 1.181
#> Level phase B 2.376 1.454 3.378 0.488 4.866 0.000 10.762 4.280 29.312
#> Slope phase B 0.038 -0.111 0.186 0.075 0.504 0.614 1.039 0.895 1.204
#> Yule's Q 2.5% 97.5%
#> Intercept -0.75 -0.88 -0.55
#> Trend mt 0.01 -0.06 0.08
#> Level phase B 0.83 0.62 0.93
#> Slope phase B 0.02 -0.06 0.09
#>
#> Formula: values/trials ~ 1 + mt + phaseB + interB
#> weights = trials
## A binomial regression (percentage as dependent variable)
exampleAB_score$Christiano |>
transform(percentage = values/trials) |>
set_dvar("percentage") |>
plm(family = "binomial", var_trials = "trials", dvar_percentage = TRUE)
#> Piecewise Regression Analysis
#>
#> Contrast model: W / level = first, slope = first
#>
#> Fitted a binomial distribution.
#> X²(3) = 240.66; p = 0.000; AIC = 120
#>
#> B 2.5% 97.5% SE t p Odds Ratio 2.5% 97.5%
#> Intercept -1.964 -2.793 -1.239 0.394 -4.991 0.000 0.140 0.061 0.290
#> Trend mt 0.023 -0.118 0.166 0.072 0.324 0.746 1.023 0.889 1.181
#> Level phase B 2.376 1.454 3.378 0.488 4.866 0.000 10.762 4.280 29.312
#> Slope phase B 0.038 -0.111 0.186 0.075 0.504 0.614 1.039 0.895 1.204
#> Yule's Q 2.5% 97.5%
#> Intercept -0.75 -0.88 -0.55
#> Trend mt 0.01 -0.06 0.08
#> Level phase B 0.83 0.62 0.93
#> Slope phase B 0.02 -0.06 0.09
#>
#> Formula: percentage ~ 1 + mt + phaseB + interB
#> weights = trials