setwd('c:/myr')
getwd()
data = read.csv(file='c:/myr/Wage.csv')
age.logit = data$age
wage.logit = vector(length = length(age.logit))
for (i in 1:length(age.logit))
{
if(data$wage[i] > 200)
wage.logit[i] = 1
else
wage.logit[i] = 0
}
#data for logistic model
data.logit = data.frame(age.logit, wage.logit)
###################################################
library(splines)
###################################################
#generate 3 pieces of spline basis , plot
###################################################
bspline = bs(age.logit, df = 3)
data.plot = cbind(age.logit, bspline)
data.plot.sort = data.plot[order(data.plot[, 1]), ]
matplot(age.logit, bspline)
matplot(data.plot.sort[,1], data.plot.sort[, 2:4], type = "l")
###################################################
#logistic cubic spline regression
###################################################
fit = glm(wage.logit ~ bs(age.logit, df = 3), family = binomial(link = "logit"))
summary(fit)
fit=update(fit)
#logit prediction
newdata = seq(20,60)
newdata1 = list(age.logit = newdata)
n<-1000
age<-50+12*rnorm(n)
pred = predict.glm(fit, newdata = list(age.logit = newdata),se.fit = TRUE)
HRpred = Predict(fit,newdata,fun=exp, ref.zero = TRUE) #报错
Error in Predict(fit, newdata, fun = exp, ref.zero = TRUE) :
predictors(s) not in model: newdata
为什么使用Predict就不行呢
1个回答
将fit = glm(wage.logit ~ bs(age.logit, df = 3), family = binomial(link = "logit")) 改成fit = lrm(wage.logit ~ bs(age.logit, df = 3), family = binomial(link = "logit"))就可以了
SofaSofa数据科学社区DS面试题库 DS面经