gmane.comp.lang.r.zelig

Re: Value for log(oratio) in "blogit" model

David Auty — 2012-05-09T22:27:46

Many thanks Olivia for the quick response, and apologies for not making myself clearer. 

I have worked through the whole 3 steps and simulated fitted values across the range of the explanatory variable in my data. That worked fine. The problem is we have been asked to supply equations by reviewers, so when we looked at the documentation where the equations are provided, it was unclear which value for psi we should use because of the x3 variable in the equation. The inverse link part was no problem, rather it was the manual calculation of joint probabilities that caused me to come unstuck. If we knew which part of the model output gives x3 (assuming log(oratio) is given by exp((intercept):3)), then I think that would solve our problem.

In no way would I want to (or be able to) reengineer the prediction process!

Thanks for any clarifications, and please excuse my lack of base knowledge about this kind of model.

Dave


-----Original Message-----
From: monkeykupo-Re5JQEeQqe8AvxtiuMwx3w< at >public.gmane.org [mailto:monkeykupo-Re5JQEeQqe8AvxtiuMwx3w< at >public.gmane.org] On Behalf Of Olivia Lau
Sent: May-09-12 4:50 PM
To: David Auty
Cc: zelig-boHuZhbfW+Y6DHxCPSz5rL5+avTOL6Rl< at >public.gmane.org
Subject: Re: [zelig] Value for log(oratio) in "blogit" model

I would suggest that you run through the whole 3 steps:  zelig, setx, and sim -- sim will output the predicted probabilities so you can examine further, no need to re-engineer the predicted probability process.  :)

You can also use the predict() function in the VGAM library, which just takes the inverse link function.  To find the inverse link function, take the output from vglm (say vout) and use:

vout< at >family< at >inverse

The inverse function takes the linear predictor (x %*% beta), but only works (well) in predict because additional objects are not called with proper lexical scoping (and are passed implicitly from the vout object itself).  But if you want to learn how predicted probabilities work, I would start here.

HTH.

Olivia

On Wed, May 9, 2012 at 1:27 PM, David Auty <auty.david.1-qLLbaQniQ8D3fQ9qLvQP4Q< at >public.gmane.org> wrote:

Re: Value for log(oratio) in "blogit" model

Olivia Lau — 2012-05-09T20:50:18

I would suggest that you run through the whole 3 steps:  zelig, setx,
and sim -- sim will output the predicted probabilities so you can
examine further, no need to re-engineer the predicted probability
process.  :)

You can also use the predict() function in the VGAM library, which
just takes the inverse link function.  To find the inverse link
function, take the output from vglm (say vout) and use:

vout< at >family< at >inverse

The inverse function takes the linear predictor (x %*% beta), but only
works (well) in predict because additional objects are not called with
proper lexical scoping (and are passed implicitly from the vout object
itself).  But if you want to learn how predicted probabilities work, I
would start here.

HTH.

Olivia

On Wed, May 9, 2012 at 1:27 PM, David Auty <auty.david.1-qLLbaQniQ8D3fQ9qLvQP4Q< at >public.gmane.org> wrote:

Value for log(oratio) in "blogit" model

David Auty — 2012-05-09T20:27:24

Hi there,

I'm fitting a bivariate logistic (blogit) model using Zelig and would like to use the equations given in the documentation to hard code some predictions, and to understand further how the model works.  The problem is that while it is easy to calculate the marginal probabilities, it is unclear from the documentation what value of psi to use in the equations for the joint probabilities. The documentation gives:

Odds ratio (psi )= exp(x3*beta3)

But we don't know what to substitute for x3 in this equation. Using exp((intercept):3) from the model does not give plausible answers, which leads us to assume we are missing something obvious and important.

Have searched the archives and also the VGAM documentation with no luck, so would be grateful for any pointers.

R2.15.0 on Windows 7.

Thanks in advance for any help.

Regards,

David Auty
Stagiaire Postdoctoral
Département des sciences du bois et de la forêt
Université Laval
Québec (QC) G1V 0A6

Gamma regression

Salvador Becerra-Munoz — 2012-04-30T20:53:45

Please I need your help!

This code sometimes predicts negative values because it is based on the inverse 
link. What is the modification to calculate de parameters using the loglink 
option?



 z.out = zelig(PT ~ Q + S + T + TDG, model = "gamma", data = PT)
 summary(z.out)
 x.out = setx(z.out)
 s.out = sim(z.out, x = x.out)
 plot(s.out)



 Thanks a lot for your advice

Conditional Frailty Models in Zelig

Brenton Peterson — 2012-04-27T01:20:25

Hi,

