연구 표본수 산출 방법(G*power test)

연구 표본수 산출 방법(G*power test)

표본이라는 용어를 이해하기 위해서는 모집단(population)과 표본집단(sample)의 특성을 알아야한다. 

모집단은 연구자가 연구하고자 하는 대상의 전체 집단을 이야기한다. 예를 들어 필자의 박사 학위 논문 제목인 '12주간의 클럽벨-스텝박스 순환운동이 남자 배드민턴 동호인의 견관절 기능, 체력, 배드민턴 기능 및 피로도에 미치는 영향' 이라는 제목을 놓고 보자. 모집단은 전세계에 있는 모든 남자 배드민턴 동호인이다. 모수(parameter)란 모집단에서 관심있는 측정값을 의미하므로 위에 나열된 종속변인(견관절 기능, 체력, 배드민턴 기능, 피로도)을 의미한다.

 

하지만 필자가 이 연구를 위해서 전세계에 있는 모든 남자 배드민턴 동호인을 데리고 연구를 할 수는 없지 않는가? 그래서 모집단의 특성을 가지고 있는 일부 인원들을 차출하여 연구를 진행하는 것이다. 이렇게 모집단의 특성을 가지고 차출된 인원을 표본 집단이라고 한다. 또한 이러한 표본 집단을 대상으로 연구를 진행하여 모집단에게 나올 결과를 추측하는 통계방법이 추리 통계(inferential statistics)라고 하며, 대부분의 인간 대상 연구에서는 이 기법을 사용한다. 

다양한 추리 통계 기법을 활용하여 표본에 대한 연구를 통해 완벽하지는 않지만 전체 모집단에 대하여 추론할 수 있다. 아래 표는 모집단의 모수와 표본의 통계치를 구분하기 위하여 일반적으로 사용하는 기호들이다. 연구에서 적절한 표본수를 구하기 위한 정보들을 공부하는 과정에 정리해서 남겨본다. 

계속해서 내용을 수정 및 업데이트 할 예정

대표적으로 추리 통계가 현실에서 사용되는 방법은 대통령 선거 이전에 조사되는 후보 지지율 조사이다. 이는 설문을 통한 표본을 사용하여 특정 후보자를 선호하는 사람들의 비율을 추정한다. 하지만 실제 투표를 했을 때와는 차이가 있을 수 있는 것 처럼 추리 통계는 약간의 오류가 있을 수 있다. 

이러한 표본의 추출 방법으로는 어떤 것이 있는지 알아보기 위해 Chat GPT, New bing 에게 도움을 청해봤다. 두 AI의 답변이 그리 다르지는 않다. Chat GPT는 대화에 특화된 모델이라 부드러운 대화를 위해서 거짓 정보를 스스로 만들어내 답변하기도 하기 때문에 개인적으로는 항상 New bing 으로 교차검증을 하여 정보를 얻고 있다. 하지만 확실히 아래 보이는 것 처럼 Chat GPT 가 더 친절한 답변을 주기는 한다. 

Chat GPT 의 답변

연구를 함에 있어서 적절한 표본수는 상당히 중요한데, 표본수가 너무 많으면 효과가 없음에도 불구하고 효과가 있다고 나오는 1종 오류가 발생할 확률이 높아지고, 반대로 표본수가 너무 적으면 효과가 있음에도 불구하고 효과가 없다고 나오는 2종 오류가 발생할 확률이 높아진다. 따라서 좋은 연구 결과를 위한 적절한 표본수가 어느 정도인지 선행연구를 통하여 검증하는 과정이 반드시 필요하다. 

예를 들어 위 표를 보면 상관관계 분석에서 r 값이 같아도 표본수가 증가하면 p 값이 유의하게 나올 수 있다.

일반적으로 실험 연구에서는 집단 당 15명, 횡단적 연구(상관분석)에서는 집단 당 30명 정도를 권장하지만 G*power test 를 하면 더 세부적으로 결정할 수 있고, 최근 IRB 에서는 표본수 산정 기준을 제시할 것을 요구하기 때문에 이는 아무쪼록 반드시 필요하다.

또한 이건 연구실마다 다르겠지만 p 값은 이와 같이 표본수에 영향을 많이 받기 때문에 95% 신뢰구간(95% CI)과 효과크기를 함께 제공하는게 좋다고 보는 곳도 있다. 95% 신뢰구간에 대해서는 다음 글에 정리하도록 하겠다. 효과크기는 이제 곧 설명이 나올거다.

