Econometric Analysis
Posted amosding
tags:
篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了Econometric Analysis相关的知识,希望对你有一定的参考价值。
## Econometric Analysis ## install.packages("wooldridge") library(wooldridge); library(glm2); library(mlogit);library(MASS) # 2. Basic method ##### ## 2.1 logistic regression #### data("mroz") ?mroz mod = glm2(inlf ~ nwifeinc + educ + exper + expersq + age + kidslt6 + kidsge6, family=binomial(link = "logit"), data =mroz); summary(mod) coef = summary(mod)$coefficients ## 2.2 multinomial model #### data("Fishing", package = "mlogit") Fish <- mlogit.data(Fishing, varying = c(2:9), shape = "wide", choice = "mode") ## a pure "conditional" model summary(mlogit(mode ~ price + catch, data = Fish)) ## a pure "multinomial model" summary(mlogit(mode ~ 0 | income, data = Fish)) ## which can also be estimated using multinom (package nnet) library("nnet") summary(multinom(mode ~ income, data = Fishing)) ## a "mixed" model m <- mlogit(mode ~ price+ catch | income, data = Fish) summary(m) ## same model with charter as the reference level m <- mlogit(mode ~ price+ catch | income, data = Fish, reflevel = "charter") ## same model with a subset of alternatives : charter, pier, beach m <- mlogit(mode ~ price+ catch | income, data = Fish, alt.subset = c("charter", "pier", "beach")) ## 2.3 poisson regresson data(crabs) satellites <- crabs$Satellites width.shifted <- crabs$Width - min(crabs$Width) dark <- crabs$Dark goodspine <- crabs$GoodSpine fit1 <- glm(satellites ~ width.shifted + factor(dark) + factor(goodspine),family = poisson(link="identity"), start = rep(1,4)) fit2 <- glm2(satellites ~ width.shifted + factor(dark) + factor(goodspine), family = poisson(link="identity"), start = rep(1,4)) fit1.eq <- glm2(satellites ~ width.shifted + factor(dark) + factor(goodspine),family = poisson(link="identity"), start = rep(1,4), method = "glm.fit") ## 2.4 negative regression fit3 <- glm.nb(satellites ~ width.shifted + factor(dark) + factor(goodspine)) coef = summary(fit3)$coefficients ## 2.5 Tobit regression library(survival) # Economists fit a model called `tobit regression‘, which is a standard # linear regression with Gaussian errors, and left censored data. tobinfit <- survreg(Surv(durable, durable>0, type=‘left‘) ~ age + quant,data=tobin, dist=‘gaussian‘) tt = summary(tobinfit) coef = summary(tobinfit)$coefficients ## 2.6 two stage model library(sampleSelection) ## Estimate a simple female wage model taking into account the labour ## force participation data(Mroz87) a <- heckit(lfp ~ huswage + kids5 + mtr + fatheduc + educ + city, log(wage) ~ educ + city, data=Mroz87) ## extract all coefficients of the model: coef( a ) ## now extract the coefficients of the outcome model only: coef( a, part="outcome") ## extract all coefficients, standard errors, t-values ## and p-values of the model: coef( summary( a ) ) ## now extract the coefficients, standard errors, t-values ## and p-values of the outcome model only: coef( summary( a ), part="outcome") ## 2.7 Multilevel Model ## linear mixed models - reference values from older code library(lme4) library(mlmRev) data(Exam);names(Exam) ### null-model mod = lmer(normexam ~ 1 + (1 | school), data=Exam) summary(mod) ### random intercept, fixed predictor in individual level mod = lmer(normexam ~ schavg + (1 | school), data=Exam) summary(mod) ### random intercept, random slope mod = lmer(normexam ~ standLRT + (standLRT | school), data=Exam, method="ML") summary(mod) ### random intercept, individual and group level predictor mod = lmer(normexam ~ standLRT + schavg + (1 + standLRT | school), data=Exam) summary(mod) ### random intercept, cross-level interaction mod = lmer(normexam ~ standLRT * schavg + (1 + standLRT | school), data=Exam) summary(mod) ## 2.7 survival regression # Fit an exponential model: the two fits are the same library(survival) mod = survreg(Surv(futime, fustat) ~ ecog.ps + rx, ovarian, dist=‘weibull‘, scale=1) mod = survreg(Surv(futime, fustat) ~ ecog.ps + rx, ovarian, dist="exponential") # A model with different baseline survival shapes for two groups, i.e., # two different scale parameters data(lung) mod = survreg(Surv(time, status) ~ ph.ecog + age + strata(sex), dist="exponential",data= lung) summary(mod) #### cox model # Create the simplest test data set test1 <- list(time=c(4,3,1,1,2,2,3), status=c(1,1,1,0,1,1,0), x=c(0,2,1,1,1,0,0), sex=c(0,0,0,0,1,1,1)) # Fit a stratified model test2 = coxph(Surv(time, status) ~ x + strata(sex), test1) ### frailty model library(frailtypack) ###--- COX proportional hazard model (SHARED without frailties) ---### ###--- estimated with penalized likelihood ---### data(kidney) frailtyPenal(Surv(time,status)~sex+age, n.knots=12,kappa=10000,data=kidney) ###--- Shared Frailty model ---### frailtyPenal(Surv(time,status)~cluster(id)+sex+age, n.knots=12,kappa=10000,data=kidney) #-- stratified analysis data(readmission) frailtyPenal(Surv(time,event)~cluster(id)+dukes+strata(sex), n.knots=10,kappa=c(10000,10000),data=readmission) ### competing risk library(cmprsk) # simulated data to test set.seed(10) ftime <- rexp(200) fstatus <- sample(0:2,200,replace=TRUE) cov <- matrix(runif(600),nrow=200) dimnames(cov)[[2]] <- c(‘x1‘,‘x2‘,‘x3‘) z <- crr(ftime,fstatus,cov, failcode=2, cencode=0) summary(z) ### mutliple state model library(mstate) data("ebmt4") ebmt = ebmt4 tmat <- transMat(x = list(c(2, 3, 5, 6), c(4, 5, 6), c(4, 5, 6), c(5, 6), c(), c()), names = c("Tx", "Rec", "AE", "Rec+AE", "Rel", "Death")) tmat msebmt <- msprep(data = ebmt, trans = tmat, time = c(NA, "rec", "ae", "recae", "rel", "srv"), status = c(NA, "rec.s", "ae.s", "recae.s", "rel.s", "srv.s"), keep = c("match", "proph", "year", "agecl")) events(msebmt) covs <- c("match", "proph", "year", "agecl") msebmt <- expand.covs(msebmt, covs, longnames = FALSE) msebmt[msebmt$id == 1, -c(9, 10, 12:48, 61:84)] msebmt[, c("Tstart", "Tstop", "time")] <- msebmt[, c("Tstart", "Tstop", "time")]/365.25 c0 <- coxph(Surv(Tstart, Tstop, status) ~ strata(trans), data = msebmt, method = "breslow") msf0 <- msfit(object = c0, vartype = "greenwood", trans = tmat) pt0 <- probtrans(msf0, predt = 0, method = "greenwood") summary(pt0, from = 1) # 4. Pandel Data ## 4.1 # Several models can be estimated with plm by filling the model argument: # the fixed effects model ("within"), # the pooling model ("pooling"), # the first-difference model ("fd"), # the between model ("between"), # the error components model ("random"). library("plm") ## https://cran.r-project.org/web/packages/plm/vignettes/plmPackage.html data("Grunfeld", package="plm") head(Grunfeld) grun.fe <- plm(inv~value+capital, data = Grunfeld, model = "within"); ## fixed effect summary(grun.fe) fixef(grun.fe, type = "dmean") summary(fixef(grun.fe, type = "dmean")) grun.re <- plm(inv~value+capital, data = Grunfeld, model = "random"); ## random effect summary(grun.re) ## Hausman test gw <- plm(inv~value+capital, data=Grunfeld, model="within") gr <- plm(inv~value+capital, data=Grunfeld, model="random") phtest(gw, gr) grun.twfe <- plm(inv~value+capital, data=Grunfeld, model="within", effect="twoways") tt = fixef(grun.twfe, type = "dmean"); ## effect = "time" summary(tt) ### Variable coefficients model data("Gasoline", package = "plm") form <- lgaspcar ~ lincomep + lrpmg + lcarpcap gasw <- plm(form, data = Gasoline, model = "within") gasp <- plm(form, data = Gasoline, model = "pooling") gasnp <- pvcm(form, data = Gasoline, model = "within") pooltest(gasw, gasnp) pooltest(gasp, gasnp) pooltest(form, data = Gasoline, effect = "individual", model = "within") pooltest(form, data = Gasoline, effect = "individual", model = "pooling")
以上是关于Econometric Analysis的主要内容,如果未能解决你的问题,请参考以下文章
MEME(Motif-based sequence analysis tools)使用说明