I am trying to estimate a conditional frailty survival model, as in
Box-Steffensmeier,  et al (2007 Political Analysis).  I have run them
previously in R but would like to implement them in Zelig for
post-estimation purposes.  Is there a way to do this?  While the Cox PH
with time-varying covariates is closest to my needs, it does not completely
fit the bill.  Is this an option in the current Zelig version?

Thanks,
B Peterson

Re: readRDS

Matt Owen — 2012-04-23T20:52:00

Hi - try installing the version found on CRAN with:
  install.packages("Zelig")

On Apr 23, 2012, at 4:38 PM, MP Sylvestre wrote:

Re: readRDS

MP Sylvestre — 2012-04-23T20:38:01

patch available tomorrow.
getting an error message:


Hello,
I am using R 15.0 and I have just downloaded Zelig from the harvard mirror.

#install.packages("Zelig", repos="http://r.iq.harvard.edu/", type="source")

However, when I try to load Zelig, I get an error similar to that described in
the initial post.


## In French (translation below)
Error : .onAttach a échoué dans attachNamespace() pour 'Zelig', détails :
  appel : NULL
  erreur : '.readRDS' est obsolète et retiré de R.
Utilisez plutôt 'readRDS'
Voir help("Defunct")
Erreur : le chargement du package / espace de noms a échoué pour ‘Zelig’

## 
Error : .onAttach failed in attachNamespace() for 'Zelig', details :
  call : NULL
  errorr : '.readRDS' is obsolete et removed from R.

Erreur : loading of package failed for ‘Zelig’


Where can I find the patch that was mentionned by  your email on April 5th?

Best,

MP





-
--
Zelig Mailing List, served by HUIT
Send messages: zelig< at >lists.gking.harvard.edu
[un]subscribe Options: http://lists.gking.harvard.edu/mailman/listinfo/zelig
Zelig program information: http://gking.harvard.edu/zelig/
Zelig mailing list
Zelig< at >lists.gking.harvard.edu
https://lists.gking.harvard.edu/mailman/listinfo/zelig

Cox Proportional Hazards Regression

Johannes Castner — 2012-04-08T03:15:31

To whom this may concern,

I am trying to run a Cox Proportional Hazards Regression in the 
following way, using Zelig:

z.out <- zelig(Surv(duration, Censored) ~ Ethnic1 + Type + Year1, model 
= "coxph", data = data)

but I get the following error:

Error in fitter(X, Y, strats, offset, init, control, weights = weights,  :
   NA/NaN/Inf in foreign function call (arg 6)
In addition: Warning message:
In fitter(X, Y, strats, offset, init, control, weights = weights,  :
   Ran out of iterations and did not converge

Does anyone know just what that means and how I can get around it?

Johannes

AUTO: Tom Aldenberg is out of the office (returning Fri04/13/2012)

Tom Aldenberg — 2012-04-06T09:48:24

I am out of the office from Sat 03/31/2012 until Fri 04/13/2012.




Note: This is an automated response to your message  "Zelig Digest, Vol 20,
Issue 1" sent on 5-4-2012 18:00:02.

This is the only notification you will receive while this person is away.

Re: readRDS

Owen, Matthew — 2012-04-04T23:23:31

It seems that R is no longer supporting the readRDS function. There will be a patch available tomorrow. Thanks for the notification.
- M
________________________________________
From: zelig-bounces-boHuZhbfW+Y6DHxCPSz5rL5+avTOL6Rl< at >public.gmane.org [zelig-bounces-boHuZhbfW+Y6DHxCPSz5rA< at >public.gmane.orgard.edu] On Behalf Of Matthias Kretschmer [mkretschmer2012-Re5JQEeQqe8AvxtiuMwx3w< at >public.gmane.org]
Sent: Wednesday, April 04, 2012 3:13 PM
To: zelig-boHuZhbfW+Y6DHxCPSz5rL5+avTOL6Rl< at >public.gmane.org
Subject: [zelig] readRDS

I tried to install Zelig on R 2.15 (on a Windows 7, 64bit machine) and keep getting an error message:

Error : .onAttach failed in attachNamespace() for 'Zelig', details:
  call: NULL
  error: '.readRDS' is defunct.
Use 'readRDS' instead.
See help("Defunct")
Error: package/namespace load failed for ‘Zelig’

Any idea what that means?

Matthias

readRDS

Matthias Kretschmer — 2012-04-04T19:13:37

I tried to install Zelig on R 2.15 (on a Windows 7, 64bit machine) and keep
getting an error message:

*Error : .onAttach failed in attachNamespace() for 'Zelig', details:
  call: NULL
  error: '.readRDS' is defunct.
Use 'readRDS' instead.
See help("Defunct")
Error: package/namespace load failed for ‘Zelig’
*
Any idea what that means?

Matthias

Re: summary with imputations bug ... any workarounds?

Valeriy — 2012-03-08T22:45:03


Is there any update or workaround? Amelia imputed data but I don't seem 
to be able to use it as with zelig logit.mixed model. 

Valeriy

Re: Plot.ci function to logit.survey

Matt Owen — 2012-03-02T22:24:21

Thanks! I'll get this into including this into the next release.

On Mar 2, 2012, at 3:16 PM, Rogério Barbosa wrote:

Plot.ci function to logit.survey

Rogério Barbosa — 2012-03-02T20:16:26

Hey all again,

I would like to inform that the most updated version of plot.ci function
doesn't have the item "logit.survey" listed on its algorithm. So it is not
possible to plot confidence intervals for predicted values calculated using
this model.

I tried to make that inclusion myself -- and it worked.

Here it goes the modified function.
(I found the original at:
http://r.iq.harvard.edu/src/contrib/Zelig/R/plot.ci.R)

Best,





plot.ci <- function(x, CI=95, qi = "ev", main = "",

                    ylab = NULL, xlab = NULL, xlim = NULL,

                    ylim = NULL, col = c("red", "blue"), ...) {

  "%w/o%" <- function(x,y) x[!x %in% y] #--  x without y

  if (class(x) != "zelig")

    stop(" plot.ci() works only for sim() output.")

  if (!(x$zelig.call$model) %in%

        c("ls", "logit",* "logit.survey"*, "normal.survey","probit", "exp",
"gamma", "lognorm",

          "weibull", "normal", "poisson", "tobit", "relogit",

          "negbin", "logit.bayes", "probit.bayes",

          "poisson.bayes", "normal.bayes", "tobit.bayes",

          "ls.mixed", "logit.mixed", "probit.mixed",

          "gamma.mixed", "poisson.mixed",

          "logit.gam", "gamma.gee", "normal.gam", "poisson.gam",

          "probit.gam", "logit.gee", "normal.gee",

          "poisson.gee", "probit.gee", "normal.survey"))

    stop("\n  plot.ci() is valid only for non-categorical, univariate
response models.")

  cip <- c((100-CI)/200, 1-(100-CI)/200)

  summarize <- function(z, cip){

    res <- NULL

    res <- cbind(res, apply(z, 2, quantile, prob=cip[1]))

    res <- cbind(res, apply(z, 2, quantile, prob=cip[2]))

    res

  }

  vv <- apply(x$x, 2, unique)

  idx <- sapply(vv, length)

  cidx <- which(idx > 1)

  if (!is.null(x$x1)) {

    vv1 <- apply(x$x1, 2, unique)

    idx1 <- sapply(vv1, length)

    cidx1 <- which(idx1 > 1)

    if (!identical(names(idx), names(idx1)))

      stop("variables in x and x1 do not match.")

    ## Checking for one dimension of variation, including interaction terms

    if (length(cidx) > length(cidx1)) {

      tmp <- names(idx)[cidx %w/o% cidx1]

      tmp1 <- names(idx)[cidx[cidx %in% cidx1]]

    }

    else {

      tmp <- names(idx1)[cidx1 %w/o% cidx]

      tmp1 <- names(idx1)[cidx1[cidx1 %in% cidx]]

    }

  check <- grep(tmp1, tmp)

  if (length(check) != length(tmp))

    stop("x and x1 vary on more than one dimension.")

  }

  var <- vv[[cidx[1]]]

  q <- pmatch(qi, names(x$qi))

  qofi <- x$qi[[q]]

  sum.qi <- summarize(qofi, cip)

  if (!is.null(x$x1) && qi == "ev") {

    fd <- x$qi$fd

    ev1 <- fd + qofi

    sum.qi1 <- summarize(ev1, cip)

  }

  else sum.qi1 <- NULL

  if (is.null(ylab))  ylab <- x$qi.name[[q]]

  if (is.null(xlab))  xlab <- paste("Range of", colnames(x$x)[cidx[1]])

  if (is.null(ylim)) {

    if (is.null(sum.qi1))  ylim <- c(min(sum.qi), max(sum.qi))

    else  ylim <- c(min(sum.qi, sum.qi1), max(sum.qi, sum.qi1))

  }

  if (is.null(xlim))  xlim <- c(min(var), max(var))

  plot.default(var, type = "n", ylab = ylab, main = main, xlab = xlab,

               xlim = xlim, ylim = ylim)

  for (i in 1:length(var)) {

    lines(c(var[i], var[i]), c(sum.qi[i,1], sum.qi[i,2]), col = col[1], ...)

    if (!is.null(x$x1) && qi == "ev")

      lines(c(var[i], var[i]), c(sum.qi1[i,1], sum.qi1[i,2]), col = col[2],
...)

  }

}





- - - - -
Rogério J. Barbosa
Researcher at Centre for Metropolitan Studies/Cebrap
São Paulo - Brazil

Re: Fwd: setx "incorrect number of dimensions" error

Matt Owen — 2012-03-02T16:33:01

Good to hear!

On Mar 2, 2012, at 11:28 AM, Rogério Barbosa wrote:

Re: Fwd: setx "incorrect number of dimensions" error

Rogério Barbosa — 2012-03-02T16:28:27

Now it worked! It was an stata database, readed using "foreign". I
converted the data to csv. Problem solved.

Thank you, Matt.

2012/3/2 Matt Owen <mowen-zvdtVsPQyYJu6/7BKvT5pw< at >public.gmane.org>

Re: Fwd: setx "incorrect number of dimensions" error

Matt Owen — 2012-03-02T15:59:56

The only other though that I have is that "pda" is not a standard data.frame object. If this is the case, you may have to restructure the data set before using it with the logit model.

On Mar 2, 2012, at 10:54 AM, Rogério Barbosa wrote:

Fwd: setx "incorrect number of dimensions" error

Rogério Barbosa — 2012-03-02T15:54:50

Hi Matt,

Thanks for your help! But, unfortunately, it didn't work. The error keep
happening even if I use just one independent variable.
For instance:

*Error in dta[complete.cases(mf), names(dta) %in% vars, drop = FALSE] : *
*  incorrect number of dimensions*

Just an observation: it stated to occur suddenly. The day before the model
was working fine.

- - - - -
Rogério J. Barbosa
Researcher at Centre for Metropolitan Studies/Cebrap
São Paulo - Brazil






2012/3/2 Matt Owen <mowen-zvdtVsPQyYJu6/7BKvT5pw< at >public.gmane.org>

Re: setx "incorrect number of dimensions" error

Matt Owen — 2012-03-02T14:35:05

My best guess is that the error occurs because the values for "gender" and "race.white" have a different length than "age" and "age2"

Try something like:

setx(z.out, age=age.sim, age2=age2.sim, gender=rep(0, length(age.sim)), race.white=rep(1, length(age.sim)))


On Mar 1, 2012, at 11:06 PM, Rogério Barbosa wrote:

setx "incorrect number of dimensions" error

Rogério Barbosa — 2012-03-02T04:06:02

Hey all,

I get this error when I run setx: "Error in dta[complete.cases(mf),
names(dta) %in% vars, drop = FALSE] : incorrect number of dimensions"
I found some people with the same problem on internet - but apparently no
one has worked it out.
Does anyone know how to solve it?

My output is bellow the message.

Thanks,
Rogério Barbosa


data = pad)
How to cite this model in Zelig:
Kosuke Imai, Gary King, and Oliva Lau. 2008. "logit: Logistic Regression
for Dichotomous Dependent Variables" in Kosuke Imai, Gary King, and Olivia
Lau, "Zelig: Everyone's Statistical Software,"
http://gking.harvard.edu/zelig


*Error in dta[complete.cases(mf), names(dta) %in% vars, drop = FALSE] : *
*  incorrect number of dimensions*



Rogério J. Barbosa
Researcher at Centre for Metropolitan Studies/Cebrap
São Paulo - Brazil

Error with ls.mixed and weights

Harden, Jeff — 2012-02-27T18:44:53

Hello,

I am having problems with the "ls.mixed" model in zelig() and the weights argument. This is similar to the problem documented here: https://lists.gking.harvard.edu/pipermail/zelig/2011-July/000183.html. However, while it seems like the issue in that thread was missing values on the dependent variable, I am having the issue with no missing data. I can estimate the model with weights using lmer(), but not using zelig(). I have an example that reproduces this error using fake data posted here: http://pastebin.com/29KbmRZD and listed below.

##
library(Zelig)
set.seed(1000)

x <- runif(1000, -1, 1)
y <- x + rnorm(1000)
wt <- rexp(1000)
id <- rep(1:20, times = 50)
d <- as.data.frame(cbind(y, x, wt, id))

m1 <- zelig(y ~ x + tag(1|id), model = "ls.mixed", weights = "wt", data = d)
m2 <- lmer(y ~ x + (1|id), weights = wt, data = d)
##

Thanks for any insight with this problem.

Jeff Harden


---
Jeffrey J. Harden
Ph.D Candidate, Department of Political Science
Program Coordinator, American Politics Research Group (APRG)
University of North Carolina at Chapel Hill
jjharden-gpQNmkgQrW4< at >public.gmane.org<mailto:jjharden-gpQNmkgQrW4< at >public.gmane.org>
http://jjharden.web.unc.edu/