표본수 산정 시 유의사항

과학적 방법(scientific method)은 지식을 습득하기 위해 추리 통계를 사용하는 것이며, 과학적 가설의 전개와 그 가설에 대한 추리 통계를 바탕으로 가설을 검정하는 것이다. 

*연구가설(H1
, 대립가설, alternative hypothesis): 연구하고자 하는 가설
*귀무가설(H0, null hypothesis): 꽝!(무로 돌아감.)

예를 들어, 위 그림이 토끼인지 병아리인지 연구를 통해서 증명하고자 한다고 가정해보자. 연구자들은 이 그림을 토끼라고 생각한다. 그렇다면 연구가설과 귀무가설은 다음과 같다.

* 연구가설(H1): 토끼 -> 증명 해야하는 것!
* 귀무가설(H0): 병아리

이와 같이 이 그림이 토끼인지 병아리인지 판단하는 과정을 가설 검정이라고 한다. 이러한 가설검정은 귀무가설을 기각하거나 채택하는 것으로 이루어진다. 쉽게 말해서 귀무가설을 채택하면 위 그림은 병아리가 되는 것이고, 귀무가설을 기각하면 위 그림은 토끼가 되는 것이다. 

모든 가설 검정은 차이가 없다는 귀무가설을 기준으로 시작하고, 통계를 통해 차이가 있다는 것을 검증하는 것이라고 보면 된다. 

가설 검정에는 단측 검정(one-sided test, H1 : P1>P2)과 양측 검정(two-sided test, H1 : P1≠P2)이 있는데, 단측 검정은 연구가설이 특정 방향으로 설정되어 있는 경우에 사용된다. 예를 들어, 우리가 어떤 새로운 치료법이 기존 치료법보다 더 효과적이라고 주장하는 경우, 귀무가설은 "새로운 치료법이 기존 치료법보다 더 효과적이지 않다"와 같은 형태로 설정된다. 이 경우, 우리는 새로운 치료법이 더 효과적이라는 방향으로 검증을 하게 되며, 따라서 단측 검정을 사용한다. 

반면에, 양측 검정은 연구가설이 "두 그룹 간 차이가 있다" 혹은 "두 그룹 간 차이가 없다"와 같이 양쪽 방향으로 설정되어 있는 경우에 사용된다. 예를 들어, 우리가 어떤 제품의 무게가 10 g이라고 주장하는 경우, 귀무 가설은 "제품의 무게는 10 g이 아니다"와 같이 양쪽 방향으로 설정된다. 이러한 경우 이 제품이 10 g 보다 높을수도 있고 낮을 수도 있기 때문에 양측 검정을 사용한다. 

따라서, 단측 검정과 양측 검정은 검증하고자 하는 연구가설이 어떻게 설정되어 있는지에 따라서 선택되어야 한다. 하지만 대부분의 스포츠의과학 쪽 연구에서는 양측 검정을 사용한다. 

G*power test 의 화면을 보면, Tail(s) 항목에 One 을 선택하면 단측 검정, Two 를 선택하면 양측 검정을 의미한다. 

1종 오류,  2종 오류

위 표에 나와있는 1종 오류, 2종 오류, 검정력에 대하여 간단하게 알아보자. 
간단한 예시를 들기 위해서 법에 대입해서 한번 알아보자.

1) 1종 오류(α): 무죄 -> 유죄. 무죄인 사람을 유죄로 잘못 구속할 확률.
2) 2종 오류(β): 유죄 -> 무죄. 유죄인 사람을 무죄로 잘못 풀어줄 확률.
3) 검정력(power, 1-β): 유죄 -> 유죄. 유죄인 사람을 유죄로 구속할 확률.

1종 오류와 2종 오류 둘다 모두 문제가 있지만 1종 오류가 더 심각한 문제이다. 법에서는 1종 오류를 막기 위해 무죄추정의 원칙이라는게 있다.

1종 오류와 2종 오류는 함께 움직인다고 볼 수 있다. 1종 오류가 증가되면 2종 오류가 감소한다. 극단적인 예시를 들어보자. 경찰서에 온 사람을 다 구속시켜버리면 아무도 풀려날 일이 없기 때문에 2종 오류가 없어진다(1종 오류↑ - 2종 오류↓).

