Régression linéaire et prédiction avec lm() et predict()

Question

Je dois appliquer lm() à un sous-ensemble élargi de ma base de données dat, tout en effectuant une prédiction pour l'observation suivante. Par exemple, je fais:

fit model  predict 
----------  ------- 
dat[1:3, ]  dat[4, ] 
dat[1:4, ]  dat[5, ] 
    .    . 
    .    . 
dat[-1, ]  dat[nrow(dat), ] 

Je sais ce que je dois faire pour un sous-ensemble particulier (lié à cette question: predict() and newdata - How does this work?). Par exemple, pour prédire la dernière ligne, je

dat1 = dat[1:(nrow(dat)-1), ] 
dat2 = dat[nrow(dat), ] 

fit = lm(log(clicks) ~ log(v1) + log(v12), data=dat1) 
predict.fit = predict(fit, newdata=dat2, se.fit=TRUE) 

Comment puis-je faire cela automatiquement pour tous les sous-ensembles, et potentiellement extraire ce que je veux dans une table?

• De fit, j'aurais besoin du summary(fit)$adj.r.squared;
• De predict.fit J'aurais besoin de la valeur predict.fit$fit.

Merci.

Answer 1

solution (efficace)

Voici ce que vous pouvez faire:

p <- 3 ## number of parameters in lm() 
n <- nrow(dat) - 1 

## a function to return what you desire for subset dat[1:x, ] 
bundle <- function(x) { 
    fit <- lm(log(clicks) ~ log(v1) + log(v12), data = dat, subset = 1:x, model = FALSE) 
    pred <- predict(fit, newdata = dat[x+1, ], se.fit = TRUE) 
    c(summary(fit)$adj.r.squared, pred$fit, pred$se.fit) 
    } 

## rolling regression/prediction 
result <- t(sapply(p:n, bundle)) 
colnames(result) <- c("adj.r2", "prediction", "se") 

note I ont fait plusieurs choses à l'intérieur de la fonction bundle:

• Je l'ai utilisé l'argument subset pour sélectionner un sous-ensemble pour répondre
• Je l'ai utilisé model = FALSE pour ne pas enregistrer cadre de modèle où nous gagnons l'espace de travail

Dans l'ensemble, il n'y a pas de boucle évidente, mais sapply est utilisé.

• commence Montage de p, le nombre minimum de données nécessaires pour ajuster un modèle avec p coefficients;
• Le montage se termine à nrow(dat) - 1, car nous avons au moins besoin de la dernière colonne pour la prédiction.

---

test

Exemple de données (avec 30 "observations")

dat <- data.frame(clicks = runif(30, 1, 100), v1 = runif(30, 1, 100), 
        v12 = runif(30, 1, 100)) 

Code Application ci-dessus donne results (27 lignes au total, sortie tronquée pour 5 lignes)

  adj.r2 prediction  se 
[1,]   NaN 3.881068  NaN 
[2,] 0.106592619 3.676821 0.7517040 
[3,] 0.545993989 3.892931 0.2758347 
[4,] 0.622612495 3.766101 0.1508270 
[5,] 0.180462206 3.996344 0.2059014 

La première colonne correspond à la valeur ajustée de R.squared pour le modèle ajusté, tandis que la deuxième colonne correspond à la prédiction. La première valeur pour adj.r2 est NaN, parce que le premier modèle que nous ajustons a 3 coefficients pour 3 points de données, donc aucune statistique sensible n'est disponible. La même chose arrive également à se, car la ligne ajustée n'a pas de résidus 0, donc la prédiction se fait sans incertitude.

Answer 2

Je viens de créer des données aléatoires à utiliser pour cet exemple. J'appelle l'objet data parce que c'était ce qu'on appelait dans la question au moment où j'ai écrit cette solution (appelez ça comme vous voulez).

(efficace) Solution

data <- data.frame(v1=rnorm(100),v2=rnorm(100),clicks=rnorm(100)) 

data1 = data[1:(nrow(data)-1), ] 
data2 = data[nrow(data), ] 

for(i in 3:nrow(data)){ 
    nam <- paste("predict", i, sep = "") 
    nam1 <- paste("fit", i, sep = "") 
    nam2 <- paste("summary_fit", i, sep = "") 

    fit = lm(clicks ~ v1 + v2, data=data[1:i,]) 
    tmp <- predict(fit, newdata=data2, se.fit=TRUE) 
    tmp1 <- fit 
    tmp2 <- summary(fit) 
    assign(nam, tmp) 
    assign(nam1, tmp1) 
    assign(nam2, tmp2) 
} 

Tous les résultats que vous voulez seront stockés dans les objets de données cela crée.

Par exemple:

> summary_fit10$r.squared 
[1] 0.3087432 

Vous avez mentionné dans les commentaires que vous souhaitez une table des résultats. Vous pouvez créer programmation tableaux de résultats des 3 types de fichiers de sortie comme ceci:

rm(data,data1,data2,i,nam,nam1,nam2,fit,tmp,tmp1,tmp2) 
frames <- ls() 

frames.fit  <- frames[1:98] #change index or use pattern matching as needed 
frames.predict <- frames[99:196] 
frames.sum  <- frames[197:294] 

fit.table <- data.frame(intercept=NA,v1=NA,v2=NA,sourcedf=NA) 
for(i in 1:length(frames.fit)){ 
    tmp <- get(frames.fit[i]) 
    fit.table    <- rbind(fit.table,c(tmp$coefficients[[1]],tmp$coefficients[[2]],tmp$coefficients[[3]],frames.fit[i])) 
} 

fit.table 

> fit.table 
      intercept     v1     v2 sourcedf 
2 -0.0647017971121678  1.34929652763687 -0.300502017324518 fit10 
3 -0.0401617893034109 -0.034750571912636 -0.0843076273486442 fit100 
4 0.0132968863522573  1.31283604433593 -0.388846211083564 fit11 
5 0.0315113918953643  1.31099122173898 -0.371130010135382 fit12 
6 0.149582794027583 0.958692838785998 -0.299479715938493 fit13 
7 0.00759688947362175 0.703525856001948 -0.297223988673322 fit14 
8 0.219756240025917 0.631961979610744 -0.347851129205841 fit15 
9  0.13389223748979 0.560583832333355 -0.276076134872669 fit16 
10 0.147258022154645 0.581865844000838 -0.278212722024832 fit17 
11 0.0592160359650468 0.469842498721747 -0.163187274356457 fit18 
12 0.120640756525163 0.430051839741539 -0.201725012088506 fit19 
13 0.101443924785995  0.34966728554219 -0.231560038360121 fit20 
14 0.0416637001406594 0.472156988919337 -0.247684504074867 fit21 
15 -0.0158319749710781 0.451944113682333 -0.171367482879835 fit22 
16 -0.0337969739950376 0.423851304105399 -0.157905431162024 fit23 
17 -0.109460218252207  0.32206642419212 -0.055331391802687 fit24 
18 -0.100560410735971 0.335862465403716 -0.0609509815266072 fit25 
19 -0.138175283219818 0.390418411384468 -0.0873106257144312 fit26 
20 -0.106984355317733 0.391270279253722 -0.0560299858019556 fit27 
21 -0.0740684978271464 0.385267011513678 -0.0548056844433894 fit28