

Dear RUsers,
without going into details I tried to prepare a simple example to show
you where I would need help.
In particular I prepare two examplestemplate for a study I'm conduction
on discretetime methods for survival analysis.
Each of this example has two datasets which are basically equal, with
the exception that in the former one has individual data and in the
latter one aggregated data.
The difference between the two examples is on a single subject: I
substituted to the first example a censored case with a subject who died
at the first timeunit.
Afterward I fitted a logistic model (Fahrmeir and Tutz, 2001) in the glm
context, but whereas there is not difference between individual and
aggregated dataset in the first example, I noted some discrepancies in
the second example. I might guess that something with weights is going
on, but I did not manage to clearly understand.
Hope that the following example will be more clear than my explanations,
Thanks in advance,
Carlo Giovanni Camarda
rm(list = ls())
# working one
timesIND < c(rep(1:4, 3), 1, rep(1:2,2), rep(1:3 , 2), rep(1:4,
2))
statusIND < c(rep(0 ,12), 1, rep(0:1,2), rep(c(0,0,1), 2),
rep(c(0,0,0,1),2))
datiIND < as.data.frame(cbind(timesIND, statusIND))
datiIND$timesIND < as.factor(datiIND$timesIND)
timesAGG < c( 1:4, 1, 1:2, 1:3, 1:4)
statusAGG < c(rep(0,4), 1, 0:1, c(0,0,1), c(0,0,0,1))
weightAGG < c(rep(3,4), 1, rep(2,2), rep(2,3), rep(2,4))
datiAGG < as.data.frame(cbind(timesAGG, statusAGG, weightAGG))
datiAGG$timesAGG < as.factor(datiAGG$timesAGG)
coef(glm(statusIND ~ timesIND, family=binomial, data=datiIND))
coef(glm(statusAGG ~ timesAGG, family=binomial, data=datiAGG,
weights=weightAGG))
# not working one
timesINDa < c(rep(1:4, 4), rep(1:2,2), rep(1:3 , 2), rep(1:4,
2))
statusINDa < c(rep(0 ,16), rep(0:1,2), rep(c(0,0,1), 2),
rep(c(0,0,0,1),2))
datiINDa < as.data.frame(cbind(timesINDa, statusINDa))
datiINDa$timesINDa < as.factor(datiINDa$timesINDa)
timesAGGa < c( 1:4, 1:2, 1:3, 1:4)
statusAGGa < c(rep(0,4), 0:1, c(0,0,1), c(0,0,0,1))
weightAGGa < c(rep(4,4), rep(2,2), rep(2,3), rep(2,4))
datiAGGa < as.data.frame(cbind(timesAGGa, statusAGGa, weightAGGa))
datiAGGa$timesAGGa < as.factor(datiAGGa$timesAGGa)
coef(glm(statusINDa ~ timesINDa, family=binomial, data=datiINDa))
coef(glm(statusAGGa ~ timesAGGa, family=binomial, data=datiAGGa,
weights=weightAGGa))
+++++
This mail has been sent through the MPI for Demographic Rese...{{dropped}}
______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide! http://www.Rproject.org/postingguide.html


