On Jul 5, 2009, at 5:15 PM, sjkimble wrote:
R Help:
When I look at the object S, I see about half of them are
zeroes.
I've address the object S being zero so often by changing the value of
alpha, which I'm allowed to do according to the author. All values
of S are
non zero and should not give Inf.
The "expectation" of a binomial model is for the LHS of the
formulat to be either a 1/0 vector or something of the form
cbind(events, non_events). You have satisfied that expectation with p
but there are only 2 such cases. It seems unreasonable to my thinking
to expect that a logistic regression model can deliver sensible
output
from only 2 events. And that holds doubly (or perhaps infinitely?)
true when you are starting out with half of your covariates equal to
log(0) = -Inf.
Object p is presence (1) or absence (0) of a species on the 12
islands, and
some of them are rare so they are mostly zeros. Nevertheless I tried
with a
more common species whose data are:
x.crd y.crd A p
1 361763 1034071 94.169 1
2 370325 1027277 127.642 1
3 370416 1027166 127.961 1
4 370471 1027148 1804.846 1
5 369050 1031312 1790.493 1
6 370908 1026354 199.103 1
7 361562 1034311 2047.637 1
8 365437 1022188 1622678.961 1
9 347047 1025334 21.169 1
10 349186 1024441 408556.801 1
11 361762 1034052 3.414 0
12 370799 1026557 103.994 0
So now you have gone from a situation with only two events to one with
only two non-events. That is not a situation that is really any better
from a statistical point of view.
The code, plus error message:
haliclona_reniera_manglaris<-read.table(file.choose(),header=T)
attach(haliclona_reniera_manglaris)
alphascan<-function(alpha,d,A,p){
+ edis<-as.matrix(exp(-alpha*d))
+ edis<-sweep(edis,2,A,"*")
+ S<-rowSums(edis[,p>0])
Hey,.... you told me you did <something> to get rid of the zeros in
S. When I take this code out of the functional wrapper I get:
> S
1 2 3 4
5 6 7 8
6.965905e-07 5.139157e-07 8.325874e-60 1.318603e-22
1.329174e-22 0.000000e+00 1.390849e-94 1.755094e-97
9 10 11 12
1.411310e-134 0.000000e+00 0.000000e+00 0.000000e+00
So now instead of 3/4 of them being zero only 1/3 of them are still
zero. That does not seem to be satisfactory.
+ mod<-glm(p~offset(2*log(S))+log(A),family=binomial)
+ deviance(mod)
+ }
(sol<-optimize(alphascan,c(0.001,10),d=d,A=A,p=p))
So maybe the problem remain the same as before and supplying the
missing d will not be enough. You still (I think) need a dataset that
supplies enough information for a statistical approach.
David Winsemius, MD
Heritage Laboratories
West Hartford, CT
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
