R16

Datos de proporcion

numbers <- read.table("sexratio.txt",header=T)
attach(numbers)
head(numbers)

##   density females males
## 1       1       1     0
## 2       4       3     1
## 3      10       7     3
## 4      22      18     4
## 5      55      22    33
## 6     121      41    80

par(mfrow=c(1,2))
p <- males/(males+females)
plot(density,p,ylab="Proportion male",pch=16,col="blue")
plot(log(density),p,ylab="Proportion male",pch=16,col="blue")

# unir ambos conteos en un solo vector

y <- cbind(males,females)
model <- glm(y~density,binomial)
# en el modelo la proprocion de sexos en función de la densidad
summary(model)

## 
## Call:
## glm(formula = y ~ density, family = binomial)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.4619  -1.2760  -0.9911   0.5742   1.8795  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 0.0807368  0.1550376   0.521    0.603    
## density     0.0035101  0.0005116   6.862 6.81e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 71.159  on 7  degrees of freedom
## Residual deviance: 22.091  on 6  degrees of freedom
## AIC: 54.618
## 
## Number of Fisher Scoring iterations: 4

# hay sobredisperción, por lo cual se hace una transforación logaritmica
model2 <- glm(y~log(density),binomial)
summary(model2)

## 
## Call:
## glm(formula = y ~ log(density), family = binomial)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9697  -0.3411   0.1499   0.4019   1.0372  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.65927    0.48758  -5.454 4.92e-08 ***
## log(density)  0.69410    0.09056   7.665 1.80e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 71.1593  on 7  degrees of freedom
## Residual deviance:  5.6739  on 6  degrees of freedom
## AIC: 38.201
## 
## Number of Fisher Scoring iterations: 4

## No hay patrones en los residuales vs ajustados, el ggrafico de normalidad es considerablemente lineal, el punto 4 tiene una distancia grande pero omitiendolo el modelo sige siendo significativo
plot(model2)

# La proporción de machos aumenta confomre aumenta la densidad poblacional

xv <- seq(0,6,0.01)
yv <- predict(model2,list(density=exp(xv)),type="response")
plot(log(density),p,ylab="Proportion male",pch=16,col="blue")
lines(xv,yv,col="red")

library(ggplot2)
library(plotly)

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

numbers <- read.table("sexratio.txt",header=T)
attach(numbers)

## The following objects are masked from numbers (pos = 5):
## 
##     density, females, males

head(numbers)

##   density females males
## 1       1       1     0
## 2       4       3     1
## 3      10       7     3
## 4      22      18     4
## 5      55      22    33
## 6     121      41    80

par(mfrow=c(1,2))
prob <- males/(males+females)
dataframalol <- data.frame(prob, numbers$density)


densa <- ggplot(data = dataframalol, aes(x=density, y= prob))+
  geom_point(colour = rgb(0.3, 0.7, 0.6))
fig1 <- ggplotly(densa)
densa2 <- ggplot(data = dataframalol, aes(x=log(density), y= prob))+
  geom_point(colour = rgb(0.2, 0.5, 0.7))

fig2 <- ggplotly(densa2)
fig3 <- subplot(fig1, fig2)
fig3

y <- cbind(males,females)
model <- glm(y~density,binomial)
model112 <- glm(y~log(density),binomial)

xvl <- seq(0,6,0.01)
yvl <- predict(model112,list(density=exp(xvl)),type="response")
dataele <- data.frame(xvl, yvl)


densa4 <- ggplot(data = dataframalol, aes(x=log(density), y= prob))+
  geom_point(colour = rgb(0.1, 0.4, 0.3))+
  geom_line(data = dataele, aes(x=xvl, y=yvl), colour = rgb(0.2, 0.7, 0.3))

fig4 <- ggplotly(densa4)
fig4

library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:plotly':
## 
##     select

data <- read.table("bioassay.txt",header=T)
attach(data)
names(data)

