2011-07-05 20 views
45

Tengo un modelo de regresión para algunos datos de series de tiempo que investigan la utilización de fármacos. El objetivo es ajustar una spline a una serie de tiempo y trabajar en 95% IC, etc. El modelo es el siguiente:Extraer los valores del coeficiente de regresión

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) 

La salida resumen de mg es:

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 

estoy usando el Pr(>|t|) valor de a2 para comprobar si los datos bajo investigación están autocorrelacionados.

¿Es posible extraer este valor de Pr(>|t|) (en este modelo 0.33329) y almacenarlo en un escalar para realizar una prueba lógica?

Alternativamente, ¿se puede resolver usando otro método?

+0

. @ John - ¿Por qué usaste 'Pr (> | t |)' valor de 'a2' y ninguna de las tres primeras columnas? –

Respuesta

53

Un objeto summary.lm almacena estos valores en un matrix llamado 'coefficients'. Por lo que el valor que está después se puede acceder con:

a2Pval <- summary(mg)$coefficients[2, 4] 

O, en términos más generales/legible, coef(summary(mg))["a2","Pr(>|t|)"]. Consulte here para saber por qué se prefiere este método.

18

El paquete broom es útil aquí (usa el formato "ordenado").

tidy(mg) dará un formato de datos muy formados con coeficientes, t estadísticas, etc. Funciona también para otros modelos (por ejemplo, plm, ...).

Ejemplo de broom 's repo github:

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 
+0

Para responder al OP a partir de esto: '' td [1, "estimate"] '' o '' td [td $ term == "(Interceptar)", "estimar"] '' o incluso '' tdt <- as .data.table (td); setkey (tdt); tdt ["(Interceptar)", "estimar"] '' – PatrickT

0

sólo tiene que pasar su modelo de regresión en la siguiente función:

plot_coeffs <- function(mlr_model) { 
     coeffs <- coefficients(mlr_model) 
     mp <- barplot(coeffs, col="#3F97D0", xaxt='n', main="Regression Coefficients") 
     lablist <- names(coeffs) 
     text(mp, par("usr")[3], labels = lablist, srt = 45, adj = c(1.1,1.1), xpd = TRUE, cex=0.6) 
    } 

Uso como sigue:

model <- lm(Petal.Width ~ ., data = iris) 

plot_coeffs(model) 

enter image description here

Cuestiones relacionadas