Lista 2 - Regressão 2

Autor

Paulo Manoel da Silva Junior

Regressão 2 - Segunda Lista de Exercícios

Resolução da Segunda lista de Exercícios da disciplina Regressão 2.

Aluno: Paulo Manoel da Silva Junior

Matrícula: 20190041314

Questão 1

O conjunto de dados descrito no arquivo brown.txt (Brown, 1980) apresenta as variáveis cancer, se o paciente possui ou não cancer de próstata (variável resposta), ap, nível de phosphatase ácida no sangue, xray, se há presença de manchas no raio-x, estagio, estágio do câncer, grau, grau do cancer e idade, idade do paciente.

setwd("C:/Users/Pessoal/Desktop/ESTATÍSTICA/UFPB/8º PERÍODO/REGRESSÃO 2/LISTA 2")

brow <- read.table("brown.txt", header = F, col.names = c("cancer","ap","xray","estagio","grau","idade"), sep = ",")

Visualizando algumas linhas do banco de dados

glimpse(brow)

Rows: 53
Columns: 6
$ cancer  <int> 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,…
$ ap      <int> 48, 56, 50, 52, 50, 49, 46, 62, 56, 55, 62, 71, 65, 67, 47, 49…
$ xray    <int> 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,…
$ estagio <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ grau    <int> 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0,…
$ idade   <int> 66, 68, 66, 56, 58, 60, 65, 60, 50, 49, 61, 58, 51, 67, 67, 51…

brow$cancer <- factor(brow$cancer, levels = c(0,1), labels = c("Não", "Sim"))
brow$xray <- as.factor(brow$xray)
brow$estagio <- as.factor(brow$estagio)
brow$grau <- as.factor(brow$grau)

glimpse(brow)

Rows: 53
Columns: 6
$ cancer  <fct> Não, Não, Não, Não, Não, Não, Não, Não, Sim, Não, Não, Não, Nã…
$ ap      <int> 48, 56, 50, 52, 50, 49, 46, 62, 56, 55, 62, 71, 65, 67, 47, 49…
$ xray    <fct> 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,…
$ estagio <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ grau    <fct> 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0,…
$ idade   <int> 66, 68, 66, 56, 58, 60, 65, 60, 50, 49, 61, 58, 51, 67, 67, 51…

brow %>%
  slice_head(n=15) %>%
  gt::gt()

cancer	ap	xray	grau	idade
Não	48	0	0	66
Não	56	0	0	68
Não	50	0	0	66
Não	52	0	0	56
Não	50	0	0	58
Não	49	0	0	60
Não	46	1	0	65
Não	62	1	0	60
Sim	56	0	1	50
Não	55	1	0	49
Não	62	0	0	61
Não	71	0	0	58
Não	65	0	0	51
Sim	67	1	1	67
Não	47	0	1	67

Estatística Descritiva

skimr::skim(brow)

Data summary
Name	brow
Number of rows	53
Number of columns	6
_______________________
Column type frequency:
factor	4
numeric	2
________________________
Group variables	None

Variable type: factor

skim_variable	complete_rate	ordered	n_unique	top_counts
cancer	1	FALSE	2	Não: 33, Sim: 20
xray	1	FALSE	2	0: 38, 1: 15
estagio	1	FALSE	2	1: 27, 0: 26
grau	1	FALSE	2	0: 32, 1: 21

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
ap	0	1	69.42	26.20	40	50	65	78	187	▇▅▁▁▁
idade	0	1	59.38	6.17	45	56	60	65	68	▁▃▅▆▇

Letra A

Realize o ajuste da regressão logística e selecione as variáveis. O modelo é adequado? Justifique sua escolha.

fit1 <- glm(cancer~xray+idade+ap+estagio+grau, data = brow, family = binomial(link = "logit"))

Verificando se o modelo é adequado

phi1 <- summary(fit1)$dispersion
desvio1 <- summary(fit1)$deviance/phi1
q.quadro1 <- qchisq(0.95,desvio1)
ifelse(desvio1<q.quadro1, TRUE, FALSE)

[1] TRUE

Como podemos ver o modelo é adequado, prosseguiremos assim para a realização da seleção de variáveis, através do método stepwise.

menor <- glm(cancer~1, data = brow, family = binomial(link = "logit"))
maior <- glm(cancer~xray+ap+idade+grau+estagio, family = binomial(link="logit"), data = brow)
step(menor, scope = list(upper = maior), direction = "both")

Start:  AIC=72.25
cancer ~ 1

          Df Deviance    AIC
+ xray     1   59.001 63.001
+ estagio  1   62.553 66.553
+ grau     1   67.089 71.089
+ ap       1   67.116 71.116
<none>         70.252 72.252
+ idade    1   69.162 73.162

Step:  AIC=63
cancer ~ xray

          Df Deviance    AIC
+ estagio  1   53.353 59.353
<none>         59.001 63.001
+ ap       1   57.059 63.059
+ grau     1   57.186 63.186
+ idade    1   57.660 63.660
- xray     1   70.252 72.252

Step:  AIC=59.35
cancer ~ xray + estagio

          Df Deviance    AIC
+ ap       1   50.660 58.660
<none>         53.353 59.353
+ idade    1   52.085 60.085
+ grau     1   53.126 61.126
- estagio  1   59.001 63.001
- xray     1   62.553 66.553

Step:  AIC=58.66
cancer ~ xray + estagio + ap

          Df Deviance    AIC
<none>         50.660 58.660
+ idade    1   49.097 59.097
- ap       1   53.353 59.353
+ grau     1   50.183 60.183
- estagio  1   57.059 63.059
- xray     1   58.613 64.613


Call:  glm(formula = cancer ~ xray + estagio + ap, family = binomial(link = "logit"), 
    data = brow)

Coefficients:
(Intercept)        xray1     estagio1           ap  
   -3.57565      2.06179      1.75556      0.02063  

Degrees of Freedom: 52 Total (i.e. Null);  49 Residual
Null Deviance:      70.25 
Residual Deviance: 50.66    AIC: 58.66

Resposta: Podemos observar que a seleção de variáveis, resultou de que o melhor modelo tem como variáveis independentes xray, estágio e ap.

Verificando a adequação do modelo final

fit2 <- glm(cancer ~ xray + estagio + ap, family = binomial(link = "logit"),data = brow)
phi11 <- summary(fit2)$dispersion
desvio11 <- summary(fit2)$deviance/phi11
q.quadr11 <- qchisq(0.95, desvio11)

	Modelo
Desvio/\phi	50.6596343
\chi^2	68.2731267

Temos do resultado acima que: O modelo em análise foi aceito.

ajuste<-summary(fit2)
knitr::kable(ajuste$coefficients)

	Estimate	Std. Error	z value	Pr(>\|z\|)
(Intercept)	-3.5756480	1.1811454	-3.027272	0.0024677
xray1	2.0617874	0.7776701	2.651236	0.0080198
estagio1	1.7555616	0.7390242	2.375513	0.0175246
ap	0.0206294	0.0126490	1.630918	0.1029076

