Hi everyone,
I tried to replicate Tina's model from the last post and I got the
following error:
library("Amelia")
library("Zelig")
data(freetrade)
a.out <- amelia(freetrade, m=5, ts="year", cs="country")
z.out <- zelig(polity~tariff+gdp.pc, model="ls", data=a.out$imputations)
summary(z.out)
Error in is.null(subset) : 'subset' is missing
Am I doing something wrong? Is there any other way to combine imputed data sets?
Thanks in advance for your attention,
DF
On 23 May 2013 18:00, <zelig-request(a)lists.gking.harvard.edu> wrote:
> Send Zelig mailing list submissions to
> zelig(a)lists.gking.harvard.edu
>
> To subscribe or unsubscribe via the World Wide Web, visit
> https://lists.gking.harvard.edu/mailman/listinfo/zelig
> or, via email, send a message with subject or body 'help' to
> zelig-request(a)lists.gking.harvard.edu
>
> You can reach the person managing the list at
> zelig-owner(a)lists.gking.harvard.edu
>
> When replying, please edit your Subject line so it is more specific
> than "Re: Contents of Zelig digest..."
>
>
> Today's Topics:
>
> 1. How to get measures of model fit (AIC, F-statistics) in zelig
> for multiply imputed data? (Tina)
>
>
> ----------------------------------------------------------------------
>
> Message: 1
> Date: Wed, 22 May 2013 16:05:13 +0000 (UTC)
> From: Tina <freyburg(a)eup.gess.ethz.ch>
> To: zelig(a)lists.gking.harvard.edu
> Subject: [zelig] How to get measures of model fit (AIC, F-statistics)
> in zelig for multiply imputed data?
> Message-ID: <loom.20130522T180408-181(a)post.gmane.org>
> Content-Type: text/plain; charset=us-ascii
>
> I am interested in learning how to get the usual measures of the relative
> quality of a statistical model in zelig for regression using multiply
> imputed data (created with Amelia).
>
>
>
> require(Zelig)
>
> require(Amelia)
>
> data(freetrade)
>
>
>
> #Imputation of missing data
>
> a.out <- amelia(freetrade, m=5, ts="year", cs="country")
>
>
>
> # Regression model
>
> z.out <- zelig(polity~tariff+gdp.pc, model="ls", data=a.out$imputations)
>
>
>
> summary(z.out)
>
>
>
> Model: ls
>
> Number of multiply imputed data sets: 5
>
> Combined results:
>
> Call:
>
> lm(formula = formula, weights = weights, model = F, data = data)
>
> Coefficients:
>
> Value Std. Error t-stat p-value
>
> (Intercept) 1.6740501340 1.0270535468 1.6299541 0.10342186
>
> tariff 0.0196015092 0.0233789523 0.8384255 0.40234214
>
> gdp.pc 0.0003296261 0.0001844909 1.7866798 0.07409327
>
> For combined results from datasets i to j, use summary(x, subset = i:j).
>
> For separate results, use print(summary(x), subset = i:j).
>
>
>
> Question: Does anyone know how to get the values of AIC, F-statistics and
> the degree of freedoms? That would be fantastic!
>
>
>
>
>
> ------------------------------
>
> _______________________________________________
> Zelig mailing list
> Zelig(a)lists.gking.harvard.edu
>
> To unsubscribe from this list or get other information:
>
> https://lists.gking.harvard.edu/mailman/listinfo/zelig
>
> End of Zelig Digest, Vol 108, Issue 1
> *************************************
>
I am interested in learning how to get the usual measures of the relative
quality of a statistical model in zelig for regression using multiply
imputed data (created with Amelia).
require(Zelig)
require(Amelia)
data(freetrade)
#Imputation of missing data
a.out <- amelia(freetrade, m=5, ts="year", cs="country")
# Regression model
z.out <- zelig(polity~tariff+gdp.pc, model="ls", data=a.out$imputations)
summary(z.out)
Model: ls
Number of multiply imputed data sets: 5
Combined results:
Call:
lm(formula = formula, weights = weights, model = F, data = data)
Coefficients:
Value Std. Error t-stat p-value
(Intercept) 1.6740501340 1.0270535468 1.6299541 0.10342186
tariff 0.0196015092 0.0233789523 0.8384255 0.40234214
gdp.pc 0.0003296261 0.0001844909 1.7866798 0.07409327
For combined results from datasets i to j, use summary(x, subset = i:j).
For separate results, use print(summary(x), subset = i:j).
Question: Does anyone know how to get the values of AIC, F-statistics and
the degree of freedoms? That would be fantastic!
Hi Matt,
I saw this on the lists as a problem in January but I don't know whether a
fix was found.
The error is replicated using the sleepstudy data included in lme4.
z.out<-zelig(Reaction ~ Days + tag(Days|Subject), data=sleepstudy,
model="ls.mixed")
x.out<-setx(z.out, Days=4.4)
s.out<-sim(z.out, x=x.out)
Fails with:
Error in array(x, c(length(x), 1L), if (!is.null(names(x))) list(names(x),
:
attempt to set an attribute on NULL
Best wishes,
Ewen
Hi,
I was just running an old script, and some Zelig commands that used to work are not working anymore (in Zelig 4.1-3).
The problem is with including a logical vector as an x-variable. The zelig() function works properly (with the logical vector treated as a 0:1 dummy). Things break down in setx(), which doesn't seem to convert the logical values into their numerical equivalents, resulting in a "model matrix" that is NA (and a "non-comformable arguments" error when fed to sim()).
A workaround is to enter the numerical equivalents of TRUE and FALSE in setx(), but it would be more intuitive to be able to enter the logical values directly, as in the past.
(Minimal working example below.)
Rod
## Minimal working example -- Zelig problem with logical x
require(Zelig)
df <- data.frame(y=rnorm(20),
xLogical=sample(c(TRUE, FALSE),
size=20,
replace=TRUE))
z.out <- zelig(y ~ xLogical, data=df, model="ls")
x.out <- setx(z.out, xLogical=FALSE)
summary(x.out) # model matrix is NA
s.out <- sim(z.out, x.out) # error: non-conformable arguments
## Workaround
x.out2 <- setx(z.out, xLogical=0)
summary(x.out2) # model matrix is OK
s.out <- sim(z.out, x.out2)
summary(s.out)
<table width="100%" border="0" cellspacing="0" cellpadding="0" style="width:100%;">
<tr>
<td align="left" style="text-align:justify;"><font face="arial,sans-serif" size="1" color="#999999"><span style="font-size:11px;">This communication is intended for the addressee only. It is confidential. If you have received this communication in error, please notify us immediately and destroy the original message. You may not copy or disseminate this communication without the permission of the University. Only authorised signatories are competent to enter into agreements on behalf of the University and recipients are thus advised that the content of this message may not be legally binding on the University and may contain the personal views and opinions of the author, which are not necessarily the views and opinions of The University of the Witwatersrand, Johannesburg. All agreements between the University and outsiders are subject to South African Law unless the University agrees in writing to the contrary. </span></font></td>
</tr>
</table
Hi Zelig list,
I'm trying to get Zelig up and running to to perform some
Clarify-like magic on my negative binomial regression. Unfortunately,
while the model runs and returns a model object, it also gives the
following warning:
In z(.function = "glm.nb", .hook = "robust.glm.hook", weights = weights, :
.hook parameter must be a function. ignoring.
In addition to making me uneasy, there appears to be no difference
between the robust and non-robust versions of the model.
About my computer:
R 2.15.2
Zelig beta (installed today, after the released version of Zelig also
installed today returned the same warning)
Below is an example that reproduces the problem seen in my real data.
Is this just an issue with how my particular copy of R is configured,
or is there something deeper going on? If the former, how do I fix
it?
Any tips you can provide would be most welcome.
Thanks,
Ari
data(turnout)
dat <- turnout
dat$count <- rpois( n=nrow(dat), lambda=ceil(dat$age/15) )
Zunrobust <- zelig(count ~ age, model = "negbinom", data = dat)
Zrobust <- zelig(count ~ age, model = "negbinom", data = dat, robust=TRUE)
> summary(Zrobust)
Call:
"glm.nb"(weights = weights, formula = formula, data = data, .hook =
"robust.glm.hook")
Deviance Residuals:
Min 1Q Median 3Q Max
-3.1882 -0.8442 -0.0998 0.5723 3.6247
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.3582217 0.0360647 9.933 <2e-16 ***
age 0.0184236 0.0006645 27.725 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for Negative Binomial(115.4597) family taken to be 1)
Null deviance: 2917.6 on 1999 degrees of freedom
Residual deviance: 2164.1 on 1998 degrees of freedom
AIC: 7954.6
Number of Fisher Scoring iterations: 1
Theta: 115
Std. Err.: 116
2 x log-likelihood: -7948.563
> summary(Zunrobust)
Call:
"glm.nb"(weights = weights, formula = formula, data = data, .hook =
"robust.glm.hook")
Deviance Residuals:
Min 1Q Median 3Q Max
-3.1882 -0.8442 -0.0998 0.5723 3.6247
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.3582217 0.0360647 9.933 <2e-16 ***
age 0.0184236 0.0006645 27.725 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for Negative Binomial(115.4597) family taken to be 1)
Null deviance: 2917.6 on 1999 degrees of freedom
Residual deviance: 2164.1 on 1998 degrees of freedom
AIC: 7954.6
Number of Fisher Scoring iterations: 1
Theta: 115
Std. Err.: 116
2 x log-likelihood: -7948.563
Dear Zelig Community,
I have just updated my R from version 2.13 to 3.0.0, and along with it,
my Zelig package. I have updated my scripts according to the newest
modifications in the Zelig and ZeligChoice package, but I stil bump into
a problem when I try to use the setx function within an ordinary
function command (even though it DID work under 2.13!).
My problem:
While the command lines zelig() and setx() work perfectly in a normal
environment, once I embed the zelig() and setx() into a function() [I
need to embed the zelig functions due to further analysis], the whole
thing collapses. I demonstrate on the example given in ZeligChoice (see
below), but I get the same error for my own data.
I would appreciate your help with this problem - it would be very
important for me to get the setx and the sim into the function, so I can
then extract the $qi$ev for some further analysis.
Thanks for your response,
Zsófia Ignácz
####This works
data(mexico)
z.out <- zelig(as.factor(vote88) ~ pristr + othcok + othsocok, model =
"mlogit",
data = mexico)
x.weak <- setx(z.out, pristr = 1, othcok=2, othsocok=3) ###I chose the
values without considering any substantial meaning
###This doesn't (even though technically identical)
dem.func <- function (a,b,c) {
z.out <- zelig(as.factor(vote88) ~ pristr + othcok + othsocok, model =
"mlogit",
data = mexico)
x.weak <- setx(z.out, pristr=a,othcok=b,othsocok=c)
}
dem.func (a=1,b=2,c=3)
###Error:
Error in eval(expr, envir, enclos) : object 'a' not found
###Technical Info:
R version 3.0.0 (2013-04-03) -- "Masked Marvel"
Copyright (C) 2013 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)
ZELIG (Versions 4.1-3, built: 2013-01-30)
--
Zsófia Ignácz
PhD Student
Humboldt Universität zu Berlin
Berlin Graduate School of Social Sciences European PhD for
Social-economic
and Social Statistics
Email: zsofia.ignacz(a)sowi.hu-berlin.de
Research Associate
Jacobs University Bremen
Campus Ring 1
D-28759 Bremen
Germany
Phone: +49 421 200-3012
Email: z.ignacz(a)jacobs-university.de
URL: http://jacobs-university.de/directory/zignacz
Dear Zelig listserv,
I would like to implement splines along one independent variable in a logistic regression model with multiples IVs. I would like to establish knots at 3 different levels of one IV (e.g., low, medium, high).
I have searched online and reviewed the Zelig reference manual and cannot find any documentation for such a procedure. I see the logit.gam procedure (12.24) incorporates knots but there is not much
description in the manual. An example is provided from mgcv gam is:
RRR> z.out <-zelig(y~s(x0,k=4,fx=TRUE,bs="tp")+s(x1,k=12)+s(x2,k=15),
model="logit.gam", data=my.data)
I've looked in the mgcv and gam sections and not found any examples that help with what I am trying to do.
Questions:
1. Does the logit.gam procedure allow an equivalent to a spline regression (e.g., mkspline in stata: http://www.stata.com/help.cgi?mkspline)? I've pasted stata code below. Or is there another Zelig procedure you'd recommend?
2. I see k= command and a knots = in the documentation. Are there any other examples of how to use these or where to read more about their use?
Much thanks for your time,
David
## stata code spline and logistic regression predicting enrollment
mkspline nw1 -8.3938 nw2 8.294 nw3 = ASSETS_LIAB_NWIHS, marginal displayknots
logistic Enroll AgeEnroll Female Married HHSize Eduy Eft Ept Estu TANFEver Home Stock CheckA SavA nw1-nw3
Hi,
I want to combine multiple amelia runs, and use the zelig function on the entire set of amelia outputs. However, I noticed that when I run zelig on the combined set, the coefficients don't look right. Here are the script and the R outputs:
> require(Amelia)
## Amelia II: Multiple Imputation
## (Version 1.7, built: 2013-02-10)
> library(Zelig)
Attaching package: ‘zoo’
The following object(s) are masked from ‘package:base’:
as.Date, as.Date.numeric
ZELIG (Versions 4.1-3, built: 2013-01-30)
> data(freetrade)
> a.out1 <-amelia(freetrade, m = 1, ts = "year", cs = "country")
-- Imputation 1 --
1 2 3 4 5 6 7 8 9 10 11 12 13 14
> a.out2 <-amelia(freetrade, m = 1, ts = "year", cs = "country")
-- Imputation 1 --
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
> a.out3 <-ameliabind(a.out1, a.out2) # combine the amelia runs
> summary(a.out3)
Amelia output with 2 imputed datasets.
Return code: 1
Message: Normal EM convergence
Chain Lengths:
--------------
Imputation 1: 14
Imputation 2: 15
Rows after Listwise Deletion: 96
Rows after Imputation: 171
Patterns of missingness in the data: 8
Fraction Missing for original variables:
-----------------------------------------
Fraction Missing
year 0.00000000
country 0.00000000
tariff 0.33918129
polity 0.01169591
pop 0.00000000
gdp.pc 0.00000000
intresmi 0.07602339
signed 0.01754386
fiveop 0.10526316
usheg 0.00000000
# run zelig on the first amelia run
> z.out1 <- zelig(tariff ~ polity + pop + gdp.pc, data = a.out1$imputations, model = "ls", cite = FALSE)
> summary(z.out1)
Call:
lm(formula = formula, weights = weights, model = F, data = data)
Residuals:
Min 1Q Median 3Q Max
-29.686 -11.067 -3.367 10.003 45.556
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.317e+01 1.891e+00 17.538 < 2e-16 ***
polity -1.060e-01 2.434e-01 -0.435 0.664
pop 2.925e-08 5.405e-09 5.411 2.15e-07 ***
gdp.pc -2.738e-03 5.235e-04 -5.231 4.99e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 16.62 on 167 degrees of freedom
Multiple R-squared: 0.3234, Adjusted R-squared: 0.3112
F-statistic: 26.6 on 3 and 167 DF, p-value: 4.07e-14
# run zelig on the second amelia run
> z.out2 <- zelig(tariff ~ polity + pop + gdp.pc, data = a.out2$imputations, model = "ls", cite = FALSE)
> summary(z.out2)
Call:
lm(formula = formula, weights = weights, model = F, data = data)
Residuals:
Min 1Q Median 3Q Max
-40.13 -11.72 -3.63 10.16 43.77
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.400e+01 1.984e+00 17.133 < 2e-16 ***
polity -6.123e-02 2.545e-01 -0.241 0.81
pop 3.289e-08 5.689e-09 5.781 3.56e-08 ***
gdp.pc -3.016e-03 5.501e-04 -5.484 1.51e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 17.47 on 167 degrees of freedom
Multiple R-squared: 0.3495, Adjusted R-squared: 0.3378
F-statistic: 29.91 on 3 and 167 DF, p-value: 1.574e-15
# run zelig on the combined amelia runs
> z.out3 <- zelig(tariff ~ polity + pop + gdp.pc, data = a.out3$imputations, model = "ls", cite = FALSE)
> summary(z.out3)
Call:
lm(formula = formula, weights = weights, model = F, data = data)
Residuals:
Min 1Q Median 3Q Max
-29.686 -11.067 -3.367 10.003 45.556
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.317e+01 1.891e+00 17.538 < 2e-16 ***
polity -1.060e-01 2.434e-01 -0.435 0.664
pop 2.925e-08 5.405e-09 5.411 2.15e-07 ***
gdp.pc -2.738e-03 5.235e-04 -5.231 4.99e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 16.62 on 167 degrees of freedom
Multiple R-squared: 0.3234, Adjusted R-squared: 0.3112
F-statistic: 26.6 on 3 and 167 DF, p-value: 4.07e-14
Can someone please tell me why summary(z.out3) is returning the exact same results as summary(z.out1)?
How can I fix this problem so that the ameliabind output is treated the same as a direct output from the amelia function?
Thanks.
Dear All,
There is a new version of the ZeligChoice package now available for testing. ZeligChoice is a Zelig module that contains limited-dependent variable choice models. This release fixes issues in the ordered probit model, and improves plot graphics for ordered logit, ordered probit and multinomial probit. If anyone has comments on these models, this package, any suggestions they would like to see, or anything related, do let me know.
You can install this version currently from:
install.packages("ZeligChoice", repos="http://r.iq.harvard.edu", type="source")
Depending on feedback, this will be posted to CRAN in about a week.
Thanks,
James
--
James Honaker, Senior Research Scientist
//// Institute for Quantitative Social Science, Harvard University