## [1] "dose"  "dead"  "batch"

y <- cbind(dead,batch-dead)
model <- glm(y~log(dose),binomial)

dose.p(model, p=c(0.5,0.9,0.95))

##               Dose         SE
## p = 0.50: 2.306981 0.07772065
## p = 0.90: 3.425506 0.12362080
## p = 0.95: 3.805885 0.15150043

germination <- read.table("germination.txt",header=T)
data2 <-as.data.frame(germination)
names(germination)

## [1] "count"     "sample"    "Orobanche" "extract"

germination

##    count sample Orobanche  extract
## 1     10     39       a75     bean
## 2     23     62       a75     bean
## 3     23     81       a75     bean
## 4     26     51       a75     bean
## 5     17     39       a75     bean
## 6      5      6       a75 cucumber
## 7     53     74       a75 cucumber
## 8     55     72       a75 cucumber
## 9     32     51       a75 cucumber
## 10    46     79       a75 cucumber
## 11    10     13       a75 cucumber
## 12     8     16       a73     bean
## 13    10     30       a73     bean
## 14     8     28       a73     bean
## 15    23     45       a73     bean
## 16     0      4       a73     bean
## 17     3     12       a73 cucumber
## 18    22     41       a73 cucumber
## 19    15     30       a73 cucumber
## 20    32     51       a73 cucumber
## 21     3      7       a73 cucumber

res <- germination$sample-germination$count
y <- cbind(germination$count, germination$sample-germination$count)

Oro <- factor(data2$Orobanche,
                 levels = c('a75','a73'))

levels(Oro)

## [1] "a75" "a73"

levels(germination$extract)

## NULL

model <- glm(y ~ germination$Orobanche * germination$extract, binomial)
summary(model)

## 
## Call:
## glm(formula = y ~ germination$Orobanche * germination$extract, 
##     family = binomial)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.01617  -1.24398   0.05995   0.84695   2.12123  
## 
## Coefficients:
##                                                      Estimate Std. Error
## (Intercept)                                           -0.4122     0.1842
## germination$Orobanchea75                              -0.1459     0.2232
## germination$extractcucumber                            0.5401     0.2498
## germination$Orobanchea75:germination$extractcucumber   0.7781     0.3064
##                                                      z value Pr(>|z|)  
## (Intercept)                                           -2.238   0.0252 *
## germination$Orobanchea75                              -0.654   0.5132  
## germination$extractcucumber                            2.162   0.0306 *
## germination$Orobanchea75:germination$extractcucumber   2.539   0.0111 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 98.719  on 20  degrees of freedom
## Residual deviance: 33.278  on 17  degrees of freedom
## AIC: 117.87
## 
## Number of Fisher Scoring iterations: 4

model1 <- glm(y ~ germination$Orobanche * germination$extract, quasibinomial)
model2 <- update(model1, ~ . - germination$Orobanche:germination$extract)
anova(model,model2,test="F")

## Warning: using F test with a 'binomial' family is inappropriate

## Analysis of Deviance Table
## 
## Model 1: y ~ germination$Orobanche * germination$extract
## Model 2: y ~ germination$Orobanche + germination$extract
##   Resid. Df Resid. Dev Df Deviance      F  Pr(>F)  
## 1        17     33.278                             
## 2        18     39.686 -1  -6.4081 6.4081 0.01136 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(model2,test="F")

## Analysis of Deviance Table
## 
## Model: quasibinomial, link: logit
## 
## Response: y
## 
## Terms added sequentially (first to last)
## 
## 
##                       Df Deviance Resid. Df Resid. Dev       F    Pr(>F)    
## NULL                                     20     98.719                      
## germination$Orobanche  1    2.544        19     96.175  1.1954    0.2887    
## germination$extract    1   56.489        18     39.686 26.5412 6.692e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model3 <- update(model2, ~ . - germination$Orobanche)
anova(model2,model3,test="F")