반대로 2종 오류가 증가하면 1종 오류가 감소한다. 역시나 극단적인 예시를 들어보자. 경찰서에 온 사람을 다 풀어줘버리면 아무도 구속될일이 없기 때문에 1종 오류가 없어진다(2종 오류↑ - 1종 오류↓)

이러한 두 상황 모두 발생하면 안되기 때문에 1종 오류와 2종 오류를 동시에 적당히 줄여주는 것이 중요하다. 이를 위해 상대적으로 더 위험한 1종 오류를 1차적으로 예방하기 위하여 1종 오류의 최대 허용 한계를 부여하고, 그것을 넘지 않도록 한다. 이것을 유의 수준
(significant level)이라고 한다. 1종 오류과 유의 수준 모두 α 라는 기호로 사용되고 있지만 명백히 말하면 이 둘은 다른 뜻 이므로 구분할수 있어야한다. 일반적으로 대부분의 연구에서는 유의수준을 5%로 정하고 있으며, 이를 기준으로 잡고 재주껏 2종 오류를 줄이기 위해 노력해봐야한다. 이것이 가설 검정의 원칙이다. 

앞서 설명한 내용을 운동의 효과를 입증하기 위한 연구라고 가정하여 적어봤다. 

1) 1종 오류(α): 운동효과가 없는데 운동효과가 있다고 판단할 확률
2) 2종 오류(β): 운동효과가 있는데 운동효과가 없다고 판단할 확률
3) 검정력(1-β): 운동효과가 있는데 운동효과가 있다고 판단할 확률

일반적으로 연구자가 연구가설을 성공적으로 입증하려면, 1종 오류는 낮아야하고 검정력은 높아야한다. 검정력은 통상적으로 80% 이상으로 설정하며, 90% 이상으로 설정하기도 한다.

검정력을 높게 설정하면 요구되는 표본수가 증가하는 것을 볼 수 있다. 검정력이 높아진다는 것은 2종 오류가 감소한다는 뜻이고, 2종 오류가 감소하려면 더 많은 표본수가 필요하다. 따라서 높은 검정력을 시행하려면 더 많은 표본수가 필요한 것이다. 

1종 오류와 2종 오류를 줄일 수 있는 방법은?

방금 앞서서도 살짝 설명했지만 조금 더 자세히 알아보자.

1종 오류의 해결책은 유의 수준을 낮추는 것이다. 유의 수준을 5% 에서 1%로 낮추는 것이 해결책이 될 수 있다. 이에 대한 예시로 1종 오류는 집단 혹은 시기가 여러개 일때, ANOVA 를 돌리지 않고, t-test 를 반복적으로 돌렸을 때 발생할 수 있다. 예를 들어 A, B, C의 집단 간 비교를 할 때, 일반적으로는 세 집단부터는 one-way ANOVA를 통해 집단간 비교를 하고, 통계적으로 유의한 차이가 있을 때 사후 검증을 통해 더 세부적인 분석을 한다.

 

하지만 두 집단간 비교를 하는 독립 t-test 를 사용한다고 하면 A vs B, B vs C, C vs A 총 3번의 독립 t-test를 돌려야하는 것이다. 세 집단까지는 이렇게 독립 t-test 를 돌려도 나름 괜찮은데 네 집단이 넘어가면서 부터는 이렇게 반복적으로 독립 ttest를 돌리게되면 유의 수준을 5%로 설정해도 1종 오류가 발생할 가능성이 높다. 따라서 유의 수준을 1%로 낮추어서 돌리면 이것을 예방할 수 있는 것이다. 하지만 굳이 이렇게 번거롭게 할 필요는 없으니 ANOVA를 돌리도록 하자.

2종 오류의 해결책은 표본수를 늘리는 것이다. 즉, 우리가 G*power test 를 하는 이유도 적절한 표본수를 설정하여 2종 오류를 최소화 하기 위함도 있다. 이것을 위해 검정력을 80% 혹은 90%로 설정한다. 

효과크기(effect size)

