Advanced linear models and classification 2024 Assessed computer practical 1
PLEASE READ CAREFULLY
This computer practical counts 25% towards your final mark on this unit. It is due on 7th March 2024, by 3pm and should be submitted on Blackboard. You should submit your computer practical answers as a single pdf file containing your knitted Markdown file and handwritten answers if applicable – please name the file after your first name and surname. You should number your answers, both in the Markdown file and your handwritten notes, if applicable.
FAILURE TO SUBMIT YOUR ASSESSED COMPUTER PRACTICAL IN THE REQUIRED FORMAT WILL RESULT IN IT NOT BEING MARKED.
The computer practical is concerned with simple, yet practically very im- portant, variations on the models covered in the course so far and you are taken step by step through the ideas.
1 One-way ANOVA (Analysis of Variance)
The one way ANOVA model corresponds to multilinear regression in the situa- tion where the covariates are categorical. Consider the following example where the effect of different fertilizers, A,B,C and D on crop yields was tested. For each fertilizer 3 independent yields are recorded (in, say, tons)
Fertilizer | Crop yields |
---|---|
A | 33.63|37.80|36.58 |
B | 35.83|38.18|37.89 |
C | 42.95|40.43|41.46 |
D | 38.02|39.62|35.99 |
Table 1: Crop yields and fertilizers
The type of fertilizer is said to be categorical as it is not a meaningful numerical value measuring a specific quantity. We would like to be able to predict the effect of fertilizers on yields using the observations (Fertilizer, Yield) given in the table, with Yield as a response and Fertilizer as covariate. Ultimately we would like to determine which fertilizer(s) is (are) best. To that effect we would like to use what we know about linear regression.
- What is the apparent difficulty here?
- A possible solution could consist of replacing A, B, C and D with numer-ical values e.g. make the substitutions A ← 1,B ← 2,C ← 3 and D ← 4.
- (a) Do you think this is a sensible idea? Why?
- (b) One can test the idea numerically on this example.
Yield <- c(33.63, 37.8, 36.58, 35.83, 38.18, 37.89, 42.95,
40.43, 41.46, 38.02, 39.62, 35.99)
x1 <- rep(1:4, c(3, 3, 3, 3))
Explain what the second line is doing and what the (statistical) aim is? Fit a linear model to the data using this substitution and predict the observed yields using the fitted model.
(c) Fit a linear model with now the substitution A ← 3,B ← 1,C ← 4 and D ← 2 and again predict the observed yields. What do you observe and what is your conclusion?
- ANOVA model
(a) Explain why using the following linear model for the yields,
for , a residual of expectation 0 and F an abbreviation of Fertilizer could be a good choice? You should justify your answer by considering the model for each row in Table 1 and provide an interpretation for the coefficients
.
(b) Can you foresee a problem with the model as defined above when trying to estimate ?This problem is called iden- tifiability. Suggest a way to address this issue and interpret your choice.
(c) Provide an interpretation of the problem above in terms of the matrix X from the lecture notes.
- Let’s investigate how the function lm can be used to perform linear regres- sion with factors (i.e. in the situation where some or all of the covariates are categorical).
(a) The function lm handles factors automatically, but we try to under- stand how, with a bit of experimentation. Explain what the aim of the following command, which uses the function gl, is? You should explain the reason for the particular choice of parameters used in connection to yield and the values contained in factor.x.
factor.x <- gl(4, 3, 12, labels = c("A", "B", "C", "D"))
(b) You have already encountered the contrasts function in the third Practice Computer Practical.
contrasts(factor.x)
## BCD
## A 0 0 0
## B 1 0 0
## C 0 1 0
## D 0 0 1
In the light of the earlier discussion concerned with identifiability, what choice does R make? You should explain your answer by refer- ring to the output of the constrasts function above.
(c) Fit the ANOVA model to (Yield,F). Check that lm treats in the way indicated by contrasts.
(d) Interpret the output of lm. What is your conclusion about the relative efficiency of the various fertilizers?
(e) Does an analysis of deviance confirm your conclusion?
2 Contingency table: Poisson sampling
Consider the longitudinal study of coronary heart disease. In Table 2, 1329 patients were cross-classified according to their level of cholesterol (below or above 260 in some unit) and the presence of a heart disease. It is of interest to determine if the two factors are linked or not.
Serum cholesterol | Heart disease Present Absent | Total |
---|---|---|
<260 >260 | 51 992 41 245 | 1043 286 |
Total | 92 1237 | 1329 |
Table 2: Cholesterol levels and heart conditions of patients
To model this problem, there are similarities with the ANOVA model, but also differences. Here the observations within each cell are assumed to be dis- tributed according to a Poisson distribution corresponding to the number of participants with presence (P) or absence (A) of a heart disease and with serum cholesterol larger than 260, (L), or smaller (S). The mean of the Poisson distribution within each cell is modelled for and as
where and model the marginal "main" dependence on C and H only, while models the interaction between the variables C and H (e.g. do C and H tend to occur simultaneously or not, in other words are the variables dependent?). The counts within each cells are assumed independent.
- In R you are given the observations
y <- c(51, 992, 41, 245)
(a) Following what you did for the ANOVA model, create the factor lists factor.C and factor.H corresponding to C and H suitable for use with y.
(b) Check the help for the function interaction and create the factor list factor.CH suitable for use with y.
Using the glm function fit a Poisson count model assuming that there is no interaction between C and H.
Repeat the above assuming now interaction between C and H.
Which model fits the data best? Confirm your observations by perform- ing a suitable deviance analysis using the anova function. What is the conclusion of this study?
Solution
Question 1
# 定义产量数据
Yield <- c(33.63, 37.8, 36.58, 35.83, 38.18, 37.89, 42.95, 40.43, 41.46, 38.02, 39.62, 35.99)
# 为化肥A, B, C, D分配数值1, 2, 3, 4
x1 <- rep(1:4, each=3)
# 拟合线性模型
model1 <- lm(Yield ~ x1)
# 查看模型摘要
summary(model1)
# 使用拟合的模型预测产量
predict(model1)
# 为化肥A, B, C, D重新分配数值3, 1, 4, 2
x2 <- rep(c(3, 1, 4, 2), each=3)
# 重新拟合线性模型
model2 <- lm(Yield ~ x2)
# 查看模型摘要
summary(model2)
# 使用拟合的模型预测产量
predict(model2)
> # 定义产量数据
> Yield <- c(33.63, 37.8, 36.58, 35.83, 38.18, 37.89, 42.95, 40.43, 41.46, 38.02, 39.62, 35.99)
> # 为化肥A, B, C, D分配数值1, 2, 3, 4
> x1 <- rep(1:4, each=3)
> # 拟合线性模型
> model1 <- lm(Yield ~ x1)
> # 查看模型摘要
> summary(model1)
Call:
lm(formula = Yield ~ x1)
Residuals:
Min 1Q Median 3Q Max
-3.698 -1.719 0.060 1.252 4.255
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 35.7150 1.7274 20.676 1.55e-09 ***
x1 0.9933 0.6308 1.575 0.146
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.443 on 10 degrees of freedom
Multiple R-squared: 0.1987, Adjusted R-squared: 0.1186
F-statistic: 2.48 on 1 and 10 DF, p-value: 0.1464
> # 使用拟合的模型预测产量
> predict(model1)
1 2 3 4 5 6 7 8 9
36.70833 36.70833 36.70833 37.70167 37.70167 37.70167 38.69500 38.69500 38.69500
10 11 12
39.68833 39.68833 39.68833
> # 为化肥A, B, C, D重新分配数值3, 1, 4, 2
> x2 <- rep(c(3, 1, 4, 2), each=3)
> # 重新拟合线性模型
> model2 <- lm(Yield ~ x2)
> # 查看模型摘要
> summary(model2)
Call:
lm(formula = Yield ~ x2)
Residuals:
Min 1Q Median 3Q Max
-5.1217 -1.1275 0.4733 1.6117 3.0917
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 35.4317 1.6749 21.154 1.24e-09 ***
x2 1.1067 0.6116 1.809 0.1
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.369 on 10 degrees of freedom
Multiple R-squared: 0.2467, Adjusted R-squared: 0.1713
F-statistic: 3.274 on 1 and 10 DF, p-value: 0.1005
> # 使用拟合的模型预测产量
> predict(model2)
1 2 3 4 5 6 7 8 9
38.75167 38.75167 38.75167 36.53833 36.53833 36.53833 39.85833 39.85833 39.85833
10 11 12
37.64500 37.64500 37.64500
1. 问题的困难所在
在题目中,肥料类型(A、B、C、D)是分类变量,不具有可量化的数值意义,这导致我们不能直接使用这些分类变量作为线性回归模型中的协变量。线性回归模型通常需要数值型数据作为输入,而分类数据的直接使用会导致模型无法正确理解和处理这些变量之间的关系。
2. 使用数值替换分类变量的可能解决方案
(a)这是否是一个明智的想法?
将分类变量(A、B、C、D)替换为数值(1、2、3、4)可能不是一个好的解决方案,因为这种替换隐含地引入了一个假设,即肥料之间存在着一种线性关系(例如,B 比 A 有更多的影响,C 比 B 有更多的影响,以此类推),这种假设在实际情况中可能并不成立。分类变量的数值替换可能会导致模型误解肥料之间的实际关系。
结论: 不是一个好的解决方案。这种做法隐含了不成立的线性关系假设,可能导致模型误解肥料之间的实际关系。
(b) 数值替换的测试
第二行代码创建了一个向量 x1
,它重复了 1 到 4 的序列,每个值重复 3 次,对应于每种肥料的三个产量观测值。这个向量用于在拟合线性模型时表示肥料类型的数值。但这种方法忽略了肥料类型之间实际的非数值性质。
结论: x1
向量的创建是为了在拟合线性模型时表示肥料类型的数值,但这种方法忽略了肥料类型之间实际的非数值性质。
在第一个线性模型中,代码使用了rep(1:4, each=3)
来为化肥 A、B、C、D 分别分配数值 1、2、3、4,并拟合了模型 Yield ~ x1
。输出显示,x1 的系数估计值为 0.9933,但其 p 值为 0.146,意味着在统计学上这个系数并不显著,我们不能拒绝零假设(即系数等于 0,表示化肥类型的变化对产量没有显著影响)。
该模型的 R 方值为 0.1987,说明模型解释了总变异的近 20%,但调整后的 R 方值降至 0.1186,提示模型的解释能力相对有限。这个结果表明,虽然模型试图通过化肥类型的数值替换来预测产量,但化肥类型作为一个序数变量对产量的解释力并不强。
使用该模型预测的产量显示了一种规律性增加,反映了数值编码的顺序性,但由于系数不显著,这种增加可能并不反映实际情况。
(c) 使用不同的数值替换进行线性模型拟合
使用不同的数值替换(A←3, B←1, C←4, D←2)进行线性模型拟合并预测产量可能会导致不同的结果,这进一步证明了将分类变量替换为数值可能会误导模型,因为不同的替换方案会导致不同的模型输出,这说明肥料类型之间的关系并非线性,但它不能准确反映分类变量的实际影响。
结论: 不同的数值替换导致不同的模型输出,说明将分类变量替换为数值可能会误导模型,因为它不能准确反映分类变量的实际影响。
在第二个模型中,化肥 A、B、C、D 被重新分配了数值 3、1、4、2(使用rep(c(3, 1, 4, 2), each=3)
),然后拟合了模型 Yield ~ x2
。这次 x2 的系数估计值为 1.1067,p 值稍微降低至 0.1,尽管这依旧在统计学上不显著,但显示了稍微强一些的趋势。
这一模型的 R 方值为 0.2467,相较于第一个模型有所提高,调整后 R 方值为 0.1713,显示了略微增强的解释能力。重新分配数值后的模型同样没有提供有力的证据支持化肥类型对产量有显著影响的假设。
预测产量的结果同样受到了新的数值编码的影响,展现了不同的规律性变化,这进一步说明了将分类变量简单替换为数值并进行线性拟合可能会导致误导性的解释,因为预测值的变化仅仅反映了数值编码的变化,而非化肥类型实际对产量的影响。
3. ANOVA 模型
(a) 使用线性模型的理由
使用 ANOVA 模型是因为它允许我们将分类变量(如肥料类型)作为模型的一部分,而不需要将其转换为数值。在这个模型中, 是基准产量, 分别代表不同肥料相对于基准产量的效果差异。这种模型能够更准确地反映不同肥料对产量的影响,因为它直接将肥料类型作为因素考虑进去,而不是试图通过数值转换来间接表示这些影响。
结论: ANOVA模型是更好的选择,因为它可以直接处理分类变量,如肥料类型,而无需将其转换为数值。这种模型能更准确地反映不同肥料对产量的影响。
(b) 模型的可识别性问题
模型定义中存在的一个问题是所有 系数可能不是唯一确定的,这是因为如果我们给所有的系数同时增加一个常数,模型预测的值不会改变,这就是所谓的不可识别性问题。为了解决这个问题,我们可以约束一个 系数(例如,将 设为 0),作为参考组,其他组的 系数代表与这个参考组的差异。这种方式可以确保模型的系数是可识别的,并且便于解释。
结论: 为了解决可识别性问题,可以通过设置一个系数(例如)为0作为参考点,其他系数表示相对于这个参考点的效果差异。这样可以确保模型的系数是可识别的。
(c) 矩阵 X 的解释问题
在题目中提到的矩阵 X 代表了设计矩阵,用于连接观测值和模型参数。在处理分类变量时,矩阵 X 包含了用于表示分类变量每个级别的指示变量(dummy variables)。如果没有适当处理可识别性问题,设计矩阵X可能会导致列线性相关,进而影响模型参数的估计。
结论: 设计矩阵X需要正确处理分类变量,以避免列线性相关问题,确保模型参数的可识别性和准确估计。
在题目中提到的矩阵 X 是设计矩阵,它包含了所有的预测变量(包括常数项和分类变量的指示变量)。可识别性问题在矩阵 X 的语境下意味着矩阵 X 的某些列线性相关,这会导致 X 的逆矩阵不存在,从而无法使用最小二乘法直接估计回归系数。
4. 使用 lm
函数进行线性回归和 ANOVA
(a) gl
函数的目的
gl
函数用于生成因子变量 factor.x
,表示肥料类型。参数选择与产量数据直接相关,确保每种肥料对应的观测值正确地被标记。这样可以在不需要手动转换分类变量为数值的情况下,直接在模型中使用肥料类型作为一个因素。
# 生成因子变量
factor.x <- gl(4, 3, 12, labels = c("A", "B", "C", "D"))
(b) R 中的默认对比选择
contrasts
函数输出展示了 R 在处理分类变量时的默认对比设置(即虚拟编码)。在这种编码下,第一个水平(A)被作为基准,其他水平(B、C、D)相对于基准水平的效果被编码为 0 或 1。这种方法解决了可识别性问题,因为它自动省略了一个基准类别,避免了全 1 的列。
# (b) 查看对比
contrasts(factor.x)
(c) 拟合 ANOVA 模型
使用 lm
函数拟合(Yield, F)的 ANOVA 模型时,R 会按照 contrasts
指示的方式处理分类变量。这确保了模型的拟合是按照预期的方式进行,考虑了分类变量的特性。
# (c) 拟合ANOVA模型
Yield <- c(33.63, 37.8, 36.58, 35.83, 38.18, 37.89, 42.95, 40.43, 41.46, 38.02, 39.62, 35.99)
lm.model <- lm(Yield ~ factor.x)
(d) lm
输出的解释
# (d) 解释输出
summary(lm.model)
# (e) 方差分析
anova(lm.model)
> # (a) 生成因子变量
> factor.x <- gl(4, 3, 12, labels = c("A", "B", "C", "D"))
> # (b) 查看对比
> contrasts(factor.x)
B C D
A 0 0 0
B 1 0 0
C 0 1 0
D 0 0 1
> # (c) 拟合ANOVA模型
> Yield <- c(33.63, 37.8, 36.58, 35.83, 38.18, 37.89, 42.95, 40.43, 41.46, 38.02, 39.62, 35.99)
> lm.model <- lm(Yield ~ factor.x)
> # (d) 解释输出
> summary(lm.model)
Call:
lm(formula = Yield ~ factor.x)
Residuals:
Min 1Q Median 3Q Max
-2.3733 -1.2550 0.3600 0.9942 1.7967
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 36.0033 0.9641 37.345 2.9e-10 ***
factor.xB 1.2967 1.3634 0.951 0.36942
factor.xC 5.6100 1.3634 4.115 0.00337 **
factor.xD 1.8733 1.3634 1.374 0.20670
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.67 on 8 degrees of freedom
Multiple R-squared: 0.7005, Adjusted R-squared: 0.5882
F-statistic: 6.237 on 3 and 8 DF, p-value: 0.01726
> # (e) 方差分析
> anova(lm.model)
Analysis of Variance Table
Response: Yield
Df Sum Sq Mean Sq F value Pr(>F)
factor.x 3 52.172 17.3907 6.237 0.01726 *
Residuals 8 22.306 2.7883
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
根据 lm
函数的输出,可以得到以下信息:
- (Intercept) 36.0033: 这是基准化肥(A)的平均产量估计值,标准误为0.9641,t值为37.345,非常显著地不同于0(几乎接近于 0的p值)。
- factor.xB 1.2967: 这表明B化肥相对于A化肥平均增加了1.2967吨的产量,但这个差异不是统计显著的(p=0.36942)。
- factor.xC 5.6100: C化肥相对于A化肥显著增加了5.61吨的产量(p=0.00337),这是一个显著的结果。
- factor.xD 1.8733: D化肥相对于A化肥增加了1.8733吨的产量,但这个差异同样不是统计显著的(p=0.20670)。
残差的标准误为1.67,反映了模型拟合数据的精度。多重决定系数(R-squared)为0.7005,说明模型解释了大约70.05%的产量变异。调整后的R-squared为0.5882,考虑到自由度调整后的解释力度略有下降。
(e) 方差分析
ANOVA 表显示:
- factor.x:化肥类型有3个自由度,总平方和为52.172,平均平方和为17.3907,F值为6.237,p值为0.01726,表明至少有一种化肥的效果与其他化肥显著不同。
- Residuals:残差有8个自由度,总平方和为22.306,平均平方和为2.7883。
方差分析的结果证实了化肥类型对产量有显著的影响,具体来说,C 化肥相比于 A 化肥有显著提高的效果。
Question 2
1. 创建因子列表
(a) 创建因子列表 factor.C
和 factor.H
首先,需要为 C
(胆固醇水平,低于或高于 260)和 H
(心脏病存在与否)创建因子列表。y
向量中的四个值对应于(低胆固醇&心脏病存在,低胆固醇&心脏病不存在,高胆固醇&心脏病存在,高胆固醇&心脏病不存在)的病例数。
factor.C <- factor(rep(c("S", "L"), each=2))
factor.H <- factor(rep(c("P", "A"), times=2))
(b) 使用 interaction
函数创建 factor.CH
interaction
函数可以用来创建两个或更多因子的交互作用的因子。对于 C
和 H
的交互,可以这样做:
factor.CH <- interaction(factor.C, factor.H)
2. 不考虑交互作用的泊松回归模型
使用 glm
函数拟合一个泊松回归模型,此时不考虑 C
和 H
之间的交互作用。模型中只包括 C
和 H
的主要效应:
model1 <- glm(y ~ factor.C + factor.H, family="poisson")
3. 考虑交互作用的泊松回归模型
现在,假设 C
和 H
之间存在交互作用,使用 glm
函数再次拟合模型,但这次包括交互项:
model2 <- glm(y ~ factor.C * factor.H, family="poisson")
4. 模型比较和结论
使用 anova
函数进行偏差分析,比较两个模型,以确定哪一个更适合数据:
anova(model1, model2, test="Chisq")
完整代码:
# 数据准备
y <- c(51, 992, 41, 245) # 观察值
# 创建因子列表
factor.C <- factor(rep(c("S", "L"), each=2)) # 胆固醇水平因子: S代表小于260, L代表大于或等于260
factor.H <- factor(rep(c("P", "A"), times=2)) # 心脏病存在与否因子: P代表存在, A代表不存在
# 交互作用因子
factor.CH <- interaction(factor.C, factor.H) # C和H的交互作用因子
# 拟合模型(不考虑交互作用)
model1 <- glm(y ~ factor.C + factor.H, family="poisson")
# 拟合模型(考虑交互作用)
model2 <- glm(y ~ factor.C * factor.H, family="poisson")
# 模型比较
anova_result <- anova(model1, model2, test="Chisq")
# 打印偏差分析结果
print(anova_result)
# 结论
# 通过查看anova_result的输出,你可以比较两个模型的偏差和卡方检验的P值。
# 如果模型2相对于模型1有显著更小的偏差,并且P值小于显著性水平(例如0.05),
# 那么可以认为包括C和H的交互作用可以更好地拟合数据,从而得出C和H是相关的结论。
> # 数据准备
> y <- c(51, 992, 41, 245) # 观察值
> # 创建因子列表
> factor.C <- factor(rep(c("S", "L"), each=2)) # 胆固醇水平因子: S代表小于260, L代表大于或等于260
> factor.H <- factor(rep(c("P", "A"), times=2)) # 心脏病存在与否因子: P代表存在, A代表不存在
> # 交互作用因子
> factor.CH <- interaction(factor.C, factor.H) # C和H的交互作用因子
> # 拟合模型(不考虑交互作用)
> model1 <- glm(y ~ factor.C + factor.H, family="poisson")
> # 拟合模型(考虑交互作用)
> model2 <- glm(y ~ factor.C * factor.H, family="poisson")
> # 模型比较
> anova_result <- anova(model1, model2, test="Chisq")
> # 打印偏差分析结果
> print(anova_result)
Analysis of Deviance Table
Model 1: y ~ factor.C + factor.H
Model 2: y ~ factor.C * factor.H
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 1 26.43
2 0 0.00 1 26.43 2.733e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
根据提供的偏差分析结果,我们可以判断哪个模型更适合数据。分析结果显示,当考虑胆固醇水平(C)和心脏病存在与否(H)之间的交互作用时(即模型2),残差偏差(Resid. Dev)从26.43减少到0.00,且与模型1相比,增加的自由度(Df)为1,偏差下降了26.43,这种改进是非常显著的(Pr(>Chi) = 2.733e-07,即p值<0.001,非常显著)。
最佳模型
模型2(考虑C和H的交互作用)明显比模型1(只考虑C和H的主效应)更适合数据。这是因为模型2的残差偏差显著减少到0,表明模型2能更准确地拟合观察到的数据。
偏差分析的确认
通过偏差分析(ANOVA),我们可以看到模型1与模型2相比,偏差显著减少,且p值远小于0.001,表明加入交互项后模型的拟合度显著提高。因此,偏差分析确认模型 2 更适合数据。
研究结论
研究的结论是,胆固醇水平与心脏病存在与否之间存在显著的交互作用。这意味着胆固醇水平对心脏病风险的影响依赖于是否存在心脏病,或者说心脏病的存在与否会改变胆固醇水平对心脏病风险的影响。因此,我们不能仅考虑这两个因素的独立效应,还必须考虑它们的交互作用,以更准确地了解心脏病的风险因素。
公众号:AI悦创【二维码】
AI悦创·编程一对一
AI悦创·推出辅导班啦,包括「Python 语言辅导班、C++ 辅导班、java 辅导班、算法/数据结构辅导班、少儿编程、pygame 游戏开发、Web、Linux」,全部都是一对一教学:一对一辅导 + 一对一答疑 + 布置作业 + 项目实践等。当然,还有线下线上摄影课程、Photoshop、Premiere 一对一教学、QQ、微信在线,随时响应!微信:Jiabcdefh
C++ 信息奥赛题解,长期更新!长期招收一对一中小学信息奥赛集训,莆田、厦门地区有机会线下上门,其他地区线上。微信:Jiabcdefh
方法一:QQ
方法二:微信:Jiabcdefh
- 0
- 0
- 0
- 0
- 0
- 0