## Analysis of Deviance Table
## 
## Model 1: y ~ germination$Orobanche + germination$extract
## Model 2: y ~ germination$extract
##   Resid. Df Resid. Dev Df Deviance      F Pr(>F)
## 1        18     39.686                          
## 2        19     42.751 -1   -3.065 1.4401 0.2457

coef(model3)

##                 (Intercept) germination$extractcucumber 
##                  -0.5121761                   1.0574031

#se calcula p = 1/(1+(1/(e^x)))

calc<- 1/(1+1/(exp(-0.5122)))
calc

## [1] 0.3746779

calc2<- 1/(1+1/(exp(-0.5122+1.0574)))
calc2

## [1] 0.6330212

calc3<-tapply(predict(model3,type="response"), germination$extract,mean)
calc3

##      bean  cucumber 
## 0.3746835 0.6330275

p <- germination$count/germination$sample

tapply(p,germination$extract,mean)

##      bean  cucumber 
## 0.3487189 0.6031824

#Este resultado a pesar de ser cercano no es adecuado, primero hay que obtener el total de conteos

tapply(germination$count,germination$extract,sum) #numero total de germinación

##     bean cucumber 
##      148      276

tapply(germination$sample,germination$extract,sum) #numero total de semillas utilizadas

##     bean cucumber 
##      395      436

as.vector(tapply(germination$count,germination$extract,sum))/as.vector(tapply(germination$sample,germination$extract,sum))

## [1] 0.3746835 0.6330275

props <- read.table("flowering.txt",header=T)
attach(props)

## The following object is masked from data:
## 
##     dose

names(props)

## [1] "flowered" "number"   "dose"     "variety"

y <- cbind(flowered,number-flowered)
pf <- flowered/number
pfc <- split(pf,variety)
dc <- split(dose,variety)
plot(dose,pf,type="n",ylab="Proportion flowered")
points(jitter(dc[[1]]),jitter(pfc[[1]]),pch=21,col="blue",bg="red")
points(jitter(dc[[2]]),jitter(pfc[[2]]),pch=21,col="blue",bg="green")
points(jitter(dc[[3]]),jitter(pfc[[3]]),pch=21,col="blue",bg="yellow")
points(jitter(dc[[4]]),jitter(pfc[[4]]),pch=21,col="blue",bg="green3")
points(jitter(dc[[5]]),jitter(pfc[[5]]),pch=21,col="blue",bg="brown")

model1 <- glm(y~dose*variety,binomial)
summary(model1)

## 
## Call:
## glm(formula = y ~ dose * variety, family = binomial)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.6648  -1.1200  -0.3769   0.5735   3.3299  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -4.59165    1.03215  -4.449 8.64e-06 ***
## dose           0.41262    0.10033   4.113 3.91e-05 ***
## varietyB       3.06197    1.09317   2.801 0.005094 ** 
## varietyC       1.23248    1.18812   1.037 0.299576    
## varietyD       3.17506    1.07516   2.953 0.003146 ** 
## varietyE      -0.71466    1.54849  -0.462 0.644426    
## dose:varietyB -0.34282    0.10239  -3.348 0.000813 ***
## dose:varietyC -0.23039    0.10698  -2.154 0.031274 *  
## dose:varietyD -0.30481    0.10257  -2.972 0.002961 ** 
## dose:varietyE -0.00649    0.13292  -0.049 0.961057    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 303.350  on 29  degrees of freedom
## Residual deviance:  51.083  on 20  degrees of freedom
## AIC: 123.55
## 
## Number of Fisher Scoring iterations: 5

