GE is dropped, alphabetically, as the intercept term. As eipi10 stated, you can interpret the coefficients for the other levels in states with GE as the baseline (statesLA = 0.1 meaning LA is, on average, 0.1x more than GE). EDIT: To respond to your updated question: If you include all of the levels in a … Read more
I have just realized that there was a closely related question Extracting orthogonal polynomial coefficients from R’s poly() function? 2 years ago. The answer there is merely explaining what predict.poly does, but my answer gives a complete picture. Section 1: How does poly represent orthogonal polynomials My understanding of orthogonal polynomials is that they take … Read more
Provided you do manual weighting correctly, you won’t see discrepancy. So the correct way to go is: X <- model.matrix(~ q + q2 + b + c, mydata) ## non-weighted model matrix (with intercept) w <- mydata$weighting ## weights rw <- sqrt(w) ## root weights y <- mydata$a ## non-weighted response X_tilde <- rw * … Read more
You can inspect the predict function with body(predict.lm). There you will see this line: if (p < ncol(X) && !(missing(newdata) || is.null(newdata))) warning(“prediction from a rank-deficient fit may be misleading”) This warning checks if the rank of your data matrix is at least equal to the number of parameters you want to fit. One way … Read more
Here are some alternatives. No packages are used in the first 3. 1) reformulate fo <- reformulate(regressors, response = “income”) lm(fo, Duncan) or you may wish to write the last line as this so that the formula that is shown in the output looks nicer: do.call(“lm”, list(fo, quote(Duncan))) in which case the Call: line of … Read more
For a single element like this, use [[ rather than [. Compare: coefficients(out)[“newx”] # newx # 1 coefficients(out)[[“newx”]] # [1] 1 More generally, use unname(): unname(coefficients(out)[c(“newx”, “(Intercept)”)]) # [1] 1.0 1.5 head(unname(mtcars)) # NA NA NA NA NA NA NA NA NA NA NA # Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 … Read more
You could subtract the explicit intercept from the regressand and then fit the intercept-free model: > intercept <- 1.0 > fit <- lm(I(x – intercept) ~ 0 + y, lin) > summary(fit) The 0 + suppresses the fitting of the intercept by lm. edit To plot the fit, use > abline(intercept, coef(fit)) P.S. The variables … Read more
I like to use ggplot2 for this because it’s usually very intuitive to add layers of data. library(ggplot2) fit <- lm(mpg ~ hp + I(hp^2), data = mtcars) prd <- data.frame(hp = seq(from = range(mtcars$hp)[1], to = range(mtcars$hp)[2], length.out = 100)) err <- predict(fit, newdata = prd, se.fit = TRUE) prd$lci <- err$fit – 1.96 … Read more
This is not an answer, but is too long for a comment The .Environment for the terms component is identical for each resulting model e1 <- attr(fit[[‘V1’]][[1]]$terms, ‘.Environment’) e2 <- attr(fit[[‘V1’]][[2]]$terms, ‘.Environment’) e3 <- attr(fit[[‘V1’]][[3]]$terms, ‘.Environment’) identical(e1,e2) ## TRUE identical(e2, e3) ## TRUE It appears that data.table is using the same bit of memory (my … Read more
Some statistical result / background (Link in the picture: Function to calculate R2 (R-squared) in R) Computational details Computations involved here is basically the computation of the variance-covariance matrix. Once we have it, results for all pairwise regression is just element-wise matrix arithmetic. The variance-covariance matrix can be obtained by R function cov, but functions … Read more