Podemos observar pelo teste aplicado que existe um forte indício de que a variável ap não seja tão relevante para o ajuste que foi realizado, pois a hipótese de que \beta = 0 não foi rejeitada. Todavia, manteremos ela no ajuste que foi selecionado.

Letra B

Construa um envelope para os resíduos. Há algum ponto que não pertence ao envelope? Se sim, qual?

par(mfrow=c(1,1))
X <- model.matrix(fit2)
n <- nrow(X)
p <- ncol(X)
w <- fit2$weights
W <- diag(w)
H <- solve(t(X)%*%W%*%X)
H <- sqrt(W)%*%X%*%H%*%t(X)%*%sqrt(W)
h <- diag(H)
td <- resid(fit2,type="deviance")/sqrt(1-h)
e <- matrix(0,n,100)
#
for(i in 1:100){
dif <- runif(n) - fitted(fit2)
dif[dif >= 0 ] <- 0
dif[dif<0] <- 1
nresp <- dif
fit <- glm(nresp ~ X, family=binomial)
w <- fit$weights
W <- diag(w)
H <- solve(t(X)%*%W%*%X)
H <- sqrt(W)%*%X%*%H%*%t(X)%*%sqrt(W)
h <- diag(H)
e[,i] <- sort(resid(fit,type="deviance")/sqrt(1-h))}
#
e1 <- numeric(n)
e2 <- numeric(n)
#
for(i in 1:n){
  eo <- sort(e[i,])
e1[i] <- (eo[2]+eo[3])/2
e2[i] <- (eo[97]+eo[98])/2}
#
med <- apply(e,1,mean)
faixa <- range(td,e1,e2)
par(pty="s")
qqnorm(td,xlab="Percentil da N(0,1)",
ylab="Componente do Desvio", ylim=faixa, pch=16)
#
par(new=T)
#
qqnorm(e1,axes=F,xlab="",ylab="",type="l",ylim=faixa,lty=1)
par(new=T)
qqnorm(e2,axes=F,xlab="",ylab="", type="l",ylim=faixa,lty=1)
par(new=T)
qqnorm(med,axes=F,xlab="", ylab="", type="l",ylim=faixa,lty=2)

Resposta: Não há pontos fora do pacote.

Letra C

Construa um intervalo de confiança de 98% para os parâmetros do modelo.

confint(fit2, level = 0.98) %>%
  knitr::kable()

Waiting for profiling to be done...

	1 %	99 %
(Intercept)	-6.7217321	-1.0541869
xray1	0.3512182	4.0831080
estagio1	0.1363791	3.6816376
ap	-0.0091789	0.0551541

Questão 2

O conjunto de dados descrito no arquivo MissAmerica2008.txt apresenta as variáveis Y, número de vezes que cada estado americano teve alguma finalista dos anos de 2000 a 2008 (variável resposta), x1, logaritmo da populalação, x2, logaritmo da média final de qualificação de concorrentes de cada estado, x3, logaritmo da área geográfica de cada estatdo, x4, latitude de cada estado e x5, longitude de cada estado. (Até o ano de 2008, ocorreram 81 concursos de beleza).

Miss <- read.table("MissAmerica2008.txt", dec = ".", header = F, col.names = c("Estado", "Vezes", "população", "média_final", "área_geográfica", "latitude", "longitude"))

glimpse(Miss)

Rows: 51
Columns: 7
$ Estado          <chr> "AL", "AK", "AZ", "AR", "CA", "CO", "CT", "DE", "DC", …
$ Vezes           <int> 6, 0, 0, 4, 5, 0, 2, 0, 2, 3, 4, 3, 0, 2, 2, 1, 1, 2, …
$ população       <dbl> 11.9249, 9.8011, 12.0543, 11.2702, 14.0005, 11.8820, 1…
$ média_final     <dbl> 3.895894, 2.708050, 2.862201, 3.766997, 3.935740, 2.94…
$ área_geográfica <dbl> 10.8670, 13.4049, 11.6439, 10.8814, 12.0058, 11.5530, …
$ latitude        <dbl> 32.3833, 58.3667, 33.4333, 34.7333, 38.5167, 39.7500, …
$ longitude       <dbl> 86.367, 134.583, 112.017, 92.233, 121.500, 104.867, 72…

Retirando a primeira coluna do banco de dados

Miss <- Miss[,-1]

Visualizando algumas linhas do banco de dados

Miss %>%
  slice_head(n=10) %>%
  gt::gt()

Vezes	população	média_final	área_geográfica	latitude	longitude
6	11.9249	3.895894	10.8670	32.3833	86.367
0	9.8011	2.708050	13.4049	58.3667	134.583
0	12.0543	2.862201	11.6439	33.4333	112.017
4	11.2702	3.766997	10.8814	34.7333	92.233
5	14.0005	3.935740	12.0058	38.5167	121.500
0	11.8820	2.944439	11.5530	39.7500	104.867
2	11.5742	2.944439	8.6203	41.7333	72.650
0	10.3397	2.852631	7.8196	39.1333	75.467
2	10.3899	2.456736	4.2195	38.8500	77.033
3	13.0882	3.725693	11.0937	30.3833	84.367

y <- 81-Miss$Vezes
Miss <- data.frame(y,Miss)

Miss%>%
  slice_head(n=10)%>%
  gt::gt()

y	Vezes	população	média_final	área_geográfica	latitude	longitude
75	6	11.9249	3.895894	10.8670	32.3833	86.367
81	0	9.8011	2.708050	13.4049	58.3667	134.583
81	0	12.0543	2.862201	11.6439	33.4333	112.017
77	4	11.2702	3.766997	10.8814	34.7333	92.233
76	5	14.0005	3.935740	12.0058	38.5167	121.500
81	0	11.8820	2.944439	11.5530	39.7500	104.867
79	2	11.5742	2.944439	8.6203	41.7333	72.650
81	0	10.3397	2.852631	7.8196	39.1333	75.467
79	2	10.3899	2.456736	4.2195	38.8500	77.033
78	3	13.0882	3.725693	11.0937	30.3833	84.367

Estatística Descritiva

skimr::skim(Miss)

Data summary
Name	Miss
Number of rows	51
Number of columns	7
_______________________
Column type frequency:
numeric	7
________________________
Group variables	None

Variable type: numeric

skim_variable	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
y	1	79.24	1.73	74.00	78.00	79.00	81.00	81.00	▁▁▃▃▇
Vezes	1	1.76	1.73	0.00	0.00	2.00	3.00	7.00	▇▃▃▁▁
população	1	11.68	1.03	9.69	10.84	11.76	12.35	14.00	▃▅▇▅▁
média_final	1	3.14	0.46	1.91	2.86	3.09	3.45	4.04	▁▃▇▅▅
área_geográfica	1	10.62	1.43	4.22	10.49	10.94	11.34	13.40	▁▁▂▇▂
latitude	1	39.40	5.67	21.33	35.99	39.75	42.77	58.37	▁▃▇▂▁
longitude	1	93.14	18.65	69.80	78.06	89.67	102.78	157.92	▇▅▂▁▁

