-
Notifications
You must be signed in to change notification settings - Fork 0
/
.Rhistory
126 lines (126 loc) · 3.89 KB
/
.Rhistory
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
library(cenGAM)
library(cenGAM)
library(cenGAM)
library(cenGAM)
library(cenGAM)
?family.mgcv
set.seed(1)
x = matrix(2*rnorm(300), 100)
yn = 2*x[,3]+4*cos(x[,1]*2)
y = yn + rnorm(100)
ycensored = pmax(y, 0)
ycensored2 = replace(y, yn<0, 0); censid2 = which(yn<0)
sample(100, 30)-> censid3
ycensored3 = replace(y, censid3, 0);
par(mfrow = c(3,3))
plot(gam(y~s(x[,1])+s(x[,2]) + s(x[, 3])), ylim=c(-5, 5))
plot(gam(ycensored~s(x[,1])+s(x[,2]) + s(x[, 3])), ylim=c(-5, 5))
gam(cbind(ycensored, 0, Inf)~s(x[,1])+s(x[,2]) + s(x[, 3]), family = tobit1) -> obj
library(cenGAM) #this might not be compatible with some versions of mgcv, check with me if this doesn't work
library(cenGAM)
set.seed(1)
x = matrix(2*rnorm(300), 100)
yn = 2*x[,3]+4*cos(x[,1]*2)
y = yn + rnorm(100)
ycensored = pmax(y, 0)
ycensored2 = replace(y, yn<0, 0); censid2 = which(yn<0)
sample(100, 30)-> censid3
ycensored3 = replace(y, censid3, 0);
par(mfrow = c(3,3))
#true
set.seed(1)
x = matrix(2*rnorm(300), 100)
yn = 2*x[,3]+4*cos(x[,1]*2)
y = yn + rnorm(100)
ycensored = pmax(y, 0)
ycensored2 = replace(y, yn<0, 0); censid2 = which(yn<0)
sample(100, 30)-> censid3
ycensored3 = replace(y, censid3, 0);
par(mfrow = c(3,3))
#true
plot(gam(y~s(x[,1])+s(x[,2]) + s(x[, 3])), ylim=c(-5, 5))
plot(gam(ycensored~s(x[,1])+s(x[,2]) + s(x[, 3])), ylim=c(-5, 5))
gam(cbind(ycensored, 0, Inf)~s(x[,1])+s(x[,2]) + s(x[, 3]), family = tobit1) -> obj
summary(obj)
obj$family$thetaest
plot(obj, ylim = c(-5, 5))
devtools::document()
devtools::use_vignette("Fitting generalized additive models to censored response variables")
m <- gam(cbind(ycensored, 0, Inf) ~ s(x[,1]) + s(x[,2]) + s(x[, 3]), family = tobit1)
summary(m)
m$family$thetaest
x <- matrix(2*rnorm(300), 100)
yn <- 2*x[,3] + 4*cos(x[,1]*2)
y <- yn + rnorm(100)
ycensored <- pmax(y, 0)
par(mfrow = c(3,3))
# True model
plot(gam(y ~ s(x[,1]) + s(x[,2]) + s(x[, 3])), ylim=c(-5, 5))
# Naive estimation
plot(gam(ycensored ~ s(x[,1]) + s(x[,2]) + s(x[, 3])), ylim=c(-5, 5))
# Tobit I estimation
m <- gam(cbind(ycensored, 0, Inf) ~ s(x[,1]) + s(x[,2]) + s(x[, 3]), family = tobit1)
summary(m)
m$family$thetaest
plot(m, ylim = c(-5, 5))
plot(m)
plot(m)
plot(m, ylim = c(-5, 5),main="Tobit I")
set.seed(1)
x <- matrix(2*rnorm(300), 100)
yn <- 2*x[,3] + 4*cos(x[,1]*2)
y <- yn + rnorm(100)
ycensored <- pmax(y, 0)
par(mfrow = c(3,3))
# True model
plot(gam(y ~ s(x[,1]) + s(x[,2]) + s(x[, 3])), ylim=c(-5, 5), main = "True")
# Naive estimation
plot(gam(ycensored ~ s(x[,1]) + s(x[,2]) + s(x[, 3])), ylim=c(-5, 5), main = "Naive")
# Tobit I estimation
m <- gam(cbind(ycensored, 0, Inf) ~ s(x[,1]) + s(x[,2]) + s(x[, 3]), family = tobit1)
summary(m)
m$family$thetaest # TO REVISE: THIS IS CALLED SIGMA IN SUMMARY OUTPUT HEADING
plot(m, ylim = c(-5, 5), main = "Tobit I") #
head(x)
?pmax
yn
ycensored
?pmax
min(5:1, pi)
pmin(5:1, pi)
pmin(1:5, pi)
set.seed(1)
x <- matrix(2*rnorm(300), 100)
yn <- 2*x[,3] + 4*cos(x[,1]*2)
y <- yn + rnorm(100)
y
ycensored <- pmax(y, 0)
ycensored
ycensored2 <- replace(y, yn < 0, 0)
ycensored2
all(ycensored==ycensored2)
all(round(ycensored,5)==round(ycensored2,5))
round(ycensored,5)==round(ycensored2,5)
?replace
sample(100, 30)
ycensored
x <- matrix(2 * rnorm(300), 100)
yn <- 2 * x[, 3] + 4 * cos(x[, 1] * 2)
y <- yn + rnorm(100)
ycensored <- pmax(y, 0) # data left-censored at zero
ycensored2 <- replace(y, yn < 0, 0) # WHAT IS THE DIFFERENCE WITH ycensored?
censid2 <- which(yn < 0)
censid3 <- sample(100, 30)
ycensored3 <- replace(y, censid3, 0)
# Fitting alternative models
m1 <- gam(cbind(ycensored, 0, Inf) ~ s(x[, 1]) + s(x[, 2]) + s(x[, 3]), family = tobit2)
plot(m1, ylim = c(-5, 5))
m2 <- gam(cbind(ycensored2, replace(rep(-Inf, 100), censid2, 0), Inf) ~ s(x[, 1]) + s(x[, 2]) + s(x[, 3]), family = tobit2)
plot(m2, ylim = c(-5, 5))
m3 <- gam(cbind(ycensored3, replace(rep(0, 100), censid3,-1)) ~ s(x[, 1]) + s(x[, 2]) + s(x[, 3]), family = tobit2)
plot(m3, ylim = c(-5, 5))
m1
library(cenGAM)
library(cenGAM)
library(mgcv)
?family.mgcv