search for: glm1

Hi - How can I 'manually' reproduce the results in 'pred1' below? My attempt is pred_manual, but is not correct. Any help is much appreciated. library(splines) set.seed(12345) y <- rgamma(1000, shape =0.5) age <- rnorm(1000, 45, 10) glm1 <- glm(y ~ ns(age, 4), family=Gamma(link=log)) dd <- data.frame(age = 16:80) mm <- model.matrix( ~ ns(dd$age, 4)) pred1 <- predict(glm1, dd, type = "response") pred_manual <- exp(coefficients(glm1)[1] * mm[,1] + coefficients(glm1)[2] * mm[,2] +...

...-value X 3 295 103.4730 <.0001 > anova(glmm2$lme, glmm1$lme) Model df AIC BIC logLik Test L.Ratio p-value glmm2$lme 1 3 4340.060 4351.172 -2167.030 glmm1$lme 2 5 4292.517 4311.036 -2141.258 1 vs 2 51.54367 <.0001 > > # a linear model also works, though no p-value is reported > glm1 <- gam(y ~ x1 + x2) > glm2 <- gam(y ~ 1) > anova(glm1) Family: gaussian Link function: identity Formula: y ~ x1 + x2 Parametric Terms: df F p-value x1 1 45.58 7.69e-11 x2 1 13.96 0.000224 > anova(glm2, glm1) Analysis of Deviance Table Model 1: y ~ 1 Model 2: y ~ x1 + x2 Resid...

...iones. Desde ya muchas gracias y disculpas si la pregunta es muy básica, adjunto los comandos que estoy utilizando. Si hiciera falta el archivo de texto puedo adjuntarlo. Desde ya mcuhas gracias library(MuMIn) archivo<-read.table("C:\\Users\\R\\arch.txt",header=T) attach(C.narino) glm1=glm(Frec~Sexo*Hora*Comportamiento, family=poisson) summary(glm1) anova(glm1,test="Chisq") dredge(glm1, rank = "AIC") [[alternative HTML version deleted]]

...iones. Desde ya muchas gracias y disculpas si la pregunta es muy básica, adjunto los comandos que estoy utilizando. Si hiciera falta el archivo de texto puedo adjuntarlo. Desde ya mcuhas gracias library(MuMIn) archivo<-read.table("C:\\Users\\R\\arch.txt",header=T) attach(C.narino) glm1=glm(Frec~Sexo*Hora*Comportamiento, family=poisson) summary(glm1) anova(glm1,test="Chisq") dredge(glm1, rank = "AIC") . [[alternative HTML version deleted]]

...iones. Desde ya muchas gracias y disculpas si la pregunta es muy básica, adjunto los comandos que estoy utilizando. Si hiciera falta el archivo de texto puedo adjuntarlo. Desde ya mcuhas gracias library(MuMIn) archivo<-read.table("C:\\Users\\R\\arch.txt",header=T) attach(C.narino) glm1=glm(Frec~Sexo*Hora*Comportamiento, family=poisson) summary(glm1) anova(glm1,test="Chisq") dredge(glm1, rank = "AIC") . [[alternative HTML version deleted]]

...8 ### sandwichGLM.R system("wget http://www.ats.ucla.edu/stat/stata/faq/eyestudy.dta") library(foreign) dat <- read.dta("eyestudy.dta") ### Ach, stata codes factor contrasts backwards dat$carrot0 <- ifelse(dat$carrot==0,1,0) dat$gender1 <- ifelse(dat$gender==1,1,0) glm1 <- glm(lenses~carrot0, data=dat, family=poisson(link=log)) summary(glm1) library(sandwich) vcovHC(glm1) sqrt(diag(vcovHC(glm1))) sqrt(diag(vcovHC(glm1, type="HC0"))) ### Result: # > summary(glm1) # Call: # glm(formula = lenses ~ carrot0, family = poisson(link = log), # d...

...iduals should all be restricted in the range 0 to 1 (since I am predicting a binary outcome). I read many posts on this list and I realized that there are four(!?) different types of residuals. I need a simple account of these four types of residuals, if anyone can help it will be great. residuals(glm1, "response") residuals(glm1, "pearson") residuals(glm1, "deviance") residuals(glm1, "working") - especially this one confuses me a lot! What is the "working" option and how is this different? Thank you Jason Dr. Iasonas Lamprianou Assistant Pr...

...teractions. Why does the second one specify a 1 degree of freedom test? The code below should help clarify my question. Thanks much, Danny ## I'm currently using: R Version 1.9.1 / Windows 2000 / P4/2.8 Ghz Z <- as.factor(rep(LETTERS[1:3],20)) Y <- rep(0:1, 30) X <- rnorm(60) glm1 <- glm(Y ~ Z + X + Z:X, family=binomial) summary(glm1)$coeff regTermTest( model=glm1 , test.terms=~Z:X) ZB <- ifelse(Z=="B",1,0) ZC <- ifelse(Z=="C",1,0) glm2 <- glm(Y ~ ZB + ZC + X + ZB:X + ZC:X, family=binomial) summary(glm2)$coeff ## Okay, same as glm1 regTe...

I am getting a repeated error when I try to run a logistic regression in R 2.8.1 >(glm(prop1~x1,data=glm1,family=binomial("logit"),weights=nt1)) Error in model.frame.default(formula = prop1 ~ x1, data = glm1, weights = nt1, : invalid type (list) for variable 'x1' x1 is multistate categorical (3 categories). 2 of the categories have 12 observation, one has 9. Is this what it is ob...

...13-233-5101 University of Leeds fax (44) 113-233-5090 Leeds LS2 9JT, England e-mail: J.T.Kent@leeds.ac.uk Input: x_1:4 # regressor variable y_c(2,6,7,8) # response binomial counts n_rep(10,4) # number of binomial trials ym_cbind(y,n-y) # response variable as a matrix glm1_glm(ym~x,binomial) # fit a generalized linear model f_fitted(glm1) rp1_(y-n*f)/sqrt(n*f*(1-f)) # direct calculation of pearson residuals rp2_residuals(glm1,type="pearson") # should be pearson residuals rd_residuals(glm1) # deviance residuals cbind(rp1,rp2,rd) # note second column is wrong...

Hello, I estimaded two logit models via glm(). A null model (called glm00) and "full" model with all accessible covariates and interactions between them (glm1). Then I tried to get even better model by step procedure. I tried the following code: > step(glm00, scope=formula(glm1), method="both") and another one: > step(glm00, scope=formula(glm1), method="forward") In both cases step procedure terminated after several...

Dear R-Users, I have a problem on extracting T-Stat and P-Value. I have written R-code below library("Matching") data("lalonde") attach(lalonde) names(lalonde) Y <- lalonde$re78 Tr <- lalonde$treat glm1 <- glm(Tr~age+educ+black+hisp+married+nodegr+re74+re75,family=binomial,data=lalonde) pscore.predicted <- predict(glm1) rr1 <- Match(Y=Y,Tr=Tr,X=glm1$fitted,estimand="ATT", M=1,ties=TRUE,replace=TRUE) summary(rr1) > summary(rr1) Estimate...  2624.3  AI SE......  802.19  T-...

Hi, all I have a question on mgcv and ns. Now I want to compare the results from glm, gam and ns. Take a simple model y~x for example. glm1 = glm(y~x, data=data1) gam1 = gam(y~s(x), data=data1) ns1 = glm(y~ns(x),data=data1) In order to confirm the result from glm1 is consistent to those from gam1 and ns1, I want to define degree=1 in mgcv and ns. I am wondering if there is somebody can give me some suggestions? I appreciate. Hellen...

Appologies if this has been addressed before, but I can't seem to find it in the help archives. I'm looking to do something like the following but it looks like save.plot is deprecated. save.plot(plot(glm1$residuals,gain,main = "Hist of residuals and gain"),file="Desktop/hist1.png") Thanks in advance, Sean Session Info: R version 2.9.1 (2009-06-26) i386-apple-darwin8.11.1 locale: en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8 attached base packages: [1] stats graph...

Dear Listers, If I have several glm objects with names glm1, glm2.... and want to apply new data to these objects. Instead of typing "predict(glm1, newdata)..." 100 times, is there way I could do so in a loop? Thank you so much! wensui [[alternative HTML version deleted]]

...yet explained by another predictor. As far as I understand Baron and Li describe how to do this, but my question is now: how do I report this in an article? Can anyone recommend a particular article that shows a concrete example of how the results from te following simple modeling can be reported: glm1 = glm(DV ~ A, family = binomial) glm2 = glm(DV ~ A + B, family = binomial) anova(glm1, glm2, test = "Chisq") Any help on how this simple kind of modeling should be reported is appreciated. Greetings, Edwin Commandeur [[alternative HTML version deleted]]

...e this is an environment problem, but I do not have a solution. Any assistance would be greatly appreciated. The function is as follows: > colonization <- function(x,y){ library(segmented) col <- subset(final,final[,x] == 2 | final[,x] == 4) col[,x] <- ifelse(col[,x] == 2,1,0) glm1 <- glm(col[,x]~col[,"PLAND"],family=binomial) glm2 <- segmented.glm(obj = glm1, seg.Z = ~PLAND, psi=y, it.max = 100) } > colonization("PIWO",48) Error in eval(expr, envir, enclos) : object "x" not found Thank you.

dear list! i have run several glm analysises to estimate a mean rate of dung decay for independent trials. i would like to compare these results statistically but can't find any solution. the glm calls are: dung.glm1<-glm(STATE~DAYS, data=o_cov, family="binomial(link="logit")) dung.glm2<-glm(STATE~DAYS, data=o_cov_T12, family="binomial(link="logit")) as all the trials have different sample sizes (around 20 each), anova(dung.glm1, dung.glm2) is not applicable. has anyone a...

...ause I don't understand what the fonction need in the argument "coeff". The two factors are Y (Ya, Yb,..., Ye) and Auto (Auto1, Auto2,...., Auto10) The responses are binomial (counts in the matrice "mat") and gaussian (variable "tot"). I use the glm procedure : glm1 <- glm (mat ~ y*auto, data, family=binomial) glm2 <- glm (tot ~ y*auto, data, family=gaussian) Subsequently I use the Type III anova : Anova (glm1, type="III") Anova (glm2, type="III") Using the contrast.treatment for both factors Y and Auto, R gives coherent results,...

...9, 8, 4, 14, 7, 3, 3, 9, 12, 8, 4, 6) failures <- c(14, 4, 9, 6, 9, 4, 7, 7, 6, 4, 8, 8, 9, 5, 8, 11, 13, 13, 3, 5, 7, 11, 11, 8, 6, 12, 11, 4, 13, 11, 6, 9, 12, 12, 7, 14, 4, 9, 11, 13, 11, 12, 16, 6, 13, 17, 17, 11, 8, 12, 16, 14) y <- cbind(successes, failures) data <- data.frame(y, x) glm1 <- glm(y ~ x, family= quasibinomial, data= data) glm2 <- glm(y ~ log(x), family=quasibinomial, data= data) # residual deviance is lower with log transformed x-value plot(x, successes) lines(x, predict(glm1, type= "response"), lwd=2) Firstly, because of the skewed distributio...