Letra A

Realize o ajuste da regressão logística e selecione as variáveis. O modelo é adequado? Justifique sua escolha.

fit3 <- glm(cbind(y,Vezes)~+população+média_final+área_geográfica+latitude+longitude, data = Miss, family = binomial(link = "logit"))

Verificando se o modelo é adequado

phi21 <- summary(fit3)$dispersion
desvio21 <- summary(fit3)$deviance/phi21
q.quadro21 <- qchisq(0.95,desvio21)
ifelse(desvio21<q.quadro21, TRUE, FALSE)

[1] TRUE

Como podemos ver o modelo é adequado, prosseguiremos assim para a realização da seleção de variáveis, através do método stepwise.

menor1 <- glm(cbind(y,Vezes)~1, data = Miss, family = binomial(link = "logit"))
maior1 <- glm(cbind(y,Vezes)~população+média_final+área_geográfica+latitude+longitude, family = binomial(link="logit"), data = Miss)
step(menor1, scope = list(upper = maior1), direction = "both")

Start:  AIC=192.87
cbind(y, Vezes) ~ 1

                  Df Deviance    AIC
+ média_final      1   60.257 158.03
+ população        1   68.530 166.31
+ latitude         1   75.863 173.64
<none>                 97.095 192.87
+ longitude        1   96.148 193.92
+ área_geográfica  1   97.000 194.78

Step:  AIC=158.03
cbind(y, Vezes) ~ média_final

                  Df Deviance    AIC
+ população        1   54.051 153.83
+ latitude         1   54.762 154.54
+ área_geográfica  1   56.217 155.99
<none>                 60.257 158.03
+ longitude        1   59.076 158.85
- média_final      1   97.095 192.87

Step:  AIC=153.83
cbind(y, Vezes) ~ média_final + população

                  Df Deviance    AIC
+ área_geográfica  1   46.752 148.53
+ latitude         1   48.462 150.24
<none>                 54.051 153.83
+ longitude        1   53.353 155.13
- população        1   60.257 158.03
- média_final      1   68.530 166.31

Step:  AIC=148.53
cbind(y, Vezes) ~ média_final + população + área_geográfica

                  Df Deviance    AIC
+ latitude         1   42.367 146.14
<none>                 46.752 148.53
+ longitude        1   46.314 150.09
- área_geográfica  1   54.051 153.83
- população        1   56.217 155.99
- média_final      1   66.531 166.31

Step:  AIC=146.14
cbind(y, Vezes) ~ média_final + população + área_geográfica + 
    latitude

                  Df Deviance    AIC
<none>                 42.367 146.14
+ longitude        1   42.163 147.94
- latitude         1   46.752 148.53
- área_geográfica  1   48.462 150.24
- média_final      1   51.353 153.13
- população        1   52.095 153.87


Call:  glm(formula = cbind(y, Vezes) ~ média_final + população + 
    área_geográfica + latitude, family = binomial(link = "logit"), 
    data = Miss)

Coefficients:
    (Intercept)      média_final        população  área_geográfica  
        7.74751         -1.06493         -0.45263          0.27347  
       latitude  
        0.05503  

Degrees of Freedom: 50 Total (i.e. Null);  46 Residual
Null Deviance:      97.09 
Residual Deviance: 42.37    AIC: 146.1

Resultado: Depois de aplicada a seleção stepWise, ele verificou que apenas as variáveis média_final, população, área_geográfica e latitude são suficientes para explicar o modelo.

fit4 <- glm(cbind(y,Vezes) ~ média_final + população + área_geográfica+latitude, family = binomial(link = "logit"),data = Miss)
phi22 <- summary(fit4)$dispersion
desvio22 <- summary(fit4)$deviance/phi22
q.quadr22 <- qchisq(0.95, desvio22)

	Modelo
Desvio/\phi	42.3671014
\chi^2	58.557273

Temos do resultado acima que: O modelo em análise foi aceito.

Letra B

Identifique se há algum ponto de alavanca.

h1<- boot::glm.diag(fit4)$h
indice <- as.vector(1:length(Miss$Vezes))
h3 <- ggplot2::ggplot(Miss, ggplot2::aes(x=indice, y=h1)) + ggplot2::geom_point(size=4, color='#e7ad52') + ggplot2::geom_abline(intercept = 2*(length(coef(fit4))/length(Miss$Vezes)), slope = 0, color='#001132') + ggplot2::theme(axis.title.x = ggplot2::element_text(colour = '#001132'), axis.text.x = ggplot2::element_text(colour = '#001132'), axis.title.y = ggplot2::element_text(colour = '#001132'), axis.text.y = ggplot2::element_text(colour = '#001132')) 
h3

Conforme o gráfico vemos que tem alguns pontos de alavanca, mais precisamente 6 pontos. Que são as observações de número:5, 9, 10, 12, 44 e 45.

Miss[which(h1>2*length(coef(fit4))/length(Miss$Vezes)),]

    y Vezes população média_final área_geográfica latitude longitude
5  76     5   14.0005    3.935740         12.0058  38.5167   121.500
9  79     2   10.3899    2.456736          4.2195  38.8500    77.033
10 78     3   13.0882    3.725693         11.0937  30.3833    84.367
12 78     3   10.6205    2.564949          9.2994  21.3333   157.917
44 74     7   13.4640    3.811097         12.5009  30.3000    97.700
45 79     2   11.5123    4.040123         11.3492  40.7667   111.967

Letra C

Identifique se há algum ponto influente.

DC<-boot::glm.diag(fit4)$cook
dc <- ggplot2::ggplot(Miss, ggplot2::aes(x=indice, y=DC)) + ggplot2::geom_point(size=4, color='#e7ad52') + ggplot2::geom_abline(intercept = qchisq(0.05,length(coef(fit4)))/5, slope = 0, color='#001132') + ggplot2::theme(axis.title.x = ggplot2::element_text(colour = '#001132'), axis.text.x = ggplot2::element_text(colour = '#001132'), axis.title.y = ggplot2::element_text(colour = '#001132'), axis.text.y = ggplot2::element_text(colour = '#001132')) 
dc

É observado que temos um único ponto de influência, que é a observação de número: 12. Sendo os valores observados:

Miss[which(DC>qchisq(0.05,length(coef(fit4)))/5),]

    y Vezes população média_final área_geográfica latitude longitude
12 78     3   10.6205    2.564949          9.2994  21.3333   157.917

Questão 3

Milecer e Szczotka (1966) investigam a idade do início da menstruação em 3918 garotas de Varsóvia. Para 25 médias de idade foram observadas a ocorrência (Y = 1) ou não (Y = 0) do início de períodos de menstruação nas adolescentes. O conjunto de dados descrito no arquivo meninas.txt apresenta as variáveis y1, garotas menstruando, m, garotas entrevistadas e idade, idade média.