xv <- seq(0,35,0.1)
vn <- rep("A",length(xv))
yv <- predict(model1,list(variety=factor(vn),dose=xv),type="response")
lines(xv,yv,col="red")
vn <- rep("B",length(xv))
yv <- predict(model1,list(variety=factor(vn),dose=xv),type="response")
lines(xv,yv,col="green")
vn <- rep("C",length(xv))
yv <- predict(model1,list(variety=factor(vn),dose=xv),type="response")
lines(xv,yv,col="yellow")
vn <- rep("D",length(xv))
yv <- predict(model1,list(variety=factor(vn),dose=xv),type="response")
lines(xv,yv,col="green3")
vn <- rep("E",length(xv))
yv <- predict(model1,list(variety=factor(vn),dose=xv),type="response")
lines(xv,yv,col="brown")
legend(c(0.5,1),legend=c("A","B","C","D","E"),title="variety",
lty=rep(1,5),col=c("red","green","yellow","green3","brown"))

tapply(pf,list(dose,variety),mean)

##            A          B          C          D         E
## 0  0.0000000 0.08333333 0.00000000 0.06666667 0.0000000
## 1  0.0000000 0.00000000 0.14285714 0.11111111 0.0000000
## 4  0.0000000 0.20000000 0.06666667 0.15789474 0.0000000
## 8  0.4000000 0.50000000 0.17647059 0.53571429 0.1578947
## 16 0.8181818 0.90000000 0.25000000 0.73076923 0.7500000
## 32 1.0000000 0.50000000 1.00000000 0.77777778 1.0000000

# el modelo logístico se ajusta solamente a dos de las variedades (A y D), ajusta medianamente para una (C) y ajusta muy pobremente para las últimas dos (B y E)

#Se evidencia en la tabla y en la gráfica.

lizards <- read.table("lizards.txt",header=T)
attach(lizards)
names(lizards)

## [1] "n"       "sun"     "height"  "perch"   "time"    "species"

head(lizards)

##    n   sun height  perch    time  species
## 1 20 Shade   High  Broad Morning opalinus
## 2 13 Shade    Low  Broad Morning opalinus
## 3  8 Shade   High Narrow Morning opalinus
## 4  6 Shade    Low Narrow Morning opalinus
## 5 34   Sun   High  Broad Morning opalinus
## 6 31   Sun    Low  Broad Morning opalinus

sorted <- lizards[order(species,sun,height,perch,time),]
head(sorted)

##    n   sun height  perch      time  species
## 41 4 Shade   High  Broad Afternoon grahamii
## 33 1 Shade   High  Broad   Mid.day grahamii
## 25 2 Shade   High  Broad   Morning grahamii
## 43 3 Shade   High Narrow Afternoon grahamii
## 35 1 Shade   High Narrow   Mid.day grahamii
## 27 3 Shade   High Narrow   Morning grahamii

short <- sorted[1:24,]


names(short)[1] <- "Ag"
names(short)

## [1] "Ag"      "sun"     "height"  "perch"   "time"    "species"

short <- short[,-6]
head(short)

##    Ag   sun height  perch      time
## 41  4 Shade   High  Broad Afternoon
## 33  1 Shade   High  Broad   Mid.day
## 25  2 Shade   High  Broad   Morning
## 43  3 Shade   High Narrow Afternoon
## 35  1 Shade   High Narrow   Mid.day
## 27  3 Shade   High Narrow   Morning

sorted$n[25:48]

##  [1]  4  8 20  5  4  8 12  8 13  1  0  6 18 69 34  8 60 17 13 55 31  4 21 12

new.lizards <- data.frame(sorted$n[25:48], short)
names(new.lizards)[1] <- "Ao"
head(new.lizards)

##    Ao Ag   sun height  perch      time
## 41  4  4 Shade   High  Broad Afternoon
## 33  8  1 Shade   High  Broad   Mid.day
## 25 20  2 Shade   High  Broad   Morning
## 43  5  3 Shade   High Narrow Afternoon
## 35  4  1 Shade   High Narrow   Mid.day
## 27  8  3 Shade   High Narrow   Morning

detach(lizards)
rm(short,sorted)
attach(new.lizards)