효과 크기란 두 집단간의 평균 차이의 정도를 의미하는 것이라고 생각하면 된다. 효과크기에는 여러 가지가 존재한다. Raw mean 
difference(D), standardized mean difference(d or g), response ratios(R), risk ratio(RR), odds ratio(OR), risk difference(RD), correlation(r) 등등.. 이 외에도 f, Ω 등 다양한 종류의 효과크기가 존재한다. 이것들은 모두 상호 교환이 가능하다. 어떤 통계 기법을 돌리느냐에 따라서 효과크기도 달라진다. 우리 연구실에서 자주 사용하는 ANOVA는 f 를 사용한다. 

효과크기는 실험 연구에서는 처치 효과가 얼마나 강력한지를 나타내는 지표가 될 수 있다. 즉, 독립변인이 종속변인을 얼마나 강하게 변화시켰는지 그 정도를 나타내는 것이다. p값은 처치 효과의 유무(효과가 있는지 없는지)를 나타내는 것이라면, 효과크기는 처치 효과의 정도(얼마나 많이 변화시키는지)를 나타내는 것이라고 생각하면 된다.


효과 크기의 여러가지 유형을 간단하게 정리하면 다음과 같다. 

연속형 데이터에서 대중적으로 많이 사용되는 것이 표준화된 평균 차이(standardized mean difference:SMD)인 cohen's d 이다. 이는 평균과 표준편차가 있으면 구할 수 있으며, 흔히 메타분석에서도 많이 사용되는 효과크기이다. RCT를 바탕으로 한 메타분석에서는 사후 검사 결과의 효과크기를 사용하며, 실험연구에서도 두 집단을 비교한 t-검증에서는 cohen's d를 사용한다. 

 

범주형 변수에 해당되는 효과 크기로는 OR, RR, RD 등이 있다. 이는 보통 특정 위험인자가 질병의 발생 유무에 미치는 영향에 대한 연구에서 주로 사용되는 보건 의료쪽 연구에서 많이 등장한다.

효과 크기를 구하는 방법은 다양하지만 d 값은 위와 같은 표처럼 계산을 진행한다. "집단 1의 평균 - 집단 2의 평균/ 평균 표준편차(pooled standard deviation)". 우리가 계산 공식을 모두 알아야하는 것은 아니다. 대상자 수, 평균, 표준편차를 알고 있으면 아래 사이트에서 효과 크기를 계산할 수 있다. 

사이트에 들어가보면 연속형 데이터에 대한 효과 크기 구하는 방법이 standardized mean difference(SMD)라는 메뉴에 다양하게 나와있다.

또 다른 연속형 데이터인 상관분석에 대한 효과 크기를 구하는 방법도 제시되어 있다.

범주형 데이터에 대한 효과크기 구하는 방법도 제시가 되어있다.
엑셀로도 계산할 수 있으니 아래 링크를 참고하자. 


엑셀에서 효과크기(Cohen's D) 구하기
두 데이터 세트, 또는 두 그룹 간에 차이가 있는지 없는지 확인하…
loadtoexcelmaster.tistory.com

G*power test 시 효과 크기를 정하는 방법

연구를 계획할 때 내가 보려고 하는 주요 종속변인(primary outcome)이 선행연구의 유사한 처치에 의해서 얼마만큼의 효과크기를 보이며 변화하였는지를 확인하고, 그에 맞는 표본수를 모집해야 좋은 연구 결과를 도출할 수 있다. cohen's d 같은 경우에는 대상자 수, 평균, 표준편차만 있으면 쉽게 구할 수 있지만 f 같은 경우에는 더 많은 데이터가 필요해 계산이 어려울 수 있다. 이때는 효과 크기가 서로 상호 교환이 된다는 사실을 기억하고 cohen's d를 사용해 f 를 구하면 된다. 효과 크기를 서로 변환해주는 바로 아래 사이트에서!

아래 링크는 다양한 효과 크기가 정리되어있는 사이트! 어느 정도 숫자가 각 효과 크기의 small, medium, large 인지 확인할 수 있다. 

위 내용을 좀 더 자세히 한번 정리해보겠다. G*power test 할 때 효과 크기를 정하는 방법은..

1. 과거의 비슷한 연구를 참고
관련된 연구가 많이 있을 경우 메타 분석을 사용해 추정하는 것이 가장 바람직하다.

 

2. 관련 연구가 거의 없다면, 작은 규모의 예비 실험(pilot study)를 실시해 추정

 

3. 예비 실험을 실시할 수 없는 경우, 가상표를 작성해 임상적/이론적 가치가 있기에 충분히 크다고 간주할 수 있는 가장 작은 효과크기 계산
-예) 어떤 금연 요법이 적어도 10%(small)의 금연 효과를 가져올 때, 효과적이고 가치 있는 것이라고 간주할 수 있다면, 효과크기를 이 값에 근거하여 추정

 

