代码改变世界

longitudinal model

2014-07-26 11:40  breece  阅读(340)  评论(0编辑  收藏  举报

1.读数据
    excel文件
        1.library(RODBC);z<-odbcConnectExcel("weather.xls");(w<-sqlFetch(z,"Sheet1"));odbcClose(z);
        选择文件:    data=odbcConnectExcel(file.choose());显示文件信息:sqlTables(data);
        2.文件另存为prn格式;w<-read.table("weather.prn",header=T)
        3.文件另存为csv格式:w<-read.csv("weather.prn",header=T)

    文本文件
        x=read.table("weather.txt");编辑文档edit(x);


2.写文件
    write.table(w,file="e:\\rtemp\\weather.txt")

   
3.longitudinal model
    建立模型
        mixed effects model:
            library(nlme)
            e.g:lme . fit1 <- lme ( distance ~age , data = Orthodont ,random = ~ age | Subject , method = "ML")

            lme.fit1<-lme(Pm~Temperature,data=w,random=~Temperature|Location,method="ML");summary(lme.fit1);
            lme.fit2<-update(lme.fit1,fixed=Pm~Temperature*Time);summary(lme.fit2);

                lme.fit1=lme(logpm~Wind+Population+Green+Weather+Road+GDP+Temperature+Zone,data=w.dat,random=~GDP|Location,method="ML")


            模型比较
                anova(lme.fit1,lme.fit2)
       
        Gee model:
            fit.gee1=gee(Pm~Zone,Location,data=w.dat);
            fit.gee2=gee(Pm~Zone,Location,data=w.dat,corstr="AR-M")
            fit.gee3=gee(Pm~Zone,Location,data=w.dat,corstr="exchangeable")
            fit.gee4=gee(Pm~Zone,Location,data=w.dat,corstr="unstructured")
           
            fit.gee1=gee(Pm~Wind+Population+Green+Weather+Road+GDP+Temperature+Zone,Location,data=w.dat)

 

4.画图
    散点图
        attach(w);plot(Temperature,Pm);plot(z$Time,z$Pm,type="l")

        plot(z$Time,z$Pm,type="b");qqnorm(Pm)
        plot(w.dat$Time,w.dat$Pm,type="b")
5.退出
    q()
    清屏:Ctrl+L

6.分类数据
    Weather=1:3;Weather=factor(Weather);Location=1:20;Location=factor(Location); Time=1:7;Time=factor(Time);
Zone=1:4;Zone=factor(Zone);


lpm=log(Pm)
w.dat=groupedData(Pm~GDP|Location,data=w);

 

7.对残差进行分析:   
    plot(lme.fit3,resid(.,type="p")~fitted(.)|Zone,id=0.05,adj=-0.3)
    qqnorm(lme.fit3,~resid(.)|Zone)