Av Erik och Tobbe
Det finns underlag (NIDA) för att 50% av risken att bli beroende av t.ex. tobak kommer av ens genetiska arv. Vår fråga är om samma genetiska disposition också ökar risken för lungcancer. Vårt experiment går ut på att skilja tvillingar vid födseln och att tvinga den ena att röka och den andra att inte röka(*). Detta för att utesluta faktorer såsom uppfostran och miljö resp. kontrollera för rökningens inverkan på risken för lungcancer.
Vår hypotes är att risken för lungcancer är densamma hos båda tvillingarna och därför oberoende av rökning. Vi har fått anslag från ett flertal stora tobaksbolag och våra resurser är således i allt väsentligt obegränsade. Med detta i åtanke satsar vi främst på att samla ett stort underlag av tvillingpar från runt om i världen för optimerad randomisering. Varje tvillingpar är en experimental unit och således är den enda variabeln som är kopplad till dessa units huruvida båda får lungcancer – en indicator variable som är vår response variable med värdet 1 om det är fallet eller 0 om det inte är fallet.
Mätningen sker då alla försökspersoner gått ur tiden. Ett signifikant resultat skulle vara att andelen par där båda haft lungcancer var högre än andelen lungcancerfall generellt i populationen, i vilket fall vi då har anledning att tro att arvsanlaget spelar en viktig roll för lungcancerrisk.
(*) Vi är mycket förvånade över att etiska nämnden godkände detta experiment!
require(mosaic)
require(parallel)
# Funktion för att använda ordinala värden som numeriska
as.numeric.factor <- function(x) {
(as.numeric(levels(x))[x])
}
# Antal försökspersoner
antal = 250
# Sannolikhetsfördelningar för huruvida man dör av lungcancer som rökare
# resp. icke-rökare
icke = rbinom(n = antal, size = 1, prob = 0.05)
rökare = rbinom(n = antal, size = 1, prob = 0.15)
# Resamplar skillnaden i dödsrisken mellan grupperna
skillnad <- do(500) * {
mean(resample(rökare)) - mean(resample(icke))
}
confint(skillnad, method = "quantile")
## name lower upper level method
## 1 result 0.06 0.156 0.95 quantile
Konfidensintervallets lägre gräns tycks hamna någonstans mellan 0.01 och 0.08 vid upprepade försök när man använder 250 försökspersoner. Färre än så ger ett konfidensintervall som täcker nollan ibland.
Vi är intresserade av att se hur sambandet mellan snabbmatskonsumtion och depression ser ut. Ger ökat snabbmatsintag depression eller vice versa? För att ta reda på detta mäter vi en persons inkomst och snabbmatsutgifter i kronor samt ordinala skalor för stress och depression. Detta genomförs genom utvecklingen av en enkel app som kopplas till personens bankkonto. Appen mäter inkomster och relevanta utgifter. Förutom detta så ger appen en antal frågor relaterade till depression och stress några gånger i veckan. Datan från appen utgör dataunderlaget för experimentet. Varje individ är en experimental unit.
I vårt kausala diagram finns det direkta korrelationer mellan inkomst och depression; stress och depression respektive snabbmatskonsumtion och depression. Förutom dessa finns “back-door pathways” från såval inkomst som stress via snabbmatskonsumtion. Detta betyder att alla variabler måste vara med i en modell eftersom vi vill isolera snabbmatskonsumtionslänken.
# Genomsnittslön i samhället (kr per månad)
medLön = 20000
# Standardavvikelsen för snabbmatskonsumtion mätt i kr per månad
snabbmatSD = 25
# En data.frame som innehåller alla försökspersoners data
personer = do(100) * {
# En normalfördelning av potentiell inkomst för en enskild person
löneFördelning = rnorm(100, mean = medLön, sd = 5000)
# Personens faktiska inkomst plockat ur fördelningen
personensInkomst = deal(löneFördelning, 1)
# Stressnivå från 1 till 7, antas vara normalfödelad
personensStress = round(rnorm(1, mean = 3.5, sd = 1), 0)
# Genomsnittlig snabbmatsutgift baserat på personens inkomst och stress.
# Normaliserad gentemot personens teoretiska maxlön för att undvika att
# talet blir negativt (pga distributionen vi plockar ur).
personensSnabbmatMean = (160 + snabbmatSD - (160 * (personensInkomst/max(löneFördelning)))) *
1.2 * personensStress
# Faktisk snabbmatsutgift plockat ur en normalfödelning baserat på ovan
personensSnabbmat = rnorm(1, mean = personensSnabbmatMean, sd = snabbmatSD)
# Personens deppression som funktion av stress, snabbmatskonsumtion och
# inkomst. (Blir ett högt, abstrakt, värde pga lönen. Normaliseras nedan.)
personensDeppression = (personensStress * personensSnabbmat * personensInkomst)
# En person i data.frame-format
person = data.frame(inkomst = personensInkomst, stress = personensStress,
snabbmat = personensSnabbmat, depression = personensDeppression)
}
# Normaliserar depressionen genom att dela in den i 7 lika stora intervall
# och ger dessa namnen 1 till 7.
personer$depression = cut(personer$depression, 7, labels = c(1:7))
# Skriver ut vår experimentdata
personer
## inkomst stress snabbmat depression .row .index
## 1 25197 4 279.26 3 1 1
## 2 23528 3 302.60 2 1 2
## 3 13280 5 726.45 4 1 3
## 4 25454 4 313.02 3 1 4
## 5 15132 4 573.16 3 1 5
## 6 19392 3 333.41 2 1 6
## 7 21534 2 219.62 1 1 7
## 8 19375 4 424.48 3 1 8
## 9 23021 5 454.95 5 1 9
## 10 21934 5 410.99 4 1 10
## 11 19092 3 310.29 2 1 11
## 12 16570 4 482.12 3 1 12
## 13 12885 3 419.14 2 1 13
## 14 13809 4 546.17 3 1 14
## 15 21874 4 361.34 3 1 15
## 16 16788 3 370.81 2 1 16
## 17 20325 3 306.91 2 1 17
## 18 16745 3 433.36 2 1 18
## 19 28132 5 268.44 3 1 19
## 20 17028 2 224.54 1 1 20
## 21 22158 3 229.17 2 1 21
## 22 18709 2 191.83 1 1 22
## 23 21003 3 283.44 2 1 23
## 24 14716 2 260.40 1 1 24
## 25 21788 4 422.18 3 1 25
## 26 15747 3 388.12 2 1 26
## 27 20748 6 693.42 7 1 27
## 28 20403 3 297.73 2 1 28
## 29 18087 4 490.15 3 1 29
## 30 31242 4 230.97 3 1 30
## 31 13600 3 416.96 2 1 31
## 32 23642 4 362.55 3 1 32
## 33 15555 3 413.32 2 1 33
## 34 24203 5 413.33 4 1 34
## 35 21974 3 260.23 2 1 35
## 36 20430 2 219.21 1 1 36
## 37 21371 3 288.67 2 1 37
## 38 15435 4 559.62 3 1 38
## 39 22476 3 218.29 2 1 39
## 40 20147 3 273.94 2 1 40
## 41 18420 5 590.75 5 1 41
## 42 18050 5 544.99 4 1 42
## 43 19988 3 267.62 2 1 43
## 44 16193 4 482.68 3 1 44
## 45 26449 3 200.98 2 1 45
## 46 26402 4 251.24 3 1 46
## 47 16174 5 653.69 5 1 47
## 48 21996 3 266.33 2 1 48
## 49 23163 5 450.85 5 1 49
## 50 25402 4 316.73 3 1 50
## 51 21500 4 368.98 3 1 51
## 52 17788 4 508.34 3 1 52
## 53 13718 2 292.59 1 1 53
## 54 19134 3 265.58 2 1 54
## 55 14044 4 596.60 3 1 55
## 56 16790 2 256.80 1 1 56
## 57 23301 3 242.83 2 1 57
## 58 21602 1 76.50 1 1 58
## 59 16822 5 574.47 4 1 59
## 60 20174 5 495.89 4 1 60
## 61 17743 4 472.17 3 1 61
## 62 18092 4 553.25 4 1 62
## 63 26100 3 228.32 2 1 63
## 64 25982 4 253.38 3 1 64
## 65 17356 5 636.38 5 1 65
## 66 18031 4 448.26 3 1 66
## 67 18351 2 244.84 1 1 67
## 68 26096 2 108.76 1 1 68
## 69 13520 3 419.73 2 1 69
## 70 20117 3 266.82 2 1 70
## 71 12154 3 449.92 2 1 71
## 72 18707 4 395.02 3 1 72
## 73 18936 4 399.50 3 1 73
## 74 11041 2 343.39 1 1 74
## 75 25385 4 237.09 2 1 75
## 76 13506 5 690.07 4 1 76
## 77 6397 5 957.18 3 1 77
## 78 25530 6 486.45 7 1 78
## 79 16769 3 338.44 2 1 79
## 80 21319 3 192.70 1 1 80
## 81 16814 3 385.73 2 1 81
## 82 22306 4 392.38 3 1 82
## 83 21239 3 271.90 2 1 83
## 84 21478 3 308.89 2 1 84
## 85 23037 3 239.02 2 1 85
## 86 25108 3 202.63 2 1 86
## 87 21365 3 258.28 2 1 87
## 88 22017 4 368.93 3 1 88
## 89 19949 4 414.42 3 1 89
## 90 31439 4 213.90 3 1 90
## 91 29961 2 94.66 1 1 91
## 92 10589 3 480.40 2 1 92
## 93 21744 4 368.68 3 1 93
## 94 27339 4 130.69 2 1 94
## 95 26913 6 284.08 4 1 95
## 96 27739 5 210.30 3 1 96
## 97 12845 3 389.20 2 1 97
## 98 19057 3 363.33 2 1 98
## 99 18061 6 706.75 7 1 99
## 100 20800 2 190.30 1 1 100
xyplot(depression ~ snabbmat, data = personer, type = c("p", "r"))
bwplot(depression ~ snabbmat, data = personer)
model = lm(as.numeric.factor(depression) ~ snabbmat, data = personer)
confint(model, method = "quantile")
## 2.5 % 97.5 %
## (Intercept) 0.111362 1.076624
## snabbmat 0.004379 0.006802
I detta experiment vill vi testa hur fysisk ansträning påverkar arbetsminnet. Försökspersonerna fick först utföra ett fitness-test där deras allmänna fitness mättes som återhämtning från maxpuls till vilopuls, samt ett IQ-test. (i) Varje person fick sedan – utvilade – utföra ett test av arbetsminnet som gick ut på att memorera så många ord som möjligt. (ii) De fick också utföra en krävande fysisk övning på cross-trainer i 10 minuter varpå de testades med ett liknande arbetsminnestest. Ordningen på (i) och (ii) varierades slumpmässigt från person till person.
Vår response variable är resultatet på arbetsminnestestet, våra explanatory variables är IQ (kvantitativ), fitness (kvantitativ) samt fysisk aktivitet under experimentet (indicator variable). Eftersom tränings-varabeln varieras hos samma försöksperson är detta ett within-group experiment. Denna experimentdesign gör att vi kan utesluta IQ som varibel eftersom det är samma person som utför båda testen. Vår hypotes är att personer med hög fitness påvisar större prestationsökning efter fysisk aktivitet på arbetsminnestestet än de med låg fitness.
# 100 för lågt
antalPersoner = 150
personer = do(antalPersoner) * {
# Återhämtning efter träning i sekunder (tid från max- till vilopuls)
personensFitness = rnorm(1, mean = 180, sd = 60)
# En poissonfördelning över hur många ord en person kommer ihåg ur en lista.
# +1 för att ingen kommer ihåg 0 ord.
personensArbetsminne = rpois(n = antalPersoner, lambda = 4) + 1
# Plocka ett värde ur fördelningen för personens arbetsminne utan
# ansträgning
personensArbetsminneUtan = round(deal(personensArbetsminne, 1), 0)
# Normalfördelning över personens arbetsminne utan ansträgning modulerat av
# personens fitness. Plockar ett värde ur denna och rundar till heltal.
personensArbetsminneMed = round(rnorm(1, mean = (personensArbetsminneUtan +
1.5 * (100/personensFitness)), sd = 1), 0)
person = data.frame(fitness = personensFitness, arbetsminneUtan = personensArbetsminneUtan,
arbetsminneMed = personensArbetsminneMed)
}
personer
## fitness arbetsminneUtan arbetsminneMed .row .index
## 1 252.02 5 6 1 1
## 2 189.90 6 6 1 2
## 3 192.20 3 3 1 3
## 4 142.44 6 9 1 4
## 5 221.01 4 6 1 5
## 6 175.87 5 5 1 6
## 7 169.49 7 7 1 7
## 8 207.34 6 6 1 8
## 9 210.58 5 7 1 9
## 10 179.90 3 2 1 10
## 11 190.33 6 6 1 11
## 12 262.52 4 6 1 12
## 13 74.69 3 4 1 13
## 14 168.59 7 8 1 14
## 15 164.67 6 7 1 15
## 16 249.91 4 7 1 16
## 17 208.92 3 4 1 17
## 18 202.14 7 8 1 18
## 19 150.55 3 3 1 19
## 20 124.75 2 3 1 20
## 21 151.43 6 7 1 21
## 22 190.66 7 8 1 22
## 23 216.76 5 7 1 23
## 24 165.69 5 5 1 24
## 25 236.57 4 3 1 25
## 26 213.46 2 3 1 26
## 27 168.76 5 5 1 27
## 28 224.22 3 2 1 28
## 29 104.17 6 7 1 29
## 30 143.07 6 7 1 30
## 31 256.64 6 5 1 31
## 32 256.66 6 6 1 32
## 33 156.56 4 3 1 33
## 34 192.96 8 9 1 34
## 35 182.62 5 6 1 35
## 36 240.53 4 4 1 36
## 37 146.25 2 3 1 37
## 38 165.13 4 4 1 38
## 39 161.27 5 6 1 39
## 40 116.22 5 6 1 40
## 41 181.38 3 4 1 41
## 42 169.09 4 7 1 42
## 43 144.06 3 3 1 43
## 44 20.41 7 14 1 44
## 45 196.97 7 6 1 45
## 46 311.32 7 8 1 46
## 47 246.28 5 6 1 47
## 48 115.98 4 7 1 48
## 49 204.19 8 8 1 49
## 50 128.31 8 8 1 50
## 51 213.82 5 7 1 51
## 52 228.98 6 8 1 52
## 53 161.88 4 5 1 53
## 54 210.14 6 5 1 54
## 55 272.97 4 5 1 55
## 56 68.20 5 7 1 56
## 57 174.33 4 5 1 57
## 58 191.69 5 7 1 58
## 59 204.68 6 8 1 59
## 60 187.97 10 12 1 60
## 61 234.67 6 7 1 61
## 62 217.61 7 9 1 62
## 63 276.91 2 3 1 63
## 64 239.87 7 8 1 64
## 65 156.48 6 7 1 65
## 66 222.38 6 8 1 66
## 67 81.59 3 6 1 67
## 68 207.06 5 6 1 68
## 69 129.70 7 10 1 69
## 70 71.32 4 7 1 70
## 71 108.35 4 6 1 71
## 72 209.05 5 6 1 72
## 73 263.07 9 10 1 73
## 74 278.08 5 5 1 74
## 75 158.39 7 9 1 75
## 76 229.68 7 7 1 76
## 77 161.17 5 7 1 77
## 78 208.55 5 5 1 78
## 79 273.73 6 7 1 79
## 80 239.91 5 6 1 80
## 81 158.28 9 10 1 81
## 82 181.09 7 8 1 82
## 83 136.82 6 5 1 83
## 84 241.22 3 3 1 84
## 85 141.19 3 4 1 85
## 86 181.09 4 4 1 86
## 87 90.86 8 8 1 87
## 88 152.90 3 5 1 88
## 89 128.75 8 9 1 89
## 90 223.24 8 7 1 90
## 91 188.06 6 6 1 91
## 92 105.59 11 10 1 92
## 93 218.79 5 5 1 93
## 94 211.97 6 6 1 94
## 95 172.67 5 5 1 95
## 96 135.65 4 5 1 96
## 97 111.44 6 7 1 97
## 98 174.54 6 8 1 98
## 99 220.98 7 6 1 99
## 100 102.00 2 2 1 100
## 101 219.65 5 4 1 101
## 102 177.05 3 3 1 102
## 103 183.68 5 7 1 103
## 104 204.61 6 5 1 104
## 105 81.66 6 9 1 105
## 106 253.69 3 4 1 106
## 107 270.91 9 8 1 107
## 108 153.24 9 10 1 108
## 109 235.32 5 5 1 109
## 110 161.79 5 5 1 110
## 111 154.86 3 4 1 111
## 112 255.79 5 6 1 112
## 113 309.64 3 3 1 113
## 114 160.97 5 8 1 114
## 115 201.92 5 5 1 115
## 116 220.02 4 5 1 116
## 117 161.90 6 7 1 117
## 118 141.88 3 4 1 118
## 119 186.46 3 3 1 119
## 120 75.55 6 8 1 120
## 121 170.06 6 7 1 121
## 122 239.49 4 4 1 122
## 123 185.17 4 5 1 123
## 124 188.38 7 9 1 124
## 125 174.99 5 5 1 125
## 126 147.91 7 11 1 126
## 127 241.50 2 2 1 127
## 128 205.12 8 9 1 128
## 129 169.78 5 7 1 129
## 130 272.82 3 4 1 130
## 131 345.92 4 4 1 131
## 132 200.13 3 3 1 132
## 133 181.98 5 7 1 133
## 134 148.06 4 6 1 134
## 135 192.51 7 8 1 135
## 136 163.70 3 2 1 136
## 137 240.03 6 7 1 137
## 138 93.56 4 6 1 138
## 139 102.25 5 7 1 139
## 140 143.59 5 6 1 140
## 141 79.31 3 7 1 141
## 142 188.53 6 8 1 142
## 143 162.35 5 6 1 143
## 144 174.13 5 6 1 144
## 145 255.81 1 0 1 145
## 146 148.04 5 6 1 146
## 147 208.77 4 4 1 147
## 148 63.09 3 6 1 148
## 149 171.91 3 3 1 149
## 150 149.45 3 2 1 150
xyplot(arbetsminneMed - arbetsminneUtan ~ fitness, data = personer, type = c("p",
"r"))
bwplot(arbetsminneMed - arbetsminneUtan ~ fitness, data = personer)
histogram(~arbetsminneMed - arbetsminneUtan, data = personer)
model = lm(arbetsminneMed - arbetsminneUtan ~ fitness, data = personer)
confint(model, method = "quantile")
## 2.5 % 97.5 %
## (Intercept) 1.85693 3.137837
## fitness -0.01205 -0.005361
För att få ett konfidensintervall som inte täcker nollan måste man använda ca 150 försökspersoner. Kanske lite jobbigt.