meninas <- read.table("meninas.txt", header = F, col.names = c("Garotas_mestruando", "Garotas_entrevistadas","Idade"), dec = ".")

glimpse(meninas)

Rows: 25
Columns: 3
$ Garotas_mestruando    <int> 0, 0, 0, 2, 2, 5, 10, 17, 16, 29, 39, 51, 47, 67…
$ Garotas_entrevistadas <int> 376, 200, 93, 120, 90, 88, 105, 111, 100, 93, 10…
$ Idade                 <dbl> 9.21, 10.21, 10.58, 10.83, 11.08, 11.33, 11.58, …

Visualizando as primeiras linhas do banco de dados

meninas %>%
  slice_head(n=10) %>%
  gt::gt()

Garotas_mestruando	Garotas_entrevistadas	Idade
0	376	9.21
0	200	10.21
0	93	10.58
2	120	10.83
2	90	11.08
5	88	11.33
10	105	11.58
17	111	11.83
16	100	12.08
29	93	12.33

Estatística Descritiva

skimr::skim(meninas)

Data summary
Name	meninas
Number of rows	25
Number of columns	3
_______________________
Column type frequency:
numeric	3
________________________
Group variables	None

Variable type: numeric

skim_variable	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
Garotas_mestruando	1	92.32	203.74	0.00	10.00	51.00	92.00	1049.00	▇▁▁▁▁
Garotas_entrevistadas	1	156.72	194.60	88.00	98.00	105.00	117.00	1049.00	▇▁▁▁▁
Idade	1	13.10	2.03	9.21	11.58	13.08	14.58	17.53	▅▇▇▇▁

Letra A

Realize o ajuste da regressão logística. O modelo é adequado? Justifique sua escolha.

y <- meninas$Garotas_entrevistadas-meninas$Garotas_mestruando
meninas <- data.frame(y,meninas)
fit <- glm(cbind(y,Garotas_mestruando)~Idade, family = binomial(link = "logit"), data = meninas)

phi22 <- summary(fit)$dispersion
desvio22 <- summary(fit)$deviance/phi22
q.quadr22 <- qchisq(0.95, desvio22)

	Modelo
Desvio/\phi	26.8055526
\chi^2	39.8748065

Temos do resultado acima que: O modelo em análise foi aceito.

Letra B

Faça um gráfico dos valores ajustados com a idade média.

plot(fitted(fit)~meninas$Idade)

Letra C

Qual a chance de uma menina menstruar se sua idade é de 13 anos?

predict(fit,list(Idade=13),type="response")

        1 
0.5026559

chance<-0.4973441/(1-0.4973441)
chance

[1] 0.9894325

Questão 4

Considere o banco de dados trees do R que fornece a medida da circunferência (Girth), altura (Height) e volume (Volume - variável resposta) da madeira de 31 cerejeiras pretas derrubadas.

Conhecendo o banco de dados

glimpse(trees)

Rows: 31
Columns: 3
$ Girth  <dbl> 8.3, 8.6, 8.8, 10.5, 10.7, 10.8, 11.0, 11.0, 11.1, 11.2, 11.3, …
$ Height <dbl> 70, 65, 63, 72, 81, 83, 66, 75, 80, 75, 79, 76, 76, 69, 75, 74,…
$ Volume <dbl> 10.3, 10.3, 10.2, 16.4, 18.8, 19.7, 15.6, 18.2, 22.6, 19.9, 24.…

Visualizando algumas linhas do banco de dados

trees %>%
  slice_head(n=10) %>%
  gt::gt()

Girth	Height	Volume
8.3	70	10.3
8.6	65	10.3
8.8	63	10.2
10.5	72	16.4
10.7	81	18.8
10.8	83	19.7
11.0	66	15.6
11.0	75	18.2
11.1	80	22.6
11.2	75	19.9

Estatística Descritiva

skimr::skim(trees)

Data summary
Name	trees
Number of rows	31
Number of columns	3
_______________________
Column type frequency:
numeric	3
________________________
Group variables	None

Variable type: numeric

skim_variable	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
Girth	1	13.25	3.14	8.3	11.05	12.9	15.25	20.6	▃▇▃▅▁
Height	1	76.00	6.37	63.0	72.00	76.0	80.00	87.0	▃▃▆▇▃
Volume	1	30.17	16.44	10.2	19.40	24.2	37.30	77.0	▇▅▁▂▁

Letra A

Faça o gráfico de dispersão da variável resposta Volume pelas variávei explicativas Girth e Height.

ggplot2::ggplot(trees, ggplot2::aes(x = Girth, y = Volume)) +
  ggplot2::geom_point()+
  ggplot2::xlab("Circunferência") + 
  ggplot2::ylab("Volume") +
  ggplot2::labs(title = "Gráfico de dispersão entre o Volume e a Circunferência", subtitle = "De 31 cerejeiras pretas derrubadas.")

ggplot2::ggplot(trees, ggplot2::aes(x = Height, y = Volume)) +
  ggplot2::geom_point()+
  ggplot2::xlab("Altura") + 
  ggplot2::ylab("Volume") +
  ggplot2::labs(title = "Gráfico de dispersão entre o Volume e a Altura", subtitle = "De 31 cerejeiras pretas derrubadas.")

Letra B

Realize o ajuste de um modelo GAM com a variável resposta Volume tendo uma distribuição Gama.

Verificando através de um histograma se a melhor distribuição é a Gama para a variável resposta

ggplot2::ggplot(trees, ggplot2::aes(Volume))+
  ggplot2::geom_histogram(fill='white', 
                          color = "black", 
                          breaks = hist(trees$Volume, plot = F)$breaks) +
  ggplot2::xlab("Volume")+
  ggplot2::ylab("Frequência")+
  ggplot2::labs(title = "Histograma do Volume", subtitle = "De 31 cerejeiras pretas derrubadas.")

É observado que podemos utilizar o ajuste de acordo com a distribuição Gama

Primeiro Ajuste

fit6 <- mgcv::gam(Volume~te(Girth)+te(Height), family = Gamma(), data = trees)
summary(fit6)


Family: Gamma 
Link function: inverse 

Formula:
Volume ~ te(Girth) + te(Height)

Parametric coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.0428631  0.0007279   58.89   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Approximate significance of smooth terms:
             edf Ref.df      F p-value    
te(Girth)  3.884  3.991 141.13 < 2e-16 ***
te(Height) 1.000  1.000  24.29 4.6e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.976   Deviance explained = 97.8%
GCV = 0.0088946  Scale est. = 0.0071462  n = 31

Segundo ajuste

fit7 <- mgcv::gam(Volume~t2(Girth)+t2(Height), family = Gamma(), data = trees)
summary(fit7)


Family: Gamma 
Link function: inverse 

Formula:
Volume ~ t2(Girth) + t2(Height)