4. 앞의 세가지 방법을 적용할 수 없을 경우, 과거의 경험적 결과로부터 효과크기가 작은지(small), 적당한지(medium), 큰지(large)만을 선택해 그에 따라 정해진 추정치를 사용 기존에 많이 이루어지지 않던 새로운 연구 주제에서는 보통 효과크기가 기대보다 작을 수 있다. 따라서 이러한 경우 처치 효과를 얻어내려면 더 많은 표본수가 필요한 것이다.

 

최근 효과 크기에 대한 관심 증가하면서 많은 분야의 학계에서 연구자들이 연구 결과를 보고할 때, 효과크기를 함께 보고해야한다고 주장한다. APA publication manual 에서는 양적 방법론을 사용했음에도 효과크기를 보고하지 않은 연구는 결점이 있다고 간주하기도 한다. 

표본수 구하기

현대시대 다양한 분야에서 연구자가 주장하려는 가설 입증은 표본을 추출하여 선행 연구결과와 현재 연구결과가 차이가 있는지 또는 대조군과 처리군이 다른지 통계분석을 실시한다.

이때 표본을 어느 정도 추출해야 통계적 신뢰도를 보장받는지가 주된 관심이다. 따라서 이 장은 연구자의 신뢰성과 타당성 보장을 위한 표본수 계산 방법을 소개한다.

표본수를 계산하는 프로그램은

• r project → 통계프로그램(무료)

• g power → 통계프로그램(무료)

• SPSS SamplePower → 통계프로그램(유료)

등이 있다.

통계적 가설검정

통계적 가설 검정 계산에 사용하는 몇 가지 용어를 소개한다.

• 귀무가설 (null hypothesis;H0 ) : 이미 전에 연구자가 입증한 가설
• 대립가설 (alternative hypothesis;H1 ) : 연구자가 연구결과 입증하려는 가설로 귀무가설이 아닌 가설
• 유의수준 (significance level;�) : 연구자가 귀무가설이 옳은데도 잘못하여 귀무가설을 기각하는 오류의 최대 허용한계의 확률로 많은 경우 0.05로 설정
• 기각역 (rejection region for H0 ) : 유의수준에서 귀무가설을 기각하는 영역
• 유의확률 (probability value) : 데이터에서 구한 검정통계량이 귀무가설이 옳음에도 잘못하여 귀무가설을 기각하는 오류의 확률
• 검정통계량 (test statistic) : 귀무가설 기각, 채택여부 판정할 때 사용하는 통계량으로 관측한 데이터에서 계산
• 제 1종의 오류 확률 (type I error probability; �) : 유의수준 과 동일한 의미
• 제 2종의 오류 확률 (type II error probability; �) : 귀무가설 이 옳지 않을 때 귀무가설을 기각하지 않는 오류 확률
• 검정력 (power;1−�) : 귀무가설이 옳지 않을 때 귀무가설을 기각하는 오류의 확률

다음 그림은 귀무가설 기각역, 채택역, 제 1종의 오류, 제 2종의 오류, 귀무가설, 대립가설을 설명한다.

단일표본 평균 검정 사례

초콜렛 제품의 무게는 평균 260g, 표준편차가 10g으로 알려져 있다. 이 초콜렛은 두 개의 공장에서 만든다고 하자. 각 공장마다 20개씩 표본을 추출한 결과, A 공장에서 만든 초콜렛의 평균이 262.9g 이었고, B 공장은 264.7g 이었다. 각 공장에서 만든 초콜렛은 평균이 260g이라고 할 수 있는가? 아니면 260g보다 크다고 할 수 있는가?

문제에서 주어진 모수와 통계량으로 계산한 결과가 아래 그림에 있다.