names(new.lizards)

## [1] "Ao"     "Ag"     "sun"    "height" "perch"  "time"

y <- cbind(Ao,Ag)

model1 <- glm(y~sun*height*perch*time,binomial)

model2 <- step(model1)

## Start:  AIC=102.82
## y ~ sun * height * perch * time

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

##                         Df   Deviance    AIC
## - sun:height:perch:time  1 2.1801e-10 100.82
## <none>                     3.5823e-10 102.82

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## 
## Step:  AIC=100.82
## y ~ sun + height + perch + time + sun:height + sun:perch + height:perch + 
##     sun:time + height:time + perch:time + sun:height:perch + 
##     sun:height:time + sun:perch:time + height:perch:time

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

##                     Df Deviance     AIC
## - sun:height:time    2   0.4416  97.266
## - sun:perch:time     2   0.8101  97.634
## - height:perch:time  2   3.2217 100.046
## <none>                   0.0000 100.824
## - sun:height:perch   1   2.7088 101.533
## 
## Step:  AIC=97.27
## y ~ sun + height + perch + time + sun:height + sun:perch + height:perch + 
##     sun:time + height:time + perch:time + sun:height:perch + 
##     sun:perch:time + height:perch:time
## 
##                     Df Deviance    AIC
## - sun:perch:time     2   1.0713 93.896
## <none>                   0.4416 97.266
## - height:perch:time  2   4.6476 97.472
## - sun:height:perch   1   3.1113 97.936
## 
## Step:  AIC=93.9
## y ~ sun + height + perch + time + sun:height + sun:perch + height:perch + 
##     sun:time + height:time + perch:time + sun:height:perch + 
##     height:perch:time
## 
##                     Df Deviance    AIC
## - sun:time           2   3.3403 92.165
## <none>                   1.0713 93.896
## - sun:height:perch   1   3.3016 94.126
## - height:perch:time  2   5.7906 94.615
## 
## Step:  AIC=92.16
## y ~ sun + height + perch + time + sun:height + sun:perch + height:perch + 
##     height:time + perch:time + sun:height:perch + height:perch:time
## 
##                     Df Deviance    AIC
## <none>                   3.3403 92.165
## - sun:height:perch   1   5.8273 92.651
## - height:perch:time  2   8.5418 93.366

model3 <- update(model2,~. - height:perch:time)
model4 <- update(model2,~. - sun:height:perch)
anova(model2,model3,test="Chi")

## Analysis of Deviance Table
## 
## Model 1: y ~ sun + height + perch + time + sun:height + sun:perch + height:perch + 
##     height:time + perch:time + sun:height:perch + height:perch:time
## Model 2: y ~ sun + height + perch + time + sun:height + sun:perch + height:perch + 
##     height:time + perch:time + sun:height:perch
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
## 1         7     3.3403                       
## 2         9     8.5418 -2  -5.2014  0.07422 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(model2,model4,test="Chi")

## Analysis of Deviance Table
## 
## Model 1: y ~ sun + height + perch + time + sun:height + sun:perch + height:perch + 
##     height:time + perch:time + sun:height:perch + height:perch:time
## Model 2: y ~ sun + height + perch + time + sun:height + sun:perch + height:perch + 
##     height:time + perch:time + height:perch:time
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1         7     3.3403                     
## 2         8     5.8273 -1   -2.487   0.1148

model5 <- glm(y~(sun+height+perch+time)^2-sun:time,binomial)

model6 <- update(model5,~. - sun:height)
anova(model5,model6,test="Chi")

## Analysis of Deviance Table
## 
## Model 1: y ~ (sun + height + perch + time)^2 - sun:time
## Model 2: y ~ sun + height + perch + time + sun:perch + height:perch + 
##     height:time + perch:time
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1        10     10.903                     
## 2        11     13.254 -1  -2.3511   0.1252