Parametric coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.0428631  0.0007279   58.89   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Approximate significance of smooth terms:
             edf Ref.df      F p-value    
t2(Girth)  3.884  3.991 143.55 < 2e-16 ***
t2(Height) 1.000  1.000  24.29 4.6e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.976   Deviance explained = 97.8%
GCV = 0.0088946  Scale est. = 0.0071462  n = 31

Terceiro ajuste

fit68 <- mgcv::gam(Volume~s(Girth)+s(Height), family = Gamma(), data = trees)
summary(fit68)


Family: Gamma 
Link function: inverse 

Formula:
Volume ~ s(Girth) + s(Height)

Parametric coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.0428076  0.0007274   58.85   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Approximate significance of smooth terms:
            edf Ref.df     F  p-value    
s(Girth)  5.632  6.579 84.13  < 2e-16 ***
s(Height) 1.000  1.000 24.42 5.42e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.976   Deviance explained =   98%
GCV = 0.0096142  Scale est. = 0.0072122  n = 31

Escolha do melhor modelo

Resultado: Conforme os ajustes acima, o terceiro ajuste, no qual utiliza o parâmetro de suavização s, se mostrou mais eficiente, tendo um coeficiente R^2 ajustado de 97.6% e um Desvio explicado de cerca de 98%.

Sendo assim, escolhemos este o modelo a ser adotado, por ter evidenciado estar bem ajustado.

Letra C

Faça uma análise de diagnósticos do modelo escolhido. O que você pode concluir do modelo?

par(mfrow = c(1,2))
plot(fit68)

Através da análise gráfica acima, podemos observar que o modelo ficou bem ajustado, vamos observar outro gráfico para reforçar.

par(mfrow = c(2,2))
mgcv::gam.check(fit68)


Method: GCV   Optimizer: outer newton
full convergence after 5 iterations.
Gradient range [-7.44552e-09,-1.588236e-09]
(score 0.009614184 & scale 0.007212205).
Hessian positive definite, eigenvalue range [7.445419e-09,0.0004204206].
Model rank =  19 / 19 

Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.

            k'  edf k-index p-value
s(Girth)  9.00 5.63    1.25    0.90
s(Height) 9.00 1.00    1.34    0.94

Resposta: De acordo com as duas análises gráficas dos resíduos podemos verificar que o modelo ficou bem ajustado e que os desvios são aproximadamente normais.

Questão 5

Considere o banco de dados mcycle do pacote MASS do R que fornece uma série de medidas da aceleração da cabeça (accel - variável resposta) em um acidente de motocicleta simulado, usado para testar o capacete em acidentes.

Conhecendo o banco de dados

glimpse(MASS::mcycle)

Rows: 133
Columns: 2
$ times <dbl> 2.4, 2.6, 3.2, 3.6, 4.0, 6.2, 6.6, 6.8, 7.8, 8.2, 8.8, 8.8, 9.6,…
$ accel <dbl> 0.0, -1.3, -2.7, 0.0, -2.7, -2.7, -2.7, -1.3, -2.7, -2.7, -1.3, …

Visualizando algumas linhas do banco de dados

MASS::mcycle %>%
  slice_head(n=10) %>%
  gt::gt()

times	accel
2.4	0.0
2.6	-1.3
3.2	-2.7
3.6	0.0
4.0	-2.7
6.2	-2.7
6.6	-2.7
6.8	-1.3
7.8	-2.7
8.2	-2.7

Estatística Descritiva

skimr::skim(MASS::mcycle)

Data summary
Name	MASS::mcycle
Number of rows	133
Number of columns	2
_______________________
Column type frequency:
numeric	2
________________________
Group variables	None

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
times	0	1	25.18	13.13	2.4	15.6	23.4	34.8	57.6	▃▇▅▃▂
accel	0	1	-25.55	48.32	-134.0	-54.9	-13.3	0.0	75.0	▃▃▅▇▂

Letra A

Faça o gráfico de dispersão da variável resposta accel pela variável explicativa tempo (times)

ggplot2::ggplot(MASS::mcycle, ggplot2::aes(x = times, y = accel)) +
  ggplot2::geom_point()+
  ggplot2::xlab("Times") + 
  ggplot2::ylab("Accel") +
  ggplot2::labs(title = "Gráfico de dispersão entre Times e Accel", subtitle = "De 133 simulações.")

Letra B

1. Realize o ajuste de um modelo GAM com a variável resposta accel tendo uma distribuição Normal.

Verificando através de um histograma se a melhor distribuição é a Normal para a variável resposta

ggplot2::ggplot(MASS::mcycle, ggplot2::aes(accel))+
  ggplot2::geom_histogram(fill='white', 
                          color = "black", 
                          breaks = hist(MASS::mcycle$accel, plot = F)$breaks) +
  ggplot2::xlab("Accel") + 
  ggplot2::ylab("Frequência Absoluta") +
  ggplot2::labs(title = "Histograma do Accel", subtitle = "De 133 simulações.")

Podemos observar de maneira assintótica que os dados podem se comportar mais a frente de maneira simétrica.

Primeiro Ajuste

fit8 <- mgcv::gam(accel~te(times), data = MASS::mcycle)
summary(fit8)


Family: gaussian 
Link function: identity 

Formula:
accel ~ te(times)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -25.546      3.191  -8.006 6.19e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Approximate significance of smooth terms:
            edf Ref.df     F p-value    
te(times) 3.813  3.978 23.46  <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =   0.42   Deviance explained = 43.7%
GCV = 1404.9  Scale est. = 1354      n = 133

Segundo ajuste

fit9 <- mgcv::gam(accel~t2(times), data =MASS::mcycle)
summary(fit9)


Family: gaussian 
Link function: identity 

Formula:
accel ~ t2(times)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -25.546      3.191  -8.006 6.19e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Approximate significance of smooth terms:
            edf Ref.df     F p-value    
t2(times) 3.813  3.978 26.35  <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =   0.42   Deviance explained = 43.7%
GCV = 1404.9  Scale est. = 1354      n = 133

Terceiro ajuste

fit10 <- mgcv::gam(accel~s(times), data = MASS::mcycle)
summary(fit10)


Family: gaussian 
Link function: identity 

Formula:
accel ~ s(times)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -25.546      1.951   -13.1   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Approximate significance of smooth terms:
           edf Ref.df     F p-value    
s(times) 8.693  8.972 53.52  <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.783   Deviance explained = 79.8%
GCV = 545.78  Scale est. = 506       n = 133

Escolha do melhor modelo

Resultado: Conforme os ajustes acima, o terceiro ajuste, no qual utiliza o parâmetro de suavização s, se mostrou mais eficiente, tendo um coeficiente R^2 ajustado de 78.3% e um Desvio explicado de cerca de 79.8%.

Sendo assim, escolhemos este o modelo a ser adotado, por ter evidenciado estar bem ajustado.

Letra C

Faça uma análise de diagnósticos do modelo escolhido. O que você pode concluir do modelo?

