①lm(Murder ~ Population + Illiteracy + Income + Frost,data=states)函數擬合多元線性回歸 ②在這次實驗中,可以從結果看出,文盲率的回歸系數為4.14(注:4.134e+00=4.134x10^0),所以說,在控制住人口、收入以及氣溫的情況下,文盲率上升1%,那么謀殺率會上升4.14%,它的系數在p<0.001的水平,顯著不為零的。
>#檢測各變量的相關性> states <- as.data.frame(state.x77[,c("Murder","Population","Illiteracy","Income","Frost")])> cor(states)Murder Population Illiteracy Income Frost
Murder 1.00000000.34364280.7029752-0.2300776-0.5388834
Population 0.34364281.00000000.10762240.2082276-0.3321525
Illiteracy 0.70297520.10762241.0000000-0.4370752-0.6719470
Income -0.23007760.2082276-0.43707521.00000000.2262822
Frost -0.5388834-0.3321525-0.67194700.22628221.0000000>#檢測二變量關系> library(car)> scatterplotMatrix(states,spread=FALSE,smooth.args=list(lty=2),main="Scatter Plot Matrix")>#多元線性回歸> fit <- lm(Murder ~ Population + Illiteracy + Income + Frost,data=states)> summary(fit)Call:
lm(formula = Murder ~ Population + Illiteracy + Income + Frost, data = states)Residuals:Min 1Q Median 3Q Max
-4.7960-1.6495-0.08111.48157.6210 Coefficients:Estimate Std. Error t value Pr(>|t|)(Intercept)1.235e+003.866e+000.3190.7510
Population 2.237e-049.052e-052.4710.0173*
Illiteracy 4.143e+008.744e-014.7382.19e-05***
Income 6.442e-056.837e-040.0940.9253
Frost 5.813e-041.005e-020.0580.9541---
Signif. codes:0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error:2.535 on 45 degrees of freedom
Multiple R-squared:0.567, Adjusted R-squared:0.5285
F-statistic:14.73 on 4 and 45 DF, p-value:9.133e-08>#有顯著交互項的多元線性回歸> fit1 <- lm(mpg ~ hp + wt + hp:wt,data = mtcars)> summary(fit1)Call:
lm(formula = mpg ~ hp + wt + hp:wt, data = mtcars)Residuals:Min 1Q Median 3Q Max
-3.0632-1.6491-0.73621.42114.5513 Coefficients:Estimate Std. Error t value Pr(>|t|)(Intercept)49.808423.6051613.8165.01e-14***
hp -0.120100.02470-4.8634.04e-05***
wt -8.216621.26971-6.4715.20e-07***
hp:wt 0.027850.007423.7530.000811***---
Signif. codes:0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error:2.153 on 28 degrees of freedom
Multiple R-squared:0.8848, Adjusted R-squared:0.8724
F-statistic:71.66 on 3 and 28 DF, p-value:2.981e-13> library(statmod)> library(carData)> library(effects)> plot(effect("hp:wt",fit1,,list(wt=c(2.2,3.2,4.2))),multiline=TRUE)
There were 50 or more warnings (use warnings() to see the first 50)