여기서 초콜렛의 무게가 260g 보다 크다고 할 수 있는 값, 즉 통계적으로 유의하다고 판단할 수 있는 기각역은 과자의 평균 무게가 264.7g보다 크면 되겠다. 따라서 A 공장은 통계적으로 유의하지 않고, B 공장은 통계적으로 유의하다. 즉 A 공장의 초콜렛 평균 262.9는 예전 평균 260g 보다 크다고 할 수 없고, B 공장의 초콜렛 평균 264.7은 예전 평균 260g 보다 크다고 할 수 있다.

단일표본 평균에 대한 표본수

엑셀에서 정규분포 분위수 계산은

NORMSINV(probability)

함수이며, 입력값은

• probability : 확률

이다. 양측검정에서 유의수준 �=0.05로 설정하였으면 �1−�/2=�0.975이므로 다음과 같이 입력한다.

=NORMSINV(0.975) 
1.95996

t–분포 분위수 계산은

TINV(probability, d.f )

사례1: 단일 표본 평균의 표본수

=TINV(A9,H9-1)

1. 데이터 → 데이터 도구 → 가상 분석 → 목표값 찾기 메뉴 선택
2. 목표값 찾기 창에

• 수식 셀에는 |왼쪽 n − 오른쪽 식|인 H9 클릭한다.
• 찾는 값에는 두 식의 차이 0을 입력한다.
• 값을 바꿀 셀에는 구할 표본수 n인 G9를 클릭한다.

1. 목표값 찾기로 구한 표본수이다.

두 번째는 통계자료분석 프로그램인 R로 표본수를 계산하였다.

> power.t.test(sd=5.1,delta=3.3,type="one.sample",
alternative="two.side",power=.8,sig.level=0.05)
One-sample t test power calculation
n = 20.74904
delta = 3.3
sd = 5.1
sig.level = 0.05
power = 0.8
alternative = two.sided

독립인 두 표본 평균 차이에 대한 표본수 - 집단 크기가 같은 경우

사례2: 독립인 두 표본 평균 차이의 표본수

 

 

\

1. 데이터 → 데이터 도구 → 가상 분석 → 목표값 찾기 메뉴 선택
2. 목표값 찾기 창에
3. 목표값 찾기로 구한 표본수이다.

계산결과 각 집단별 176명 이상 표본을 추출한다.

다음은 통계자료분석 프로그램인 R로 계산한 결과이다.

> power.t.test(sd=10,delta=3,type="two.sample",
alternative="two.side",power=.8,sig.level=0.05)
Two-sample t test power calculation
n = 175.3851
delta = 3
sd = 10
sig.level = 0.05
power = 0.8
alternative = two.sided
NOTE: n is number in *each* group

독립인 두 표본 평균 차이에 대한 표본수 - 집단크기가 다른 경우

사례3: 독립인 두 표본 평균 차이의 표본수

• 
• 데이터 → 데이터 도구 → 가상 분석 → 목표값 찾기 메뉴 선택
• 목표값 찾기 창에
• 목표값 찾기로 구한 표본수이다.

식 (7)에 대하여 100번 반복 계산한 결과 만족하는 해를 찾지 못했지만 근사적으로 대조군 118명, 처리군 118×3=354 이상 표본을 추출한다.

단일표본 비율에 대한 표본수

독립인 두 표본 비율 차이에 대한 표본수 - 집단크기가 같은 경우

사례5: 독립인 두 표본 비율 차이의 표본수

R 프로그램으로 표본수를 계산한 결과

> power.prop.test(p1=0.56,p2=0.30,power=0.8,alternative = "two.sided")

     Two-sample comparison of proportions power calculation 

              n = 55.72232
             p1 = 0.56
             p2 = 0.3
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

NOTE: n is number in *each* group

독립인 두 표본 비율 차이에 대한 표본수 - 집단크기가 다른 경우

사례6: 독립인 두 표본 비율 차이의 표본수

여차 저차해서 표본수를 구했는데.. 중도 탈락율 20%를 고려해야한다고 했을때 흔히 잘못 계산하는 오류는 다음과 같다.

* 중도 탈락율을 고려하는데 100명이 필요한데 중도 탈락율이 20% 예상된다면?

100*1.2=120명(x), 100/(1-0.2)=125명(o)

표본수 선정에 대한 기본 개념 영상 강의

 

 

A7NOVA G0 *Power test 돌리는 방법

 

 

효과 크기 계산 하는 방법

 

 

케어무브 수원점 장소의 다른 글

#논문통계

 

반응형

 

SMALL