par(mfrow = c(1,1))
plot(fit10)

Através da análise gráfica acima, podemos observar que o modelo ficou bem ajustado, vamos observar outro gráfico para reforçar.

par(mfrow = c(2,2))
mgcv::gam.check(fit10)


Method: GCV   Optimizer: magic
Smoothing parameter selection converged after 8 iterations.
The RMS GCV score gradient at convergence was 0.0009645413 .
The Hessian was positive definite.
Model rank =  10 / 10 

Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.

           k'  edf k-index p-value
s(times) 9.00 8.69    1.15    0.94

Resposta: De acordo com as duas análises gráficas dos resíduos podemos verificar que o modelo ficou bem ajustado e que os desvios são aproximadamente normais. Parece que teríamos nos dados problema com a heterocedasticidade.

Questão 6

Considere o banco de dados rent do pacote gamlss.data do R. A pesquisa foi realizada em abril de 1993 pela Infratest Sozialforschung. Uma amostra aleatória de alojamento com novos contratos de arrendamento ou aumentos de rendas nos últimos quatro anos em Munique foi selecionado incluindo: i) quartos individuais , ii) pequenos apartamentos , iii ) apartamentos, iv) casas de dois familiares. Alojamento sujeitas a rendas de preços de controle , casas de uma família e casas especiais, como coberturas, foram excluídos porque eles são um pouco diferente do resto e são considerados um mercado distinto. Para o propósito deste estudo, foram utilizadas 1967 observações das variáveis listadas abaixo, isto é, a variável R alugar variável resposta seguido por as variáveis explanatórias que se verificou serem apropriados para uma abordagem de análise de regressão por Fahrmeir et ai. (1994, 1995): Fl metro quadrado do espaço, A ano de construção, Sp indica se o local é acima da média 1 ou não 0, Sm indica se o local é abaixo da média 1 ou não 0, B indica se tem banheiro 1 ou não 0, H indica se tem aquecimento central 1 ou não 0, L indica se os equipamentos da cozinha são acima da média 1 ou não 0 e loc combinação de Sp e Sm indicando se o local é abaixo 1, médio 2 ou acima da média 3.

Conhecendo o banco de dados

glimpse(gamlss.data::rent)

Rows: 1,969
Columns: 9
$ R   <dbl> 693.3, 422.0, 736.6, 732.2, 1295.1, 1195.9, 395.0, 1285.6, 1175.4,…
$ Fl  <dbl> 50, 54, 70, 50, 55, 59, 46, 94, 93, 65, 36, 33, 57, 76, 75, 83, 83…
$ A   <dbl> 1972, 1972, 1972, 1972, 1893, 1893, 1957, 1972, 1972, 1967, 1972, …
$ Sp  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ Sm  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ B   <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ H   <fct> 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ L   <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, …
$ loc <fct> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …

Visualizando algumas linhas do banco de dados

gamlss.data::rent %>%
  slice_head(n=10) %>%
  gt::gt()

R	Fl	A	H	loc
693.3	50	1972	0	2
422.0	54	1972	0	2
736.6	70	1972	0	2
732.2	50	1972	0	2
1295.1	55	1893	0	2
1195.9	59	1893	0	2
395.0	46	1957	1	2
1285.6	94	1972	0	2
1175.4	93	1972	0	2
860.0	65	1967	0	2

Estatística Descritiva

skimr::skim(gamlss.data::rent)

Data summary
Name	gamlss.data::rent
Number of rows	1969
Number of columns	9
_______________________
Column type frequency:
factor	4
numeric	5
________________________
Group variables	None

Variable type: factor

skim_variable	complete_rate	ordered	n_unique	top_counts
B	1	FALSE	2	0: 1925, 1: 44
H	1	FALSE	2	0: 1580, 1: 389
L	1	FALSE	2	0: 1808, 1: 161
loc	1	FALSE	3	2: 1247, 3: 550, 1: 172

Variable type: numeric

skim_variable	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
R	1	811.88	379.00	101.7	544.2	737.8	1022	3000	▇▇▂▁▁
Fl	1	67.73	20.86	30.0	52.0	67.0	82	120	▆▇▇▅▂
A	1	1948.48	29.02	1890.0	1934.0	1957.0	1972	1988	▃▁▃▇▇
Sp	1	0.28	0.45	0.0	0.0	0.0	1	1	▇▁▁▁▃
Sm	1	0.09	0.28	0.0	0.0	0.0	0	1	▇▁▁▁▁

Letra A

Realize o ajuste de um modelo GAMLSS com a variável resposta R tendo uma distribuição Gama.

Verificando através de um histograma se a melhor distribuição é a Gama para a variável resposta

ggplot2::ggplot(gamlss.data::rent, ggplot2::aes(R))+
  ggplot2::geom_histogram(fill='white', 
                          color = "black", 
                          breaks = hist(gamlss.data::rent$R, plot = F)$breaks) +
  ggplot2::xlab("R")+
  ggplot2::ylab("Frequência")+
  ggplot2::labs(title = "Histograma da Variável R", subtitle = "")

Podemos observar que por se tratar de uma variável que está definida nos reais positivos, podemos utilizar a família Gama.

DescTools::Desc(gamlss.data::rent$R, plotit = F, digits = 3)

------------------------------------------------------------------------------ 
gamlss.data::rent$R (numeric)

     length        n      NAs   unique         0s       mean     meanCI'
      1'969    1'969        0    1'762          0    811.880    795.130
              100.0%     0.0%                0.0%               828.631
                                                                       
        .05      .10      .25   median        .75        .90        .95
    321.000  394.980  544.200  737.800  1'022.000  1'315.540  1'525.480
                                                                       
      range       sd    vcoef      mad        IQR       skew       kurt
  2'898.300  379.003    0.467  330.323    477.800      1.068      1.902
                                                                       
lowest : 101.700, 125.500, 127.100 (2), 130.700, 132.100
highest: 2'445.200, 2'537.900, 2'764.500, 2'869.900, 3'000.000

' 95%-CI (classic)

Modelando

Seleção de variáveis para mu

fit1_mu <- gamlss::gamlss(R~1, sigma.formula = ~1, data = gamlss.data::rent, family = GA)

GAMLSS-RS iteration 1: Global Deviance = 28611.58 
GAMLSS-RS iteration 2: Global Deviance = 28611.58

step_mu <- gamlss::stepGAIC(fit1_mu, scope = list(lower = ~1, upper = ~ Fl+A+Sp+Sm+B+H+L+loc))

Distribution parameter:  mu 
Start:  AIC= 28615.58 
 R ~ 1 

       Df   AIC
+ Fl    1 28104
+ H     1 28405
+ L     1 28526
+ B     1 28544
+ loc   2 28545
+ Sp    1 28570
+ A     1 28576
+ Sm    1 28578
<none>    28616

Step:  AIC= 28104.43 
 R ~ Fl 

       Df   AIC