model7 <- update(model5,~. - sun:perch)
anova(model5,model7,test="Chi")

## Analysis of Deviance Table
## 
## Model 1: y ~ (sun + height + perch + time)^2 - sun:time
## Model 2: y ~ sun + height + perch + time + sun:height + height:perch + 
##     height:time + perch:time
##   Resid. Df Resid. Dev Df  Deviance Pr(>Chi)
## 1        10     10.903                      
## 2        11     10.927 -1 -0.023597   0.8779

model8 <- update(model5,~. - height:perch)
anova(model5,model8,test="Chi")

## Analysis of Deviance Table
## 
## Model 1: y ~ (sun + height + perch + time)^2 - sun:time
## Model 2: y ~ sun + height + perch + time + sun:height + sun:perch + height:time + 
##     perch:time
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1        10     10.903                     
## 2        11     11.143 -1 -0.24006   0.6242

model9 <- update(model5,~. - time:perch)
anova(model5,model9,test="Chi")

## Analysis of Deviance Table
## 
## Model 1: y ~ (sun + height + perch + time)^2 - sun:time
## Model 2: y ~ sun + height + perch + time + sun:height + sun:perch + height:perch + 
##     height:time
##   Resid. Df Resid. Dev Df   Deviance Pr(>Chi)
## 1        10     10.903                       
## 2        12     10.909 -2 -0.0058263   0.9971

model10 <- update(model5,~. - time:height)
anova(model5,model10,test="Chi")

## Analysis of Deviance Table
## 
## Model 1: y ~ (sun + height + perch + time)^2 - sun:time
## Model 2: y ~ sun + height + perch + time + sun:height + sun:perch + height:perch + 
##     perch:time
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1        10     10.903                     
## 2        12     11.760 -2 -0.85679   0.6516

model11 <- glm(y~sun+height+perch+time,binomial)
summary(model11)

## 
## Call:
## glm(formula = y ~ sun + height + perch + time, family = binomial)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.66015  -0.37800   0.04488   0.62644   1.48717  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.2079     0.3536   3.416 0.000634 ***
## sunSun       -0.8473     0.3224  -2.628 0.008585 ** 
## heightLow     1.1300     0.2571   4.395 1.11e-05 ***
## perchNarrow  -0.7626     0.2113  -3.610 0.000306 ***
## timeMid.day   0.9639     0.2816   3.423 0.000619 ***
## timeMorning   0.7368     0.2990   2.464 0.013730 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 70.102  on 22  degrees of freedom
## Residual deviance: 14.205  on 17  degrees of freedom
## AIC: 83.029
## 
## Number of Fisher Scoring iterations: 4

timef <- factor(time,
                levels = c('Afternoon', 'Mid.day', 'Morning'))
t2 <- timef
levels(t2)[c(2,3)] <- "other"
levels(t2)

## [1] "Afternoon" "other"

model12 <- glm(y~sun+height+perch+t2,binomial)
anova(model11,model12,test="Chi")

## Analysis of Deviance Table
## 
## Model 1: y ~ sun + height + perch + time
## Model 2: y ~ sun + height + perch + t2
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1        17     14.205                     
## 2        18     15.023 -1 -0.81863   0.3656

summary(model12)

## 
## Call:
## glm(formula = y ~ sun + height + perch + t2, family = binomial)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.59707  -0.37407   0.06965   0.64616   1.53004  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.1595     0.3484   3.328 0.000874 ***
## sunSun       -0.7872     0.3159  -2.491 0.012722 *  
## heightLow     1.1188     0.2566   4.360  1.3e-05 ***
## perchNarrow  -0.7485     0.2104  -3.557 0.000375 ***
## t2other       0.8717     0.2611   3.338 0.000844 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 70.102  on 22  degrees of freedom
## Residual deviance: 15.023  on 18  degrees of freedom
## AIC: 81.847
## 
## Number of Fisher Scoring iterations: 4