You've found a region of infinite extent over which the likelihood
function is for all practical purposes flat. This means that the
maximum likelihood estimates (MLEs) are not unique. To see this
consider the following properties of your datiINDa:
> with(datiINDa, table(statusINDa, timesINDa))
timesINDa
statusINDa 1 2 3 4
0 10 8 6 4
1 0 2 2 2
> sapply(datiINDa, class)
timesINDa statusINDa
"factor" "numeric"
You are estimating 4 parameters, an intercept plus one parameter for
each level of the factor "timesInda". The first level occurs only with
statusINDa = 0, never with statusINDa = 1. Therefore, the theoretical
MLE for that level of timesINDa would have slope = +/Inf (and the
intercept would also be adjusted to +/Inf to compensate). However, glm
doesn't bother pushing it that far, and gives up with still moderately
small values for the parameters. To understand this better, first
modify your example to store the glm fitted object as follows:
fit.a < glm(statusINDa ~ timesINDa, family=binomial, data=datiINDa)
Then apply "predict" to that object:
predict(fit.a, type="response")
The result is that the 10 cases with timesInda = 1 all have a
Pr{statusINDa = 1} = 3e9, which glm thinks is essentially 0 and quits.
Now let's do the same with your weighted version:
fit.wa < glm(statusAGGa ~ timesAGGa, family=binomial, data=datiAGGa,
weights=weightAGGa)
sort(predict(fit.wa, type="response"))
Those 10 cases now have Pr{statusINDa = 1} = 5.4e9. This is
essentially the issue of "complete separation". We can request more
precision as follows:
> fit.a3 < glm(statusINDa ~ timesINDa, family=binomial, data=datiINDa,
+ control=glm.control(epsilon=1e13,
+ maxit=250))
Warning message:
fitted probabilities numerically 0 or 1 occurred in: glm.fit(x = X, y =
Y, weights = weights, start = start, etastart = etastart,
In this case, we get a warning. For more on this, try
RSiteSearch("complete separation with logistic regression").
Sehr interessant, nicht?
hope this helps.
spencer graves
Camarda, Carlo Giovanni wrote:
> Dear RUsers,
> without going into details I tried to prepare a simple example to show
> you where I would need help.
> In particular I prepare two examplestemplate for a study I'm conduction
> on discretetime methods for survival analysis.
> Each of this example has two datasets which are basically equal, with
> the exception that in the former one has individual data and in the
> latter one aggregated data.
> The difference between the two examples is on a single subject: I
> substituted to the first example a censored case with a subject who died
> at the first timeunit.
> Afterward I fitted a logistic model (Fahrmeir and Tutz, 2001) in the glm
> context, but whereas there is not difference between individual and
> aggregated dataset in the first example, I noted some discrepancies in
> the second example. I might guess that something with weights is going
> on, but I did not manage to clearly understand.
> Hope that the following example will be more clear than my explanations,
> Thanks in advance,
> Carlo Giovanni Camarda
>
>
> rm(list = ls())
> # working one
>
> timesIND < c(rep(1:4, 3), 1, rep(1:2,2), rep(1:3 , 2), rep(1:4,
> 2))
> statusIND < c(rep(0 ,12), 1, rep(0:1,2), rep(c(0,0,1), 2),
> rep(c(0,0,0,1),2))
> datiIND < as.data.frame(cbind(timesIND, statusIND))
> datiIND$timesIND < as.factor(datiIND$timesIND)
>
> timesAGG < c( 1:4, 1, 1:2, 1:3, 1:4)
> statusAGG < c(rep(0,4), 1, 0:1, c(0,0,1), c(0,0,0,1))
> weightAGG < c(rep(3,4), 1, rep(2,2), rep(2,3), rep(2,4))
> datiAGG < as.data.frame(cbind(timesAGG, statusAGG, weightAGG))
> datiAGG$timesAGG < as.factor(datiAGG$timesAGG)
>
> coef(glm(statusIND ~ timesIND, family=binomial, data=datiIND))
> coef(glm(statusAGG ~ timesAGG, family=binomial, data=datiAGG,
> weights=weightAGG))
>
> # not working one
>
> timesINDa < c(rep(1:4, 4), rep(1:2,2), rep(1:3 , 2), rep(1:4,
> 2))
> statusINDa < c(rep(0 ,16), rep(0:1,2), rep(c(0,0,1), 2),
> rep(c(0,0,0,1),2))
> datiINDa < as.data.frame(cbind(timesINDa, statusINDa))
> datiINDa$timesINDa < as.factor(datiINDa$timesINDa)
>
> timesAGGa < c( 1:4, 1:2, 1:3, 1:4)
> statusAGGa < c(rep(0,4), 0:1, c(0,0,1), c(0,0,0,1))
> weightAGGa < c(rep(4,4), rep(2,2), rep(2,3), rep(2,4))
> datiAGGa < as.data.frame(cbind(timesAGGa, statusAGGa, weightAGGa))
> datiAGGa$timesAGGa < as.factor(datiAGGa$timesAGGa)
>
> coef(glm(statusINDa ~ timesINDa, family=binomial, data=datiINDa))
> coef(glm(statusAGGa ~ timesAGGa, family=binomial, data=datiAGGa,
> weights=weightAGGa))
>
> +++++
> This mail has been sent through the MPI for Demographic Rese...{{dropped}}
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/rhelp> PLEASE do read the posting guide! http://www.Rproject.org/postingguide.html______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide! http://www.Rproject.org/postingguide.html

