covdesign-20170826

x1<-xcovpart[,1] 
x2<-xcovpart[,2] 
x3<-xcovpart[,3] 
x4<-xcovpart[,4] 
x5<-xcovpart[,5] 
x6<-xcovpart[,6] 
x7<-xcovpart[,7] 
x8<-xcovpart[,8] 
x9<-xcovpart[,9] 
x10<-xcovpart[,10] 
x11<-xcovpart[,11] 
x12<-xcovpart[,12] 
x13<-xcovpart[,13] 
x14<-xcovpart[,14] 
x15<-xcovpart[,15] 
x16<-xcovpart[,16]
x1<-unlist(x1)
x2<-unlist(x2)
x3<-unlist(x3)
x4<-unlist(x4)
x5<-unlist(x5)
x6<-unlist(x6)
x7<-unlist(x7)
x8<-unlist(x8)
x9<-unlist(x9)
x10<-unlist(x10)
x11<-unlist(x11)
x12<-unlist(x12)
x13<-unlist(x13)
x14<-unlist(x14)
x15<-unlist(x15)
x16<-unlist(x16)

library(car)
a1=aov(ypart~0+x1+x2+as.factor(x3)+as.factor(x4)+as.factor(x5)+as.factor(x6)+as.factor(x7)+as.factor(x8)
        +as.factor(x9)+as.factor(x10)+as.factor(x11)+as.factor(x12)+as.factor(x13)+as.factor(x14)
        +as.factor(x15)+as.factor(x16),data=data.frame(xcovpart))
beta=data.frame(coef(a1))#系数
xcovpart=apply(xcovpart,2,as.numeric) #标准化系数
beta_sd=beta*sd(xcovpart)/sd(ypart)
coef(a1)
summary(a1)
aov(a1)
Anova(a1, type="III")
ypart.pre=data.frame(predict(a1))
resi=abs(ypart-ypart.pre)/ypart
apply(resi,2,mean) #误差
yall.pre=data.frame(predict(a1,data.frame(x=xcovall)))
write.table(ypart.pre, file = "/Users/vicky/Documents/code/R/ypartpre.csv")
#bartlett.test(a1,data=xcovpart)

a1.lm = lm(ypart~0+x1+x2+as.factor(x3)+as.factor(x4)+as.factor(x5)+as.factor(x6)+as.factor(x7)+as.factor(x8)
           +as.factor(x9)+as.factor(x10)+as.factor(x11)+as.factor(x12)+as.factor(x13)+as.factor(x14)
           +as.factor(x15)+as.factor(x16),data=data.frame(xcovpart))
summary(a1.lm)
beta1.lm=data.matrix(coef(a1.lm)) #系数
beta1_sta=matrix(0,2,1)
for (i in 1:2){
  beta1_sta[i]=beta1.lm[i]*sd(xcovpart[,i])/sd(ypart)
}

x17<-unlist(x17) #加入房价因素
xcovpart2=cbind(x1,x2,x17,xcovpart[,3:16]) #加入房价因素
a2.lm = lm(ypart~0+x1+x2+x17+as.factor(x3)+as.factor(x4)+as.factor(x5)+as.factor(x6)+as.factor(x7)+as.factor(x8)
           +as.factor(x9)+as.factor(x10)+as.factor(x11)+as.factor(x12)+as.factor(x13)+as.factor(x14)
           +as.factor(x15)+as.factor(x16),data=data.frame(xcovpart2))
summary(a2.lm)
beta2.lm=data.matrix(coef(a2.lm)) #系数
beta2_sta=matrix(0,3)
for (i in 1:3){
 beta2_sta[i]=beta2.lm[i]*sd(xcovpart2[,i])/sd(ypart)
}
summary(a2.lm)
coef1=cor(xcovpart,ypart)
coef2=cor(xcovpart2,ypart)

# xcovpart2=apply(xcovpart2,2,as.numeric)
# xcovpart2_sta=(xcovpart2-apply(xcovpart2,2,mean))/sd(xcovpart)
# a2_sd=aov(ypart~0+xcovpart2_sta[,1]+xcovpart2_sta[,2]+xcovpart2_sta[,3]+xcovpart2_sta[,4]+
#             xcovpart2_sta[,5]+xcovpart2_sta[,6]+as.factor(xcovpart2_sta[,7])+
#             as.factor(xcovpart2_sta[,8])+as.factor(xcovpart2_sta[,9])+as.factor(xcovpart2_sta[,10])+
#             as.factor(xcovpart2_sta[,11])+as.factor(xcovpart2_sta[,12])+as.factor(xcovpart2_sta[,13])+
#             as.factor(xcovpart2_sta[,14])+as.factor(xcovpart2_sta[,15])+as.factor(xcovpart2_sta[,16])+
#             as.factor(xcovpart2_sta[,17]))
# beta2_sd=data.frame(coef(a2_sd))
# summary(a2_sd)

x17all<-unlist(x17all)#加入房价因素
ypart.pre2=data.frame(predict(a2.lm))
resi2=abs(ypart-ypart.pre2)/ypart
apply(resi2,2,mean) #误差
xcovall=data.matrix(xcovall)
xcovall2=cbind(xcovall[,1],xcovall[,2],x17all,xcovall[,3:16])
yall.pre2=predict(a2.lm,newdata=data.frame(xcovall2))
write.table(xcovpart, file = "/Users/vicky/Documents/code/R/xcovpart.csv")
write.table(xcovpart2, file = "/Users/vicky/Documents/code/R/xcovpart2.csv")
write.table(ypart, file = "/Users/vicky/Documents/code/R/ypart.csv")

#分train和test集讨论