+ H     1 27859
+ A     1 27996
+ loc   2 28026
+ B     1 28044
+ Sm    1 28050
+ Sp    1 28066
+ L     1 28071
<none>    28104
- Fl    1 28616

Step:  AIC= 27859.09 
 R ~ Fl + H 

       Df   AIC
+ loc   2 27799
+ Sm    1 27812
+ L     1 27835
+ B     1 27835
+ A     1 27836
+ Sp    1 27836
<none>    27859
- H     1 28104
- Fl    1 28405

Step:  AIC= 27799.17 
 R ~ Fl + H + loc 

       Df   AIC
+ B     1 27772
+ L     1 27777
+ A     1 27779
<none>    27799
- loc   2 27859
- H     1 28026
- Fl    1 28354

Step:  AIC= 27771.72 
 R ~ Fl + H + loc + B 

       Df   AIC
+ L     1 27750
+ A     1 27757
<none>    27772
- B     1 27799
- loc   2 27835
- H     1 27962
- Fl    1 28316

Step:  AIC= 27750.26 
 R ~ Fl + H + loc + B + L 

       Df   AIC
+ A     1 27736
<none>    27750
- L     1 27772
- B     1 27777
- loc   2 27811
- H     1 27934
- Fl    1 28248

Step:  AIC= 27735.96 
 R ~ Fl + H + loc + B + L + A 

       Df   AIC
<none>    27736
- A     1 27750
- L     1 27757
- B     1 27758
- loc   2 27794
- H     1 27863
- Fl    1 28250

Podemos observar, que para \mu as variáveis selecionadas foram Fl, H, loc, B, L e A.

Seleção de variáveis para sigma

fit2_sigma <- gamlss::gamlss(R~Fl+A+Sp+Sm+B+H+L+loc, data = gamlss.data::rent, family = GA)

GAMLSS-RS iteration 1: Global Deviance = 27717.96 
GAMLSS-RS iteration 2: Global Deviance = 27717.96

step_sigma <- gamlss::stepGAIC(fit2_sigma, what = "sigma", scope = list(lower= ~ 1, upper = ~ Fl+A+Sp+Sm+B+H+L+loc))

Distribution parameter:  sigma 
Start:  AIC= 27735.96 
 ~1 

       Df   AIC
+ A     1 27704
+ H     1 27726
+ loc   2 27730
+ Sp    1 27731
+ Sm    1 27733
+ Fl    1 27734
<none>    27736
+ B     1 27737
+ L     1 27738

Step:  AIC= 27703.93 
 ~A 

       Df   AIC
+ loc   2 27700
+ Sm    1 27701
+ Sp    1 27701
+ H     1 27703
<none>    27704
+ Fl    1 27705
+ B     1 27706
+ L     1 27706
- A     1 27736

Step:  AIC= 27699.56 
 ~A + loc 

       Df   AIC
<none>    27700
+ H     1 27700
+ Fl    1 27700
+ L     1 27702
+ B     1 27702
- loc   2 27704
- A     1 27730

Resposta: Depois de realizado o procedimento de seleção de variáveis para \sigma, obtivemos como resultado de que apenas as variáveis A e loc devem ser utilizadas.

Ajuste do modelo

fit11 <- gamlss::gamlss(R~Fl+H+loc+B+L+A, sigma.formula = ~ A + loc, data = gamlss.data::rent, family = GA)

GAMLSS-RS iteration 1: Global Deviance = 27677.96 
GAMLSS-RS iteration 2: Global Deviance = 27675.57 
GAMLSS-RS iteration 3: Global Deviance = 27675.56 
GAMLSS-RS iteration 4: Global Deviance = 27675.56

summary(fit11)

Warning in summary.gamlss(fit11): summary: vcov has failed, option qr is used instead

******************************************************************
Family:  c("GA", "Gamma") 

Call:  
gamlss::gamlss(formula = R ~ Fl + H + loc + B + L + A, sigma.formula = ~A +  
    loc, family = GA, data = gamlss.data::rent) 

Fitting method: RS() 

------------------------------------------------------------------
Mu link function:  log
Mu Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.3549442  0.6832552   3.447  0.00058 ***
Fl           0.0103633  0.0004077  25.421  < 2e-16 ***
H1          -0.2692798  0.0244981 -10.992  < 2e-16 ***
loc2         0.2082311  0.0331883   6.274 4.31e-10 ***
loc3         0.2738120  0.0348701   7.852 6.67e-15 ***
B1          -0.2789793  0.0652614  -4.275 2.00e-05 ***
L1           0.1360224  0.0299905   4.536 6.09e-06 ***
A            0.0017667  0.0003463   5.101 3.70e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

------------------------------------------------------------------
Sigma link function:  log
Sigma Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  5.2401707  1.0436990   5.021 5.61e-07 ***
A           -0.0031467  0.0005364  -5.866 5.22e-09 ***
loc2        -0.1049001  0.0559875  -1.874  0.06113 .  
loc3        -0.1676184  0.0603189  -2.779  0.00551 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

------------------------------------------------------------------
No. of observations in the fit:  1969 
Degrees of Freedom for the fit:  12
      Residual Deg. of Freedom:  1957 
                      at cycle:  4 
 
Global Deviance:     27675.56 
            AIC:     27699.56 
            SBC:     27766.58 
******************************************************************

Letra B

Faça uma análise de diagnósticos do modelo escolhido. O que você pode concluir do modelo?

plot(fit11)

******************************************************************
          Summary of the Quantile Residuals
                           mean   =  0.001228464 
                       variance   =  1.001448 
               coef. of skewness  =  -0.1341275 
               coef. of kurtosis  =  3.19966 
Filliben correlation coefficient  =  0.9989839 
******************************************************************

Resposta: De acordo com a análise gráfica podemos verificar que o modelo ficou bem ajustado e que, os resíduos aparentam estar distribuídos normalmente.

Questão 7

Seja o banco de dados Chwirut2 da biblioteca NISTnls. As variáveis são: y: valor resposta ultrassônicas e x: valor para a distância do metal. A função de relação não linear de y em x é dada por:

f(x,(\beta_1, \beta_2, \beta_3)) = \frac{\exp (-\beta_1x)}{\beta_2 + \beta_3}

Conhecendo o banco de dados

glimpse(NISTnls::Chwirut2)

Rows: 54
Columns: 2
$ y <dbl> 92.9000, 57.1000, 31.0500, 11.5875, 8.0250, 63.6000, 21.4000, 14.250…
$ x <dbl> 0.500, 1.000, 1.750, 3.750, 5.750, 0.875, 2.250, 3.250, 5.250, 0.750…

Visualizando algumas linhas do banco de dados

NISTnls::Chwirut2 %>%
  slice_head(n=10) %>%
  gt::gt()

y	x
92.9000	0.500
57.1000	1.000
31.0500	1.750
11.5875	3.750
8.0250	5.750
63.6000	0.875
21.4000	2.250
14.2500	3.250
8.4750	5.250
63.8000	0.750

