CODEWITHSUNDEEP
rpredict
B. Davis
B. Davis

Reputation: 3441

manually predict with different functional forms

I have a data frame with coefficients from a glm (betas below). The data frame contains the covariate label, covariate form, and the estimate. The forms are linear (Li), squared/quadratic (Sq), and log (Ps).

betas <- structure(list(CovGen = c("A", "B", "C", "D", "E", "F", "G", 
                                   "G", "H"), Form = c("Li", "Li", "Li", "Li", "Li", "Li", "Li", 
                                                       "Sq", "Ps"), Estimate = c(0.0294573176934061, 0.0100315121169383, 
                                                                                 -0.0155864186367343, -0.00871344935814372, 0.0362538988332902, 
                                                                                 -0.0263072916746069, 0.0865742118052235, 0.0614689145750204, 
                                                                                 0.00229745713752781)), .Names = c("CovGen", "Form", "Estimate"
                                                                                 ), row.names = c(NA, 9L), class = "data.frame")

betas
  CovGen Form     Estimate
1      A   Li  0.029457318
2      B   Li  0.010031512
3      C   Li -0.015586419
4      D   Li -0.008713449
5      E   Li  0.036253899
6      F   Li -0.026307292
7      G   Li  0.086574212
8      G   Sq  0.061468915
9      H   Ps  0.002297457

I trying to apply the coefficient estimates to manually predict values for a new data frame (dat included here with dput). 

dat <- structure(list(B = c(-1.47218074669544, -1.46929972689195, -1.46641870708846, 
                            -1.46353768728497, -1.46065666748148, -1.45777564767799), C = c(-1.09847692593512, 
                                                                                            -1.09375316152745, -1.08902939711978, -1.08430563271211, -1.07958186830444, 
                                                                                            -1.07485810389677), D = c(-1.0109875688763, -1.00407851818141, 
                                                                                                                      -0.997169467486518, -0.990260416791627, -0.983351366096736, -0.976442315401845
                                                                                            ), E = c(-3.19632050296668, -3.19041566990116, -3.18451083683563, 
                                                                                                     -3.17860600377011, -3.17270117070458, -3.16679633763906), F = c(-2.81211918021003, 
                                                                                                                                                                     -2.80673925496675, -2.80135932972346, -2.79597940448018, -2.7905994792369, 
                                                                                                                                                                     -2.78521955399362), G = c(-2.32916817000267, -2.32368219245727, 
                                                                                                                                                                                               -2.31819621491187, -2.31271023736647, -2.30722425982107, -2.30173828227567
                                                                                                                                                                     ), H = c(0.442067970883549, 0.417909464459238, 0.393750958034926, 
                                                                                                                                                                              0.369592451610615, 0.345433945186303, 0.321275438761992)), .Names = c("B", 
                                                                                                                                                                                                                                                    "C", "D", "E", "F", "G", "H"), row.names = c(NA, 6L), class = "data.frame")                                                                                                                                                                                                                                   "C", "D", "E", "F", "G", "H"), row.names = c(NA, 6L), class = "data.frame")



> dat
          B         C          D         E         F         G         H
1 -1.472181 -1.098477 -1.0109876 -3.196321 -2.812119 -2.329168 0.4420680
2 -1.469300 -1.093753 -1.0040785 -3.190416 -2.806739 -2.323682 0.4179095
3 -1.466419 -1.089029 -0.9971695 -3.184511 -2.801359 -2.318196 0.3937510
4 -1.463538 -1.084306 -0.9902604 -3.178606 -2.795979 -2.312710 0.3695925
5 -1.460657 -1.079582 -0.9833514 -3.172701 -2.790599 -2.307224 0.3454339
6 -1.457776 -1.074858 -0.9764423 -3.166796 -2.785220 -2.301738 0.3212754

I am trying to multiply the new data values in the dat df by the respective betas and account for functional forms. More specifically, in the example included here, I want to apply the Sq form of the G beta to dat$G^2 and the Ps H beta to log(dat$H). All other betas and values can simply be multiplied directly without accounting for functional forms. Note that the A beta is not applied to a new value in the dat df.

I may need to work through an award ifelse statement but am wondering if there are other ideas and/or suggestions.

I am working within a larger loop and will not have a consistent form for each covariate.

The desired result would be a matrix or df with a column containing the predicted values for each beta-form combination. For example, there would be a single column containing predicted values for all betas except for G which would have predicted values for G and G^2.

Thanks in advance.

Upvotes: 0

Views: 61

Answers (2)

Gregor Thomas
Gregor Thomas

Reputation: 145965

I'd try to construct the formula and then use model.matrix and matrix multiplication, something like this:

betas$term = with(betas, ifelse(
    Form == "Li", CovGen,
    ifelse(Form == "Sq", sprintf("I(%s^2)", CovGen),
           ifelse(Form == "Ps", sprintf("log(%s)", CovGen), NA)
)))
betas
#   CovGen Form     Estimate   term
# 1      A   Li  0.029457318      A
# 2      B   Li  0.010031512      B
# 3      C   Li -0.015586419      C
# 4      D   Li -0.008713449      D
# 5      E   Li  0.036253899      E
# 6      F   Li -0.026307292      F
# 7      G   Li  0.086574212      G
# 8      G   Sq  0.061468915 I(G^2)
# 9      H   Ps  0.002297457 log(H)
(my_formula = as.formula(paste("~", paste(betas$term, collapse = " + "))))
#~A + B + C + D + E + F + G + I(G^2) + log(H)

X = model.matrix(my_formula, data = dat)
prediction = X %*% betas$Estimate

As MrFlick says, this won't work with your current sample data.

Upvotes: 2

MrFlick
MrFlick

Reputation: 206401

You might try a solution like this

trans <- list(
  Li=identity,
  Sq=function(x) x^2,
  Ps=function(x) log(x)
)

cpredict<-function(betas, datas) {
  Map(function(var, fun, coef) {
    trans[[fun]](datas[[var]])*coef
  }, betas$CovGen, betas$Form, betas$Estimate)
}

cpredict(betas, dat)

But this won't work with your current data because there is no dat$A and you can't take the log of a negative number.

Upvotes: 2

Related Questions