Menu
I have a regression model for some time series data investigating drug utilisation. The purpose is to fit a spline to a time series and work out 95% CI etc. The model goes as follows:
id <- ts(1:length(drug$Date)) a1 <- ts(drug$Rate) a2 <- lag(a1-1) tg <- ts.union(a1,id,a2) mg <-lm (a1~a2+bs(id,df=df1),data=tg) The summary output of mg is:
Call: lm(formula = a1 ~ a2 + bs(id, df = df1), data = tg) Residuals: Min 1Q Median 3Q Max -0.31617 -0.11711 -0.02897 0.12330 0.40442 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.77443 0.09011 8.594 1.10e-11 *** a2 0.13270 0.13593 0.976 0.33329 bs(id, df = df1)1 -0.16349 0.23431 -0.698 0.48832 bs(id, df = df1)2 0.63013 0.19362 3.254 0.00196 ** bs(id, df = df1)3 0.33859 0.14399 2.351 0.02238 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 I am using the Pr(>|t|) value of a2 to test if the data under investigation are autocorrelated.
Is it possible to extract this value of Pr(>|t|) (in this model 0.33329) and store it in a scalar to perform a logical test?
Alternatively, can it be worked out using another method?
266
Extracting coefficients from glm()The coef() function extracts coefficients and confint() extracts the confidence intervals. These functions work the same for both linear and generalized linear models.
The parameter β (the regression coefficient) signifies the amount by which change in x must be multiplied to give the corresponding average change in y, or the amount y changes for a unit increase in x. In this way it represents the degree to which the line slopes upwards or downwards.
coef is a generic function which extracts model coefficients from objects returned by modeling functions.
A summary.lm object stores these values in a matrix called 'coefficients'. So the value you are after can be accessed with:
a2Pval <- summary(mg)$coefficients[2, 4] Or, more generally/readably, coef(summary(mg))["a2","Pr(>|t|)"]. See here for why this method is preferred.
196
The package broom comes in handy here (it uses the "tidy" format).
tidy(mg) will give a nicely formated data.frame with coefficients, t statistics etc. Works also for other models (e.g. plm, ...).
Example from broom's github repo:
lmfit <- lm(mpg ~ wt, mtcars) require(broom) tidy(lmfit) term estimate std.error statistic p.value 1 (Intercept) 37.285 1.8776 19.858 8.242e-19 2 wt -5.344 0.5591 -9.559 1.294e-10 is.data.frame(tidy(lmfit)) [1] TRUE If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With