Estatística Descritiva

skimr::skim(NISTnls::Chwirut2)

Data summary
Name	NISTnls::Chwirut2
Number of rows	54
Number of columns	2
_______________________
Column type frequency:
numeric	2
________________________
Group variables	None

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
y	0	1	31.56	26.31	3.75	11.56	20.10	58.68	92.9	▇▂▁▂▂
x	0	1	2.55	1.70	0.50	0.78	2.38	3.75	6.0	▇▅▅▂▃

Letra A

Construa o gráfico de dispersão das variáveis.

ggplot2::ggplot(NISTnls::Chwirut2, ggplot2::aes(x = x, y = y)) +
  ggplot2::geom_point()+
  ggplot2::xlab("Distância do metal") + ggplot2::ylab("Resposta Ultrassônica") + 
  ggplot2::labs(title = "Gráfico de dispersão da resposta ultrassônica com relação a distância do metal", subtitle = "")

Letra B

Utilizando o chute inicial: \beta = (0.1,0.01,0.1) qual a estimativa dos parâmetros do modelo? (acrescente a curva estimada no gráfico do item a).

expFct <- function(x, beta1, beta2,
 beta3) {
 exp(-beta1 * x)/(beta2 + beta3 *
 x)
 }

fit12 <- nls(y ~ expFct(x,
 beta1, beta2, beta3), data = NISTnls::Chwirut2,start = list(beta1 = 0.1, beta2 = 0.01,beta3 = 0.1))

summary(fit12)


Formula: y ~ expFct(x, beta1, beta2, beta3)

Parameters:
       Estimate Std. Error t value Pr(>|t|)    
beta1 0.1665779  0.0383034   4.349 6.56e-05 ***
beta2 0.0051653  0.0006662   7.753 3.54e-10 ***
beta3 0.0121500  0.0015304   7.939 1.81e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.172 on 51 degrees of freedom

Number of iterations to convergence: 6 
Achieved convergence tolerance: 4.429e-06

plot(y ~ x, data = NISTnls::Chwirut2, ylim = c(0, 100), ylab = "y", xlab = "x")
concVal <- with(NISTnls::Chwirut2, seq(min(x), max(x), length.out = 100))
lines(concVal, predict(fit12, newdata = data.frame(x = concVal)))

Letra C

Verifique a validação de f e das hipóteses de homocedasticidade e normalidade dos erros. Qual sua conclusão sobre o modelo?

plot(y~fitted(fit12), data = NISTnls::Chwirut2)

Homocedasticidade

e<-nlstools::nlsResiduals(fit12)
plot(e$resi2[,2]~fitted(fit12))

Resposta: Podemos ver que os dados estão bem distribuidos no intervalo, tendo uma fraca evidência de possível violação da homocedasticidade.

Normalidade dos Resíduos

qqnorm(e$resi2[,2])  
abline(a=0,b=1,col="red")

Questão 8

Seja o banco de dados USPop da biblioteca carData. As variáveis são: y: população (em milhões) dos Estados Unidos da América e x: ano do censo da população de 1790 a 2000. A função de relação não linear de y em x é dada por:

f(x,(\theta_1, \theta_2, \theta_3)) = \frac{\theta_1}{1+\exp(-(\theta_2 + x\theta_3))}

Conhecendo o banco de dados

glimpse(carData::USPop)

Rows: 22
Columns: 2
$ year       <int> 1790, 1800, 1810, 1820, 1830, 1840, 1850, 1860, 1870, 1880,…
$ population <dbl> 3.929214, 5.308483, 7.239881, 9.638453, 12.860702, 17.06335…

Visualizando algumas linhas do banco de dados

carData::USPop %>%
  slice_head(n=10) %>%
  gt::gt()

year	population
1790	3.929214
1800	5.308483
1810	7.239881
1820	9.638453
1830	12.860702
1840	17.063353
1850	23.191876
1860	31.443321
1870	38.558371
1880	50.189209

Estatística Descritiva

skimr::skim(carData::USPop)

Data summary
Name	carData::USPop
Number of rows	22
Number of columns	2
_______________________
Column type frequency:
numeric	2
________________________
Group variables	None

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
year	0	1	1895.00	64.94	1790.00	1842.5	1895.0	1947.50	2000.00	▇▆▆▆▇
population	0	1	94.68	87.26	3.93	18.6	69.6	146.54	281.42	▇▃▂▂▂

Letra A

Construa o gráfico de dispersão das variáveis.

ggplot2::ggplot(carData::USPop, ggplot2::aes(x = year, y = population)) +
  ggplot2::geom_point()+
  ggplot2::xlab("Ano do Censo") + ggplot2::ylab("População em Milhões") + 
  ggplot2::labs(title = "Gráfico de dispersão entre o tamanho da população com os anos do censo", subtitle = "Estados Unidos da América 1790-2000")

Letra B

Utilizando o chute inicial: \theta = (400,-49,0.025) qual a estimativa dos parâmetros do modelo? (acrescente a curva estimada no gráfico do item a).

fit13 <- nls(population ~ theta1/(1 + exp(-(theta2 + theta3*year))), start=list(theta1 = 400, theta2 = -49, theta3 = 0.025),data=carData::USPop, trace=TRUE)

3060.786    (2.47e+00): par = (400 -49 0.025)
558.5357    (4.69e-01): par = (426.062 -42.30786 0.02142146)
457.9746    (1.97e-02): par = (438.4147 -42.8369 0.02167713)
457.8071    (1.85e-03): par = (440.8903 -42.69866 0.02160152)
457.8056    (1.50e-04): par = (440.8168 -42.70805 0.02160649)
457.8056    (1.39e-05): par = (440.8345 -42.70688 0.02160586)
457.8056    (1.24e-06): par = (440.8333 -42.70698 0.02160591)

summary(fit13)


Formula: population ~ theta1/(1 + exp(-(theta2 + theta3 * year)))

Parameters:
         Estimate Std. Error t value Pr(>|t|)    
theta1 440.833344  35.000138   12.60 1.14e-10 ***
theta2 -42.706977   1.839138  -23.22 2.08e-15 ***
theta3   0.021606   0.001007   21.45 8.87e-15 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.909 on 19 degrees of freedom

Number of iterations to convergence: 6 
Achieved convergence tolerance: 1.239e-06

Letra C

Verifique a validação de f e das hipóteses de homocedasticidade e normalidade dos erros. Qual sua conclusão sobre o modelo?

plot(population~fitted(fit13), data = carData::USPop)

Homocedasticidade

e1<-nlstools::nlsResiduals(fit13)
plot(e1$resi2[,2]~fitted(fit13))

Resposta: De acordo com a análise gráfica podemos observar a violação da homocedasticidade.

Normalidade dos Resíduos

qqnorm(e1$resi2[,2])  
abline(a=0,b=1,col="red")

Resposta: Conforme é visto de acordo com a visualização gráfica, é observado que os resíduos não estão distribuidos normalmente.