← Back to branch summary

~ubuntu-branches/ubuntu/wily/r-cran-gss/wily-proposed

« back to all changes in this revision

Viewing changes to R/family.R

  • Committer: Package Import Robot
  • Author(s): Dirk Eddelbuettel
  • Date: 2014-12-06 13:11:29 UTC
  • mfrom: (1.1.8)
  • Revision ID: package-import@ubuntu.com-20141206131129-3crdxrzw4j6s62wq
Tags: 2.1-4-1
* New upstream release

* debian/control: Set Build-Depends: to current R version 
* debian/control: Set Standards-Version: to current version 

Show diffs side-by-side

added added

removed removed

Lines of Context:
R/family.R
212
212
            stop("gss error: negative binomial response should be nonnegative")
213
213
        if (min(y[,2])<=0)
214
214
            stop("gss error: negative binomial size should be positive")
215
 
        p <- 1-1/(1+exp(eta))
216
 
        u <- (y[,1]+y[,2])*p-y[,2]
217
 
        w <- y[,2]*(1-p)
 
215
        odds <- exp(eta)
 
216
        p <- odds/(1+odds)
 
217
        q <- 1/(1+odds)
 
218
        u <- y[,1]*p-y[,2]*q
 
219
        w <- y[,2]*q
218
220
        ywk <- eta-u/w-offset
219
221
        wt <- w*wt
220
 
        list(ywk=ywk,wt=wt)
 
222
        list(ywk=ywk,wt=wt,u=u*wt)
221
223
    }
222
224
    else {
223
225
        if (min(y)<0)
224
226
            stop("gss error: negative binomial response should be nonnegative")
225
 
        p <- 1-1/(1+exp(eta))
226
 
        if (is.null(nu)) log.nu <- log(mean(y*exp(eta)))
 
227
        odds <- exp(eta)
 
228
        p <- odds/(1+odds)
 
229
        q <- 1/(1+odds)
 
230
        if (is.null(nu)) log.nu <- log(mean(y*odds))
227
231
        else log.nu <- log(nu)
228
232
        repeat {
229
233
            nu <- exp(log.nu)
233
237
            if (abs(log.nu-log.nu.new)/(1+abs(log.nu))<1e-7) break
234
238
            log.nu <- log.nu.new
235
239
        }
236
 
        u <- (y+nu)*p-nu
237
 
        w <- nu*(1-p)
 
240
        u <- y*p-nu*q
 
241
        w <- nu*q
238
242
        ywk <- eta-u/w-offset
239
243
        wt <- w*wt
240
244
        list(ywk=ywk,wt=wt,nu=nu,u=u*wt)
245
249
dev.resid.nbinomial <- function(y,eta,wt)
246
250
{
247
251
    if (is.null(wt)) wt <- rep(1,dim(y)[1])
248
 
    p <- 1-1/(1+exp(eta))
249
 
    as.vector(2*wt*(y[,1]*log(ifelse(y[,1]==0,1,y[,1]/(y[,1]+y[,2])/(1-p)))
 
252
    odds <- exp(eta)
 
253
    p <- odds/(1+odds)
 
254
    q <- 1/(1+odds)
 
255
    as.vector(2*wt*(y[,1]*log(ifelse(y[,1]==0,1,y[,1]/(y[,1]+y[,2])/q))
250
256
                    +y[,2]*log(y[,2]/(y[,1]+y[,2])/p)))
251
257
}
252
258
 
258
264
    if (!is.null(offset)) {
259
265
        eta <- log(p/(1-p)) - mean(offset)
260
266
        repeat {
261
 
            p <- 1-1/(1+exp(eta+offset))
262
 
            u <- (y[,1]+y[,2])*p-y[,2]
263
 
            w <- y[,2]*(1-p)
 
267
            odds <- exp(eta+offset)
 
268
            p <- odds/(1+odds)
 
269
            q <- 1/(1+odds)
 
270
            u <- y[,1]*p-y[,2]*q
 
271
            w <- y[,2]*q
264
272
            eta.new <- eta-sum(wt*u)/sum(wt*w)
265
273
            if (abs(eta-eta.new)/(1+abs(eta))<1e-7) break
266
274
            eta <- eta.new    
267
275
        }
268
276
    }
269
 
    sum(2*wt*(y[,1]*log(ifelse(y[,1]==0,1,y[,1]/(y[,1]+y[,2])/(1-p)))
 
277
    sum(2*wt*(y[,1]*log(ifelse(y[,1]==0,1,y[,1]/(y[,1]+y[,2])/q))
270
278
              +y[,2]*log(y[,2]/(y[,1]+y[,2])/p)))
271
279
}